Use FEA for beams

Use the pychrono.fea module to simulate flexible beams. Uses PyChrono.

Learn how to:

  • use the python.fea module
  • create beams with constraints
1 #------------------------------------------------------------------------------
2 # Name: pychrono example
3 # Purpose:
4 #
5 # Author: Alessandro Tasora
6 #
7 # Created: 1/01/2019
8 # Copyright: (c) ProjectChrono 2019
9 #------------------------------------------------------------------------------
10 
11 
12 import pychrono as chrono
13 import pychrono.fea as fea
14 import pychrono.mkl as mkl
15 import pychrono.irrlicht as chronoirr
16 
17 print ("Example: PyChrono using beam finite elements");
18 
19 # Change this path to asset path, if running from other working dir.
20 # It must point to the data folder, containing GUI assets (textures, fonts, meshes, etc.)
21 chrono.SetChronoDataPath("../../../data/")
22 
23 
24 # ---------------------------------------------------------------------
25 #
26 # Create the simulation system and add items
27 #
28 
29 
30 # Create a Chrono::Engine physical system
31 my_system = chrono.ChSystemNSC()
32 
33 ## Create a mesh, that is a container for groups
34 ## of elements and their referenced nodes.
35 my_mesh = fea.ChMesh();
36 
37 ## Create a section, i.e. thickness and material properties
38 ## for beams. This will be shared among some beams.
39 
40 msection = fea.ChBeamSectionAdvanced()
41 
42 beam_wy = 0.012
43 beam_wz = 0.025
44 msection.SetAsRectangularSection(beam_wy, beam_wz)
45 msection.SetYoungModulus(0.01e9)
46 msection.SetGshearModulus(0.01e9 * 0.3)
47 msection.SetBeamRaleyghDamping(0.000)
48 #msection.SetCentroid(0,0.02)
49 #msection.SetShearCenter(0,0.1)
50 #msection.SetSectionRotation(45*chrono.CH_C_RAD_TO_DEG)
51 
52 # Add some EULER-BERNOULLI BEAMS:
53 
54 beam_L = 0.1
55 
56 hnode1 = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(0, 0, 0)))
57 hnode2 = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(beam_L, 0, 0)))
58 hnode3 = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(beam_L * 2, 0, 0)))
59 
60 my_mesh.AddNode(hnode1)
61 my_mesh.AddNode(hnode2)
62 my_mesh.AddNode(hnode3)
63 
64 belement1 = fea.ChElementBeamEuler()
65 
66 belement1.SetNodes(hnode1, hnode2)
67 belement1.SetSection(msection)
68 
69 my_mesh.AddElement(belement1)
70 
71 belement2 = fea.ChElementBeamEuler()
72 
73 belement2.SetNodes(hnode2, hnode3)
74 belement2.SetSection(msection)
75 
76 my_mesh.AddElement(belement2);
77 
78 # Apply a force or a torque to a node:
79 hnode2.SetForce(chrono.ChVectorD(4, 2, 0))
80 hnode3.SetTorque(chrono.ChVectorD(0, -0.04, 0))
81 
82 # Fix a node to ground:
83 # hnode1.SetFixed(True)
84 # otherwise fix it using constraints:
85 
86 mtruss = chrono.ChBody()
87 mtruss.SetBodyFixed(True)
88 my_system.Add(mtruss)
89 
90 constr_bc = chrono.ChLinkMateGeneric()
91 constr_bc.Initialize(hnode3, mtruss, False, hnode3.Frame(), hnode3.Frame())
92 my_system.Add(constr_bc)
93 constr_bc.SetConstrainedCoords(True, True, True, # x, y, z
94  True, True, True) # Rx, Ry, Rz
95 
96 constr_d = chrono.ChLinkMateGeneric()
97 constr_d.Initialize(hnode1, mtruss, False, hnode1.Frame(), hnode1.Frame())
98 my_system.Add(constr_d)
99 constr_d.SetConstrainedCoords(False, True, True, # x, y, z
100  False, False,False) # Rx, Ry, Rz
101 
102 
103 # Add some EULER-BERNOULLI BEAMS (the fast way!)
104 
105 # Shortcut!
106 # This ChBuilderBeam helper object is very useful because it will
107 # subdivide 'beams' into sequences of finite elements of beam type, ex.
108 # one 'beam' could be made of 5 FEM elements of ChElementBeamEuler class.
109 # If new nodes are needed, it will create them for you.
110 
111 builder = fea.ChBuilderBeam()
112 
113 # Now, simply use BuildBeam to create a beam from a point to another:
114 builder.BuildBeam(my_mesh, # the mesh where to put the created nodes and elements
115  msection, # the ChBeamSectionAdvanced to use for the ChElementBeamEuler elements
116  5, # the number of ChElementBeamEuler to create
117  chrono.ChVectorD(0, 0, -0.1), # the 'A' point in space (beginning of beam)
118  chrono.ChVectorD(0.2, 0, -0.1), # the 'B' point in space (end of beam)
119  chrono.ChVectorD(0, 1, 0)) # the 'Y' up direction of the section for the beam
120 
121 ## After having used BuildBeam(), you can retrieve the nodes used for the beam,
122 ## For example say you want to fix the A end and apply a force to the B end:
123 builder.GetLastBeamNodes().back().SetFixed(True)
124 builder.GetLastBeamNodes().front().SetForce(chrono.ChVectorD(0, -1, 0))
125 
126 # Again, use BuildBeam for creating another beam, this time
127 # it uses one node (the last node created by the last beam) and one point:
128 builder.BuildBeam(my_mesh, msection, 5,
129  builder.GetLastBeamNodes().front(), # the 'A' node in space (beginning of beam)
130  chrono.ChVectorD(0.2, 0.1, -0.1), # the 'B' point in space (end of beam)
131  chrono.ChVectorD(0, 1, 0)); # the 'Y' up direction of the section for the beam
132 
133 
134 # We do not want gravity effect on FEA elements in this demo
135 my_mesh.SetAutomaticGravity(False);
136 
137 # Remember to add the mesh to the system!
138 my_system.Add(my_mesh)
139 
140 # ==Asset== attach a visualization of the FEM mesh.
141 # This will automatically update a triangle mesh (a ChTriangleMeshShape
142 # asset that is internally managed) by setting proper
143 # coordinates and vertex colors as in the FEM elements.
144 
145 mvisualizebeamA = fea.ChVisualizationFEAmesh(my_mesh)
146 mvisualizebeamA.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_ELEM_BEAM_MZ)
147 mvisualizebeamA.SetColorscaleMinMax(-0.4, 0.4)
148 mvisualizebeamA.SetSmoothFaces(True)
149 mvisualizebeamA.SetWireframe(False)
150 my_mesh.AddAsset(mvisualizebeamA)
151 
152 mvisualizebeamC = fea.ChVisualizationFEAmesh(my_mesh)
153 mvisualizebeamC.SetFEMglyphType(fea.ChVisualizationFEAmesh.E_GLYPH_NODE_CSYS)
154 mvisualizebeamC.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_NONE)
155 mvisualizebeamC.SetSymbolsThickness(0.006)
156 mvisualizebeamC.SetSymbolsScale(0.01)
157 mvisualizebeamC.SetZbufferHide(False)
158 my_mesh.AddAsset(mvisualizebeamC)
159 
160 
161 # ---------------------------------------------------------------------
162 #
163 # Create an Irrlicht application to visualize the system
164 #
165 
166 
167 # Create the Irrlicht visualization (open the Irrlicht device,
168 # bind a simple user interface, etc. etc.)
169 myapplication = chronoirr.ChIrrApp(my_system, 'Test FEA beams', chronoirr.dimension2du(1024,768))
170 
171 #application.AddTypicalLogo()
172 myapplication.AddTypicalSky()
173 myapplication.AddTypicalLogo(chrono.GetChronoDataPath() + 'logo_pychrono_alpha.png')
174 myapplication.AddTypicalCamera(chronoirr.vector3df(0.1,0.1,0.2))
175 myapplication.AddTypicalLights()
176 
177 # ==IMPORTANT!== Use this function for adding a ChIrrNodeAsset to all items
178 # in the system. These ChIrrNodeAsset assets are 'proxies' to the Irrlicht meshes.
179 
180 myapplication.AssetBindAll()
181 
182 # ==IMPORTANT!== Use this function for 'converting' into Irrlicht meshes the assets
183 # that you added to the bodies into 3D shapes, they can be visualized by Irrlicht!
184 
185 myapplication.AssetUpdateAll()
186 
187 # Mark completion of system construction
188 my_system.SetupInitial();
189 
190 
191 # THE SOFT-REAL-TIME CYCLE
192 
193 
194 # Change the solver form the default SOR to the MKL Pardiso, more precise for fea.
195 msolver = mkl.ChSolverMKLcsm()
196 my_system.SetSolver(msolver)
197 
198 
199 myapplication.SetTimestep(0.001);
200 
201 # application.GetSystem().DoStaticLinear()
202 
203 while(myapplication.GetDevice().run()):
204  myapplication.BeginScene()
205  myapplication.DrawAll()
206  myapplication.DoStep()
207  myapplication.EndScene()
208 
209 
210 
211 
212 
213