Units, trying to get correct results with Calculix

[kwahoo]

The units
Calculix does not use units itself, we have to use consistent unit system. I think, mm,N,s,K (Calculix 2.7 manual, page 18) should be the best choice:
Input:
Force [N]
Length [mm]
Young Modulus [MPa or N/mm^2]

Output:
Displacement [mm]
Stress [MPa or N/mm^2]

Meanwhile FreeCAD uses kPa in the material library, and forces are multiplied x1000 inside MechanicalAnalysis.py - I have no idea why. I changed, inside MechanicalAnalysis.py:

Code: Select all

        Force = (ForceObject.Force * 1000.0) / NbrForceNods

to

Code: Select all

        Force = ForceObject.Force / NbrForceNods

and

Code: Select all

        inpfile.write('{0:.3f}, '.format(YM.Value) )

to

Code: Select all

        inpfile.write('{0:.3f}, '.format(YM.Value/1000) ) 

Additionally I had to do small change inside Steel.FCMat, because of parser error

Code: Select all

YoungsModulus=2e+08 

to

Code: Select all

YoungsModulus=2.0e+08

Next, when you choose a face to add force, the displayed unit is N/mm^2. I think it may be incorrect, because it is not a pressure, but force, later divided by number of nodes. The unit should be just N. Am I right?


The problem
No matter if vanilla FreeCAD or FreeCAD with my changes, I cannot get correct stress and displacement results. Check the file https://app.box.com/s/14giy87qvm9c0291urq0
It is an I-beam HEB100, 1000mm length. One end fixed, other end 10 000N force along Y. The results (confirmed by Z88 FEM tool) should be:

Name Ix Iy Wx Wy ix iy
cm4 cm3 cm
100 HEB 450 167 90 33 4,2 2,53

Wx=90000mm^3
M=F*l=10000N*1000mm=10000000Nmm
Sig=M/W=10000000Nmm/90000mm^3=111N/mm^2=111MPa
f=F*l^3/(3*E*J)
f=10000N*(1000mm)^3/(3*206000N/mm^2*4.5*10^6mm^4)=3.6mm

FreeCAD + Calculix gives me results different by magnitudes.


Final question, why

Code: Select all

*INITIAL CONDITIONS, TYPE=STRESS, USER

is necessary in the Calculix *.inp file?

[dubstar-04]

Should the length units be in Meters? (N, M and N/m^2)

[kwahoo]

Any consistent unit system should be ok.

[ickby]

I have not looked into this, but in general in FEM a force [N] can only be applied to a vertex. If you apply a force to a edge it needs to be a lineforce [N/m] and this can then be applied to all nodes on the edge via the mesh element lengths. The same holds for forces applied to faces: They must be [N/m^2] to be well defined as you need to calculate the per mesh node force via the surface of the elements. So the freecad notation seems right in this regard.

[jriegel]

Hi,
I had no time so far to use the new unit parser on the FEM workbench or the material module. So at the moment all units have to be mm/kg/deg. That involves a good deal of calculation (you can use the unit calculator under the tools menu). In the long run you can use all supported units in the materials or the postprocessing.

[kwahoo]

The same holds for forces applied to faces: They must be [N/m^2] to be well defined as you need to calculate the per mesh node force via the surface of the elements.

So it have to be calculated inside Calculix. Here are two inp files (vanilla FreeCAD):
for 10x10x10 box (1N/mm^2(?) force on the 9 node 10x10 face):

Code: Select all

*Node , NSET=Nall
1,10,10,0
2,10,10,10
3,10,0,0
4,10,0,10
5,0,10,0
6,0,10,10
7,0,0,0
8,0,0,10
9,0,5,10
10,5,10,10
11,10,5,10
12,5,0,10
13,0,5,0
14,5,10,0
15,10,5,0
16,5,0,0
17,10,10,5
18,0,10,5
19,10,0,5
20,0,0,5
21,5,5,10
22,5,5,0
23,5,10,5
24,5,0,5
25,10,5,5
26,0,5,5
*Element, TYPE=C3D10, ELSET=Eall
25,7,4,6,8,24,21,26,20,12,9,
26,7,6,4,1,26,21,24,22,23,25,
27,5,7,1,6,13,22,14,18,26,23,
28,4,7,1,3,24,22,25,19,16,15,
29,4,2,1,6,11,17,25,21,10,23,


*NSET,NSET=FemConstraintFixed
3,
4,
7,
8,
12,
16,
19,
20,
24,


*NSET,NSET=FemConstraintForce
1,
2,
5,
6,
10,
14,
17,
18,
23,


*MATERIAL, Name=Steel
*ELASTIC 
200000000.000, 0.300
*SOLID SECTION, Elset=Eall, Material=Steel
*INITIAL CONDITIONS, TYPE=STRESS, USER
*STEP
*STATIC
*BOUNDARY
FemConstraintFixed,1,3,0.0
*CLOAD
FemConstraintForce,1,0.0
FemConstraintForce,2,111.11111111111111
FemConstraintForce,3,0.0
*NODE FILE
U
*EL FILE
S, E
*NODE PRINT , NSET=Nall 
U 
*EL PRINT , ELSET=Eall 
S 
*END STEP 

and for 15x10x10 box (1N/mm^2(?) force on the 9 node 15x10 face):

Code: Select all

*Node , NSET=Nall
1,15,10,0
2,15,10,10
3,15,0,0
4,15,0,10
5,0,10,0
6,0,10,10
7,0,0,0
8,0,0,10
9,0,5,10
10,7.5,10,10
11,15,5,10
12,7.5,0,10
13,0,5,0
14,7.5,10,0
15,15,5,0
16,7.5,0,0
17,15,10,5
18,0,10,5
19,15,0,5
20,0,0,5
21,7.5,5,10
22,7.5,5,0
23,7.5,10,5
24,7.5,0,5
25,15,5,5
26,0,5,5
*Element, TYPE=C3D10, ELSET=Eall
25,7,4,6,8,24,21,26,20,12,9,
26,7,6,4,1,26,21,24,22,23,25,
27,4,2,1,6,11,17,25,21,10,23,
28,1,4,7,3,25,24,22,15,19,16,
29,5,7,1,6,13,22,14,18,26,23,


*NSET,NSET=FemConstraintFixed
3,
4,
7,
8,
12,
16,
19,
20,
24,


*NSET,NSET=FemConstraintForce
1,
2,
5,
6,
10,
14,
17,
18,
23,


*MATERIAL, Name=Steel
*ELASTIC 
200000000.000, 0.300
*SOLID SECTION, Elset=Eall, Material=Steel
*INITIAL CONDITIONS, TYPE=STRESS, USER
*STEP
*STATIC
*BOUNDARY
FemConstraintFixed,1,3,0.0
*CLOAD
FemConstraintForce,1,0.0
FemConstraintForce,2,111.11111111111111
FemConstraintForce,3,0.0
*NODE FILE
U
*EL FILE
S, E
*NODE PRINT , NSET=Nall 
U 
*EL PRINT , ELSET=Eall 
S 
*END STEP 

FemConstraintForce is equal in the both cases, despite different face area. On the other hand, stress for 15x10 is not lower than for 10x10 face.

Hi,
I had no time so far to use the new unit parser on the FEM workbench or the material module. So at the moment all units have to be mm/kg/deg. That involves a good deal of calculation (you can use the unit calculator under the tools menu). In the long run you can use all supported units in the materials or the postprocessing.

What is the force unit? kgf (kG)?

[ickby]

I may be wrong here as my experiance with structural fem and mechanical calculations in general is limited, but I would say your analysis is faulty due to a way too coarse mesh. What I would expect is a nearly uniform stress level over the whole cube, and it must be the almost the same stress for both cubes (the average value equal the given [N/m^2] ). This is not the result of the fem calculation you have shown. What you get is an stress distribution over the elements as is to be expected with only 4 of them. Could you increase the number of elements to something around 100 and redo this analysis? I would be interested in the result!

[kwahoo]

Few thousands tet elements. 1N/mm^2.

10x10x10
Stress

Displacement


10x15x10
Stress

Displacement


Now, I'm really surprised. Why the average von Mises-Huber stresses are so different? And why they are bigger in the "10x15" case?

[ickby]

hm strange. maybe I have a missunderstanding in vanMisses stress as opossed to standart stress. displacement seems correct...
Maybe we need input from someone with more knowledge in this matter...

[ulrich1a]

I made two or three test calculations with netgen and calculix.
I just used the geometry from FreeCAD, as at that time the FEM-workbench had some limitations.
I used this webpage as starter: http://www.libremechanics.com/?q=node/19
The problem was similar, to what I wanted calculate.

I used *DLOAD and not *CLOAD
*DLOAD is distributed load and applied to faces. With *DLOAD a pressure can be simulated.
*CLOAD means concentrated load and is applied as a force to a node. I have no experience with it. But I think, depending on the node distribution, the forces at each node may have to be calculated individually for an equal distributed stress.

For an eqally distributed stress over the cross-section, I would use *DLOAD.

Stretching the material also involves some contraction. Depending on the degree of freedom for the fixed nodes, also forces in other directions may occur. The nature of tets as elements do also create forces in other directions.
I did only see equally distributed stress, with a look at stress in z-axis. Overall stress had also some components from contraction of the elements.

For the problem, I wanted to solve, I had to use C3D20 elements and had to do the meshing with cgx, which was more time consuming. I could not get a solution from tetrahedrals due to low memory on my system.

Ulrich