FEM Test Cases

[sgrogan]

I think once the ccxInpWriter migration is finished we should set up a few test cases for the FEM workbench. A test model, with all forces and constraints that we could run and the results could be compares with something that has been verified. A simple beam with a load could be a good start and can be easily verified. I'd like to see that, because now we have a "working" FEM wb that might in some cases produce a misleading output.

I think we have come to this point. THANK YOU to all that have advanced the FEM module so actively! The GUI is becoming quite polished.
A few test cases would be extremely useful at this point. Ideally there would be an analytical or accepted solution and maybe some results from other FEM solutions.

A cantilever by bernd is a must.
A pipe by ulrich1a.
An arch also by ulrich1a

[fandaL]

I tried to simulate ulrich1a's example of the ring viewtopic.php?f=18&t=10692&start=40 in the Nastran/Patran.
Model geometry is same, but I wasn't able to read material properties from the file, so I used E = 210 000 MPa, mu = 0.3.
Mesh in Nastran is similar: 2nd order tetrahedrons, 341 elements and 794 nodes. Ulrich1a's example has 301 elements and 754 nodes.
Loading vertical force 1000 N is distributed equaly to nodes on inner surface, respectively pressure 1000 MPa in the second loadcase.
Analytical approach for pressure is below. Formulas are from the literature Shigley, Mischke, Budynas: Mechanical Engineering Design.

analytical_approach-tube.png (11.42 KiB) Viewed 3208 times

Displacement and stress of the ring loaded by pressure:

ring_pressure_displacement.png (105.98 KiB) Viewed 3208 times

ring_pressure_stress.png (108.38 KiB) Viewed 3208 times

Ulrich1a wrote that stress in his DLOAD case was 8869,6 MPa. But nearer to Nastran 6500 MPa is 8869,6 / sqrt(2) = 6272 MPa according to kwahoo's post viewtopic.php?f=18&t=7172&start=40

but IMO, there is missing division by two. Something like

Code: Select all
mstress.append( sqrt( 0.5 * ( pow( i[0] - i[1] ,2) + pow( i[1] - i[2] ,2) + pow( i[2] - i[0] ,2) + 6 * (pow(i[3],2)+pow(i[4],2)+pow(i[5],2) ) ) ) )

Second loadcase is for the vertical force 1000 N. Displacement and stress

ring_force_displacement.png (107.7 KiB) Viewed 3208 times

ring_force_stress.png (116.05 KiB) Viewed 3208 times

Similarly Ulrich1a's stress for CLOAD 20,6 / sqrt(2) = 14,6 MPa is nearer to Nastran 16,1 MPa in neighbourhood of fix support.

[sgrogan]

Can someone look at CalculixLib.py
Line 163

Code: Select all

mstress.append( sqrt( pow( i[0] - i[1] ,2) + pow( i[1] - i[2] ,2) + pow( i[2] - i[0] ,2) + 6 * (pow(i[3],2)+pow(i[4],2)+pow(i[5],2)  )  ) )

I think we are missing a 1/2 inside the sqrt
http://en.wikipedia.org/wiki/Von_Mises_yield_criterion

[drei]

Yes, we are missing that division.

[PrzemoF]

To avoid confusion: CalculixLib.py --> ccxFrdReader.py since [1]

That formula should be split into something verifiable, possibly using variables like s11, s22, s23 (s for sigma) i.e.

Code: Select all

  s11 = i[0]
  s22 = i[1]
  s33 = i[2]
  s12 = i[3]
  s23 = i[4]
  s31 = i[5]
  s11s22 = pow(s11 - s22, 2) 
  s22s33 = pow(s22 - s33, 2) 
  s33s11 = pow(s33 - s11, 2) 
  s12s23s31 = 6 * (pow(s12, 2) + pow(s23, 2) * pow(s31, 2))
  mstress.append(sqrt(0.5 * (s11s22 + s22s33 + s33s11 + s12s23s31)))

[1] https://github.com/FreeCAD/FreeCAD_sf_m ... 1685a5c444

[PrzemoF]

viewtopic.php?f=17&t=10926

[kwahoo]

Can someone look at CalculixLib.py
Line 163

Code: Select all

mstress.append( sqrt( pow( i[0] - i[1] ,2) + pow( i[1] - i[2] ,2) + pow( i[2] - i[0] ,2) + 6 * (pow(i[3],2)+pow(i[4],2)+pow(i[5],2)  )  ) )

I think we are missing a 1/2 inside the sqrt
http://en.wikipedia.org/wiki/Von_Mises_yield_criterion

Confirmed earlier viewtopic.php?f=18&t=7172&start=40#p87668 and fixed now by Przemo.

[fandaL]

To avoid confusion: CalculixLib.py --> ccxFrdReader.py since [1]

That formula should be split into something verifiable, possibly using variables like s11, s22, s23 (s for sigma) i.e.

Code: Select all

  s11 = i[0]
  s22 = i[1]
  s33 = i[2]
  s12 = i[3]
  s23 = i[4]
  s31 = i[5]
  s11s22 = pow(s11 - s22, 2) 
  s22s33 = pow(s22 - s33, 2) 
  s33s11 = pow(s33 - s11, 2) 
  s12s23s31 = 6 * (pow(s12, 2) + pow(s23, 2) * pow(s31, 2))
  mstress.append(sqrt(0.5 * (s11s22 + s22s33 + s33s11 + s12s23s31)))

[1] https://github.com/FreeCAD/FreeCAD_sf_m ... 1685a5c444

I didn't noted it earlier. I think in von Mises formula should be + instead of * in

Code: Select all

  s12s23s31 = 6 * (pow(s12, 2) + pow(s23, 2) + pow(s31, 2))

https://github.com/fandaL/FreeCAD/commi ... 699ef750ff

[PrzemoF]

You're right... I hope that's the last error in Von Mises formula... Can you send pull request from your branch?

Edit: I'll do it as it requires test results update.
Edit2: pull request viewtopic.php?f=17&t=13203