Use FEA for a rotor

Use the pychrono.fea module to simulate the Jeffcott rotor Uses PyChrono.

Learn how to:

  • use the python.fea module
  • create a flexible rotor passing through instability
  • tweak the integrator and solver settings for higher precision
  • create an ad-hoc motion function by python-side inheritance from ChFunction
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 math as m
13 import pychrono as chrono
14 import pychrono.fea as fea
15 import pychrono.mkl as mkl
16 import pychrono.irrlicht as chronoirr
17 
18 print ("Example: FEA of the Jeffcott rotor passing through resonance.");
19 
20 # Change this path to asset path, if running from other working dir.
21 # It must point to the data folder, containing GUI assets (textures, fonts, meshes, etc.)
22 chrono.SetChronoDataPath("../../../data/")
23 
24 
25 # ---------------------------------------------------------------------
26 #
27 # Create the simulation system and add items
28 #
29 
30 
31 # Create a Chrono::Engine physical system
32 my_system = chrono.ChSystemNSC()
33 
34 my_mesh = fea.ChMesh()
35 my_system.Add(my_mesh)
36 
37 my_mesh.SetAutomaticGravity(True,2) # for max precision in gravity of FE, at least 2 integration points per element when using cubic IGA
38 my_system.Set_G_acc(chrono.ChVectorD(0,-9.81, 0));
39 
40 beam_L = 6
41 beam_ro = 0.050
42 beam_ri = 0.045
43 CH_C_PI = 3.1456
44 
45 # Create a section, i.e. thickness and material properties
46 # for beams. This will be shared among some beams.
47 
48 melasticity = fea.ChElasticityCosseratSimple()
49 melasticity.SetYoungModulus(210e9)
50 melasticity.SetGwithPoissonRatio(0.3)
51 melasticity.SetIyy( (CH_C_PI / 4.0) * (pow(beam_ro, 4) - pow(beam_ri, 4)) )
52 melasticity.SetIzz( (CH_C_PI / 4.0) * (pow(beam_ro, 4) - pow(beam_ri, 4)) )
53 melasticity.SetJ ( (CH_C_PI / 2.0) * (pow(beam_ro, 4) - pow(beam_ri, 4)) )
54 
55 msection = fea.ChBeamSectionCosserat(melasticity)
56 msection.SetDensity(7800)
57 msection.SetCircular(True)
58 msection.SetDrawCircularRadius(beam_ro) # SetAsCircularSection(..) would overwrite Ixx Iyy J etc.
59 msection.SetArea(CH_C_PI * (pow(beam_ro, 2)- pow(beam_ri, 2)))
60 
61 # Use the ChBuilderBeamIGA tool for creating a straight rod
62 # divided in Nel elements:
63 
64 builder = fea.ChBuilderBeamIGA()
65 builder.BuildBeam(my_mesh, # the mesh to put the elements in
66  msection, # section of the beam
67  20, # number of sections (spans)
68  chrono.ChVectorD(0, 0, 0), # start point
69  chrono.ChVectorD(beam_L, 0, 0), # end point
70  chrono.VECT_Y, # suggested Y direction of section
71  1) # order (3 = cubic, etc)
72 
73 node_mid = builder.GetLastBeamNodes()[m.floor(builder.GetLastBeamNodes().size()/2.0)]
74 
75 # Create the flywheel and attach it to the center of the beam
76 
77 mbodyflywheel = chrono.ChBodyEasyCylinder(0.24, 0.05, 7800) # R, h, density
78 mbodyflywheel.SetCoord(
79  chrono.ChCoordsysD(node_mid.GetPos() + chrono.ChVectorD(0,0.05,0), # flywheel initial center (plus Y offset)
80  chrono.Q_from_AngAxis(CH_C_PI/2.0, chrono.VECT_Z)) # flywheel initial alignment (rotate 90° so cylinder axis is on X)
81 )
82 my_system.Add(mbodyflywheel)
83 
84 myjoint = chrono.ChLinkMateFix()
85 myjoint.Initialize(node_mid, mbodyflywheel)
86 my_system.Add(myjoint)
87 
88 # Create the truss
89 truss = chrono.ChBody()
90 truss.SetBodyFixed(True)
91 my_system.Add(truss)
92 
93 # Create the end bearing
94 bearing = chrono.ChLinkMateGeneric(False, True, True, False, True, True)
95 bearing.Initialize(builder.GetLastBeamNodes().back(),
96  truss,
97  chrono.ChFrameD(builder.GetLastBeamNodes().back().GetPos())
98 )
99 my_system.Add(bearing)
100 
101 # Create the motor that rotates the beam
103 
104 # Connect the rotor and the stator and add the motor to the system:
105 rotmotor1.Initialize(builder.GetLastBeamNodes().front(), # body A (slave)
106  truss, # body B (master)
107  chrono.ChFrameD(builder.GetLastBeamNodes().front().GetPos(),
108  chrono.Q_from_AngAxis(CH_C_PI/2.0, chrono.VECT_Y)) # motor frame, in abs. coords
109 )
110 my_system.Add(rotmotor1)
111 
112 # use a custom function for setting the speed of the motor
113 class ChFunction_myf (chrono.ChFunction):
114  def __init__(self):
115  chrono.ChFunction.__init__(self)
116  def Get_y(self,x):
117  A1 = 0.8
118  A2 = 1.2
119  T1 = 0.5
120  T2 = 1.0
121  T3 = 1.25
122  w = 60
123  if x < T1:
124  return A1 * w * (1. - m.cos(CH_C_PI*x / T1)) / 2.0
125  elif (x > T1 and x <= T2):
126  return A1 * w
127  elif (x > T2 and x <= T3):
128  return A1 * w + (A2 - A1) * w * (1.0 - m.cos(CH_C_PI*(x - T2) / (T3 - T2))) / 2.0
129  else:
130  return A2 * w
131 
132 f_ramp = ChFunction_myf()
133 #f_ramp = chrono.ChFunction_Sine(0,0.2,40)
134 rotmotor1.SetMotorFunction(f_ramp)
135 
136 # Attach a visualization of the FEM mesh.
137 
138 mvisualizebeamA = fea.ChVisualizationFEAmesh(my_mesh)
139 mvisualizebeamA.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_SURFACE)
140 mvisualizebeamA.SetSmoothFaces(True)
141 my_mesh.AddAsset(mvisualizebeamA)
142 
143 mvisualizebeamC = fea.ChVisualizationFEAmesh(my_mesh)
144 mvisualizebeamC.SetFEMglyphType(fea.ChVisualizationFEAmesh.E_GLYPH_NODE_CSYS)
145 mvisualizebeamC.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_NONE)
146 mvisualizebeamC.SetSymbolsThickness(0.006)
147 mvisualizebeamC.SetSymbolsScale(0.01)
148 mvisualizebeamC.SetZbufferHide(False)
149 my_mesh.AddAsset(mvisualizebeamC)
150 
151 
152 # ---------------------------------------------------------------------
153 #
154 # Create an Irrlicht application to visualize the system
155 #
156 
157 
158 # Create the Irrlicht visualization (open the Irrlicht device,
159 # bind a simple user interface, etc. etc.)
160 myapplication = chronoirr.ChIrrApp(my_system, 'Test FEA: the Jeffcott rotor with IGA beams', chronoirr.dimension2du(1024,768))
161 
162 myapplication.AddTypicalLogo(chrono.GetChronoDataPath() + 'logo_pychrono_alpha.png')
163 myapplication.AddTypicalSky()
164 myapplication.AddTypicalCamera(chronoirr.vector3df(0,1,4), chronoirr.vector3df(beam_L/2, 0, 0))
165 myapplication.AddTypicalLights()
166 
167 # This is needed if you want to see things in Irrlicht 3D view.
168 myapplication.AssetBindAll()
169 myapplication.AssetUpdateAll()
170 
171 
172 # ---------------------------------------------------------------------
173 #
174 # Run the simulation
175 #
176 
177 
178 my_system.SetupInitial();
179 
180 # Set to a more precise HHT timestepper if needed
181 # my_system.SetTimestepperType(chrono.ChTimestepper.Type_HHT)
182 
183 # Change the solver form the default SOR to the MKL Pardiso, more precise for fea.
184 msolver = mkl.ChSolverMKLcsm()
185 my_system.SetSolver(msolver)
186 
187 myapplication.SetTimestep(0.002)
188 
189 my_system.DoStaticLinear()
190 
191 while(myapplication.GetDevice().run()):
192  myapplication.BeginScene()
193  myapplication.DrawAll()
194  myapplication.DoStep()
195  myapplication.EndScene()
196