Static FEA of a Flanged Pressure Vessel

About the development of the FEM module/workbench.

Moderator: bernd

Post Reply
Doug
Posts: 46
Joined: Sun Feb 26, 2017 8:51 pm

Static FEA of a Flanged Pressure Vessel

Post by Doug »

Hello fellow FreeCAD'ers,

Using FreeCAD 0.16, I wrote a Python script that built a simple model, meshed it, solved it with CalculiX and displayed the results. However, as I made the model more complex, I began having trouble with huge meshes. Following the guidance from this forum, I switched to v0.17 and GMSH. I am able to get as far as solving with CalculiX, but the following error appears when trying to display the results:
FEM: Result node numbers are not equal to FEM Mesh NodeCount.

Upon inspection, I can see the result object is populated. Here is the script I am using:

Code: Select all

import sys
import ObjectsFem
from Part import makeSphere as Sphere
	
d = App.newDocument("SphereFEA")
g = FreeCADGui.ActiveDocument

d.addObject('Part::Feature', 'Chamber').Shape = Sphere(100).cut(Sphere(90))
m = ObjectsFem.makeAnalysis('MechanicalAnalysis')

ChamberMesh = ObjectsFem.makeMeshGmsh('ChamberMesh')
d.ChamberMesh.Part = d.Chamber
d.ChamberMesh.CharacteristicLengthMax = 10
d.ChamberMesh.CharacteristicLengthMin =  1
d.ChamberMesh.FemMesh
d.ChamberMesh.ElementDimension = '3D'
m.Member += [ChamberMesh]

d.addObject("Fem::ConstraintFixed","FemConstraintFixed")
d.FemConstraintFixed.Scale = 1
d.FemConstraintFixed.References = [(d.Chamber,"Vertex1")]
m.Member += [d.FemConstraintFixed]

d.addObject("Fem::ConstraintPressure","FemConstraintPressure")
d.FemConstraintPressure.Pressure = 10.0
d.FemConstraintPressure.Reversed = False
d.FemConstraintPressure.References = [(d.Chamber,"Face2"),]
m.Member += [d.FemConstraintPressure]

Metal1 = ObjectsFem.makeMaterialSolid('SolidMaterial')
mat = Metal1.Material
mat['Name'] = "Steel-Generic"
mat['YoungsModulus'] = "210000 MPa"
mat['PoissonRatio'] = "0.30"
mat['Density'] = "7900 kg/m^3"
Metal1.Material = mat
m.Member += [Metal1]

Solver = ObjectsFem.makeSolverCalculix('CalculiX')
m.Member += [Solver]

d.recompute()
g.ActiveView.fitAll()
I have to perform 2 actions in the GUI after running the script because I cannot find the appropriate Python commands:
  • * Generate the mesh
    * Run Calculix
I get the same error message in 0.17.10490. There is a reference to this error here: https://forum.freecadweb.org/viewtopic. ... 7&start=30, but in the context of using Fenics as a solver. Can someone offer a solution?

Cheers,
Doug

PS Thank you to all the developers for making this powerful tool!

_________________________________________________________________________________________________________________________
OS: Windows 7
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.17.10640 (Git)
Build type: Release
Branch: master
Hash: 2365dd6f8b424ecf01c16b70b3217ccc53098070
Python version: 2.7.8
Qt version: 4.8.7
Coin version: 4.0.0a
OCC version: 7.0.0
Last edited by Doug on Sat Mar 31, 2018 6:44 pm, edited 1 time in total.
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Error Displaying CalculiX Results in v0.17

Post by bernd »

Doug wrote: I have to perform 2 actions in the GUI after running the script because I cannot find the appropriate Python commands:
  • * Generate the mesh
    * Run Calculix
FEM_Tutorial_Python is your friend :)
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Error Displaying CalculiX Results in v0.17

Post by bernd »

Doug wrote: FEM: Result node numbers are not equal to FEM Mesh NodeCount.
Attached the example file and some code to reproduce ...

Code: Select all

len(App.ActiveDocument.ChamberMesh.FemMesh.Nodes)
len(App.ActiveDocument.CalculiX_static_results.NodeNumbers)
output:

Code: Select all

>>> 
>>> len(App.ActiveDocument.ChamberMesh.FemMesh.Nodes)
4358
>>> len(App.ActiveDocument.CalculiX_static_results.NodeNumbers)
4354
>>> 
node-number-error.fcstd
(590.12 KiB) Downloaded 77 times
For some reason CalculiX seam to return less nodes in the frd file than there are in the inp file. Question to the CalculiX cracks ! Does some of you know whats the reason for this in ccx is?
ulrich1a
Veteran
Posts: 1957
Joined: Sun Jul 07, 2013 12:08 pm

Re: Error Displaying CalculiX Results in v0.17

Post by ulrich1a »

The calculix manual says at *node file: https://www.feacluster.com/CalculiX/ccx ... de265.html

Code: Select all

Notice that only values in nodes belonging to elements are stored. Values in nodes not belonging to any element (e.g. the rotational node in a *RIGID BODY option) can only be obtained using *NODE PRINT. 
The node 5 in the example file from Bernd seems not be used in any element definition. The coordinates are nearly identical to node 1. So I would expect, that node 5 is not contained in the results and causes this kind of problem.


Ulrich
fandaL
Posts: 440
Joined: Thu Jul 24, 2014 8:29 am

Re: Error Displaying CalculiX Results in v0.17

Post by fandaL »

There seems to be some kind of a bug like here https://forum.freecadweb.org/viewtopic.php?f=18&t=21349
When I export this mesh to .unv and open in Gmsh or Salome, volume elements have mismash node numbering. Same if I create a new mesh in fc fem 3D example (attached). But it is correct for 1st order mesh.
box_2nd_order.unv.txt
(46.91 KiB) Downloaded 77 times
OS: Windows 7
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.17.10659 (Git)
Build type: Release
Branch: master
Hash: 0ca5ebe78416aa48aeb56f21699c4c7d762ccd51
Python version: 2.7.8
Qt version: 4.8.7
Coin version: 4.0.0a
OCC version: 7.0.0
Doug
Posts: 46
Joined: Sun Feb 26, 2017 8:51 pm

Re: Error Displaying CalculiX Results in v0.17

Post by Doug »

Hi Guys,

Thanks very much for the speedy response! :)
bernd wrote:
Doug wrote: I have to perform 2 actions in the GUI after running the script because I cannot find the appropriate Python commands:
  • * Generate the mesh
    * Run Calculix
FEM_Tutorial_Python is your friend :)
Bernd, indeed, the FEM_Tutorial_Python is most helpful! I used it to make my first FEM script that worked fine with 0.16 (some of the code I posted is copied directly from the tutorial). Unfortunately, when I run FEM_Tutorial_Python in 0.17, I get the following errors:

<type 'exceptions.ImportError'>: No module named FemAnalysis
<type 'exceptions.ImportError'>: No module named FemSolverCalculix
<type 'exceptions.ImportError'>: No module named FemMaterial


These modules do not appear to be included in this build. I stepped through the actions required in the GUI and using the python console output I was able to find most of the commands I needed. I believe these commands are the ones you refer to in this post: https://forum.freecadweb.org/viewtopic.php?f=18&t=20826. However, two actions (building the mesh and running CalculiX) did not send any output to the console, so I'm stumped on these commands.

Once I get this script working, I was going to offer to update the tutorial to work with your new commands, if that is desirable.

Cheers,
Doug
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Error Displaying CalculiX Results in v0.17

Post by bernd »

@Doug: thanks for the hint, just fixed FEM_Tutorial_Python and yes https://forum.freecadweb.org/viewtopic.php?f=18&t=20826 was the culprit ...
Last edited by bernd on Fri Mar 31, 2017 8:00 pm, edited 1 time in total.
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Error Displaying CalculiX Results in v0.17

Post by bernd »

make all in python ... 0.17.10664

Code: Select all

import Part, FreeCADGui, ObjectsFem, FemGui, FemGmshTools, FemToolsCcx

doc = App.newDocument("Scripted_CalculiX_Cantilever3D")
box_obj = doc.addObject('Part::Box', 'Box')
box_obj.Height = box_obj.Width = 1000
box_obj.Length = 8000
FreeCADGui.ActiveDocument.activeView().viewAxonometric()
FreeCADGui.SendMsgToActiveView("ViewFit")
FreeCADGui.ActiveDocument.activeView().viewAxonometric()
FreeCADGui.SendMsgToActiveView("ViewFit")
analysis_object = ObjectsFem.makeAnalysis("Analysis")
solver_object = ObjectsFem.makeSolverCalculix("CalculiX")
solver_object.GeometricalNonlinearity = 'linear'
solver_object.ThermoMechSteadyState = True
solver_object.MatrixSolverType = 'default'
solver_object.IterationsControlParameterTimeUse = False
doc.Analysis.Member = doc.Analysis.Member + [solver_object]
material_object = ObjectsFem.makeMaterialSolid("SolidMaterial")
mat = material_object.Material
mat['Name'] = "Steel-Generic"
mat['YoungsModulus'] = "210000 MPa"
mat['PoissonRatio'] = "0.30"
mat['Density'] = "7900 kg/m^3"
material_object.Material = mat
doc.Analysis.Member = doc.Analysis.Member + [material_object]
fixed_constraint = ObjectsFem.makeConstraintFixed("FemConstraintFixed")
fixed_constraint.References = [(doc.Box, "Face1")]
doc.Analysis.Member = doc.Analysis.Member + [fixed_constraint]
force_constraint = ObjectsFem.makeConstraintForce("FemConstraintForce")
force_constraint.References = [(doc.Box, "Face2")]
force_constraint.Force = 9000000.0
force_constraint.Direction = (doc.Box, ["Edge5"])
force_constraint.Reversed = True
doc.Analysis.Member = doc.Analysis.Member + [force_constraint]
femmesh_obj_gmsh = ObjectsFem.makeMeshGmsh('MyGMSHMeshObj')
doc.Analysis.Member = doc.Analysis.Member + [femmesh_obj_gmsh]
femmesh_obj_gmsh.Part = box_obj
App.ActiveDocument.recompute()
gmsh_mesh = FemGmshTools.FemGmshTools(femmesh_obj_gmsh)
error = gmsh_mesh.create_mesh()
print error
App.ActiveDocument.recompute()
FemGui.setActiveAnalysis(FreeCAD.ActiveDocument.Analysis)
fea = FemToolsCcx.FemToolsCcx()
fea.update_objects()
message = fea.check_prerequisites()
if not message:
    fea.reset_all()
    fea.run()
    fea.load_results()
else:
    FreeCAD.Console.PrintError("Houston, we have a problem! {}\n".format(message))  # in report view
    print("Houston, we have a problem! {}\n".format(message))  # in python console

User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Error Displaying CalculiX Results in v0.17

Post by bernd »

fandaL wrote:There seems to be some kind of a bug like here https://forum.freecadweb.org/viewtopic.php?f=18&t=21349
When I export this mesh to .unv and open in Gmsh or Salome, volume elements have mismash node numbering. Same if I create a new mesh in fc fem 3D example (attached). But it is correct for 1st order mesh.
box_2nd_order.unv.txt
OS: Windows 7
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.17.10659 (Git)
Build type: Release
Branch: master
Hash: 0ca5ebe78416aa48aeb56f21699c4c7d762ccd51
Python version: 2.7.8
Qt version: 4.8.7
Coin version: 4.0.0a
OCC version: 7.0.0
see https://forum.freecadweb.org/viewtopic.php?f=18&t=21349
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Error Displaying CalculiX Results in v0.17

Post by bernd »

ulrich1a wrote:The calculix manual says at *node file: https://www.feacluster.com/CalculiX/ccx ... de265.html

Code: Select all

Notice that only values in nodes belonging to elements are stored. Values in nodes not belonging to any element (e.g. the rotational node in a *RIGID BODY option) can only be obtained using *NODE PRINT. 
The node 5 in the example file from Bernd seems not be used in any element definition. The coordinates are nearly identical to node 1. So I would expect, that node 5 is not contained in the results and causes this kind of problem.
load the file, and run the code ...

Code: Select all

for n in App.ActiveDocument.ChamberMesh.FemMesh.Nodes:
    node_is_in_element = False
    for e in App.ActiveDocument.ChamberMesh.FemMesh.Volumes:
        if n in App.ActiveDocument.ChamberMesh.FemMesh.getElementNodes(e):
            node_is_in_element = True
            break
    if not node_is_in_element:
        print n
output:

Code: Select all

>>>
>>> for n in App.ActiveDocument.ChamberMesh.FemMesh.Nodes:
...     node_is_in_element = False
...     for e in App.ActiveDocument.ChamberMesh.FemMesh.Volumes:
...         if n in App.ActiveDocument.ChamberMesh.FemMesh.getElementNodes(e):
...             node_is_in_element = True
...             break
...     if not node_is_in_element:
...         print n
... 
5
37
38
68
>>> 
Thanks ulrich. The question is, where do we catch this?
Post Reply