Gear mesh stiffness

[cfernandes]

I did a static non linear simulation of a spur gear in order to evaluate the mesh stiffness. I am using CCX 2.17. To do that I created a rigid body to the hub of the pinion and the wheel. The simulation has 4 steps:

1. imposed rotation on the pinion to initialize the contact (keeping wheel fixed)

2. Impose a small torque on the wheel (keeping pinion fixed)

3. Impose the desired torque on the wheel (keeping pinion fixed)

4. Impose the desired resisting torque on the wheel and impose a rotation on the pinion corresponding to 2 mesh cycles.

A surface to surface penalty contact without friction is used, everything converges and the results are (apparently) similar to previous works. The problem relies on the results of Reaction Force (RF) for the rotation nodes used for the rigid body (step 4):

• the reference nodes don’t belong to the node set;
• the imposed torque is correct on the wheel;
• the imposed rotation is correct on the by pinion;
• but the torque evaluated on the pinion rotation node is not constant and also depends on the penalty coefficient. I expected to have influence on the rotation angles only…

It would be helpful to have your feedback. My experience with CalculiX is mostly doing thermal analysis of gears, so probably I am missing some important points.

Best regards,
Carlos Fernandes

[JuanP74]

not sure if applicable to you, but please be aware of this remark in the manual:

6.11.5 Forces obtained by selecting RF
This section has been included because the output when selecting RF on a *NODE PRINT, *NODE FILE or *NODE OUTPUT card is not always what the user expects. With RF you get the sum of all external forces in a node. The external forces can be viewed as the sum of the loading forces and the reaction forces.

[cfernandes]

Juan,

Thanks for the hint. I will verify in a different way. My main doubt relies on the change of the RF as the penalty coefficient increases.

[jbr]

@cfernandes

My main doubt relies on the change of the RF as the penalty coefficient increases.

I recall a topic on contact stiffness setting. Based on that discussion, I would expect the torque evaluated on the pinion rotation node not to be constant and will depend on the penalty coefficient. Are you using linear elastic materials?

[cfernandes]

Yes, I am using linear elastic material just to validate the model with the standards/literature.

The shape of the mesh stiffness curve is exactly as expected, but its value is quite dependent on the penalty coefficient and much lower than the standards (my penalty coefficient is high, if the penalty decreases, the value of mesh stiffness also decreases even more due to penetration. Penetration changes the relation between pinion and wheel angles). Maybe I have to compute the mesh stiffness in a different way…

Thank you for the recall of the previous topic.

[cfernandes]

Thank you all for the help.

I decided to ask a colleague to run my input deck on Abaqus. It runs perfectly just needed to change the reference nodes/rotation nodes. I decided to share my observations in the case that it can be useful…

I verified that the values at the reference nodes on the CalculiX are different just in a small part of the step. If I increase the penalty coefficient it deviates in a larger part of the step.

The zones where CalculiX doesn’t fit the expected behaviour is coincident with wrong torque evaluation.

The mesh stiffness is evaluated only considering the node rotations (the torque is not needed for this analysis, but is expected to be approximately constant). I used a smaller penalty on CalculiX just to have a better shape for the curve (with 10^14 it fits partially the Abaqus results, but for less points).

I am now wondering what is the reason for this…

[JuanP74]

Did you use quadratic elements? Remember that ccx linear elements do not have locking prevention features. From the manual 4 Golden Rules:
3. USE QUADRATIC ELEMENTS (C3D10, C3D15, C3D20(R), S8, CPE8,CPS8, CAX8, B32). The standard shape functions for quadratic elements are very good. Most finite element programs use these standard functions.For linear elements this is not the case: linear elements exhibit all kind of weird behavior such as shear locking and volumetric locking.

[cfernandes]

JuanP74,

thank you for the warning. The stresses were “as expected” and I couldn’t figure easily this possible weird behaviour with the first order elements. Maybe it is exactly the cause. I have to use a coarser mesh with 2nd order elements otherwise I need a week of computation. But I will check it for sure. Thanks.

[cfernandes]

I did a trial with 2nd order elements C3D10:

The behaviour is not better.

[cfernandes]

Hi,

I managed to get a good solution both with first order and 2nd order elements. I had to use much much larger (10 times bigger than in Abaqus) number of steps and increments to get the expected result.

Thank you all.

[Matej]

Great work. The results are really in good agreement and confirm the Calculix as a great alternative to commercial codes.

May I ask, how it comes to a 10x larger number of increments. What is the reason for that. Did you not use automatic incrementation?

[cfernandes]

Hi Matej,

Maybe I was not precise in my text. I had to divide the rotation of the gear in a much larger number of divisions. In reality, if the rotation tends to a very small value the result tends to the expected one. At the same time I was increasing the number of increments of STEP calculation, since it seems to help. However, since now I found the cause, I have room to explore and see what is really making the difference.

[tian]

hi! carlos. I use currently another model, but have the same problem as you: my displacement of elements is strongly affected by changing the contact coefficient. Do you find the reason why this happens? Thank you!

best regards,
Tian

[Calc_em]

Can you say more about your model ? If you provide the necessary details, we will try to help you find the cause of this problem.

[tian]

Below is my .inp file. I have an airfoil structure and try to bend it downwards using *CLOAD to create a moment with a value of 0.273. This moment is about the Rot point inside of a rigid body (defined as a surface). There is a sliding contact condition between the stingers and the skin outside. SO I used here the *Contact pair card. I am using C3D8I elements.

My problem is: If I change my contact coefficient from 1e7 to 1e6, and the stick slope from 1e5 to 1e4 (because the stick slope is 100 times smaller than the coefficient). The max. displacement of the airfoil tip will decrease from 1.08e-2 to 9.94e-3. Thank you a lot!

*CONTACT PAIR, INTERACTION= PCONT, TYPE=SURFACE TO SURFACE, ADJUST=0.01
BACtoSkin , skinBot
*CONTACT PAIR, INTERACTION= PCONT, TYPE=SURFACE TO SURFACE, ADJUST=0.01
AirfoiltoSkin , skinBot
*SURFACE INTERACTION, NAME = PCONT

*SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=tied
1.e7

*FRICTION
0.75, 1e5

**182745, 0.1683200514098, 0.0 , 3.521669270E-04
*node, nset=ref
3000000, 0.1683200514098, 0.005 , 3.521669270E-04
3000001, 0.1683200514098, 0, 2

*RIGID BODY, NSET=ROTATE, REF NODE=3000001, ROT NODE=3000000

*NSET, NSET=ROTNODE
3000000

*STEP, INC=10000
*STATIC

• CLOAD
  3000000, 2, 0.273

NODE FILE
U
EL FILE
S
NODE PRINT, NSET=DeltaU
U
NODE PRINT, NSET=SpineCurve
U
NODE PRINT, NSET=ROTNODE
RF, U
END STEP

[Calc_em]

The difference in results is not so big considering the fact that you change the contact parameters by the order of magnitude. I would rather think about the choice of contact formulation - why did you choose tied contact ? Do you have some reference data that makes you use contact parameters other than defaults ?

[tian]

Thanks for your help! I got a code from my supervisor. He used K as 1e9 and stick slope of 0.001 (which does not work well).
(I have different skin thicknesses of this airfoil structure and make them achieve the same Tip displacement. In the end, I compare different moments for different skin thicknesses. The moments should increase with thicker skin. But up to now, the moments do not show an increased tendency but first increase then decrease with thicker skin ). I just changed his parameters to try to get better results and to solve the problem above.

I used tied because this way converges most quickly. I also tried with linear pressure overclosure, but the calculation did not converge. Another reason for tied is that the skin must be attached on top of the stringers, otherwise, the skin will fall off.

my supervisor’s code:

*SURFACE INTERACTION, NAME=IntPROP-1
*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=tied
1.e9
*friction
0.1, 0.001

[Calc_em]

How are these parts connected in real life ? Have you tried using tie constraint (*TIE) ? Maybe such simplification will be fine for your case.

[tian]

In order to hold the skin on structure, it is under 10% pretension. I used the temperature to simulate this pretension (the skin shrinks first 10% of its length and then comes the step of the moment). Because the skin is not parallel to the global coordinate. I cannot use the *initial condition, type=stress. SO in real lift, the tensed skin just covers the stringers.

I have had the TIE condition in other interfaces. But between skin and stringers, I have to use sliding contact conditions to reduce the skin out-of-plane displacement (upwards or downwards), which is also the goal of my research. If I use TIE between them, it shows a big enough out-of-plane displacement under aerodynamic pressure, which causes an uneven surface of the airfoil. I want to solve this out-of-plane displacement problem.

---

**Pretension

*INITIAL CONDITIONS, TYPE=TEMPERATURE
SkinN, 0

*STEP, INC=10000
*STATIC
*TEMPERATURE
SkinN, -0.1

*END STEP

---