Create a slider-crank system

Create a slider-crank, visulize and simulate it, and finally plot results using Python plotting tools. Uses PyChrono.

Learn how to:

  • create a simple slider-crank mechanism
  • add a motor with imposed angular speed
  • plot results using python's matplotlib library
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.core as chrono
13 import pychrono.irrlicht as chronoirr
14 import matplotlib.pyplot as plt
15 import numpy as np
16 
17 print ("Example: create a slider crank and plot results");
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 # Create the simulation system and add items
26 #
27 
28 mysystem = chrono.ChSystemNSC()
29 
30 # Some data shared in the following
31 crank_center = chrono.ChVectorD(-1,0.5,0)
32 crank_rad = 0.4
33 crank_thick = 0.1
34 rod_length = 1.5
35 
36 
37 # Create four rigid bodies: the truss, the crank, the rod, the piston.
38 
39 # Create the floor truss
40 mfloor = chrono.ChBodyEasyBox(3, 1, 3, 1000)
41 mfloor.SetPos(chrono.ChVectorD(0,-0.5,0))
42 mfloor.SetBodyFixed(True)
43 mysystem.Add(mfloor)
44 
45 # Create the flywheel crank
46 mcrank = chrono.ChBodyEasyCylinder(crank_rad, crank_thick, 1000)
47 mcrank.SetPos(crank_center + chrono.ChVectorD(0, 0, -0.1))
48 # Since ChBodyEasyCylinder creates a vertical (y up) cylinder, here rotate it:
49 mcrank.SetRot(chrono.Q_ROTATE_Y_TO_Z)
50 mysystem.Add(mcrank)
51 
52 # Create a stylized rod
53 mrod = chrono.ChBodyEasyBox(rod_length, 0.1, 0.1, 1000)
54 mrod.SetPos(crank_center + chrono.ChVectorD(crank_rad+rod_length/2 , 0, 0))
55 mysystem.Add(mrod)
56 
57 # Create a stylized piston
58 mpiston = chrono.ChBodyEasyCylinder(0.2, 0.3, 1000)
59 mpiston.SetPos(crank_center + chrono.ChVectorD(crank_rad+rod_length, 0, 0))
60 mpiston.SetRot(chrono.Q_ROTATE_Y_TO_X)
61 mysystem.Add(mpiston)
62 
63 
64 # Now create constraints and motors between the bodies.
65 
66 # Create crank-truss joint: a motor that spins the crank flywheel
68 my_motor.Initialize(mcrank, # the first connected body
69  mfloor, # the second connected body
70  chrono.ChFrameD(crank_center)) # where to create the motor in abs.space
71 my_angularspeed = chrono.ChFunction_Const(chrono.CH_C_PI) # ang.speed: 180°/s
72 my_motor.SetMotorFunction(my_angularspeed)
73 mysystem.Add(my_motor)
74 
75 # Create crank-rod joint
76 mjointA = chrono.ChLinkLockRevolute()
77 mjointA.Initialize(mrod,
78  mcrank,
79  chrono.ChCoordsysD( crank_center + chrono.ChVectorD(crank_rad,0,0) ))
80 mysystem.Add(mjointA)
81 
82 # Create rod-piston joint
83 mjointB = chrono.ChLinkLockRevolute()
84 mjointB.Initialize(mpiston,
85  mrod,
86  chrono.ChCoordsysD( crank_center + chrono.ChVectorD(crank_rad+rod_length,0,0) ))
87 mysystem.Add(mjointB)
88 
89 # Create piston-truss joint
91 mjointC.Initialize(mpiston,
92  mfloor,
93  chrono.ChCoordsysD(
94  crank_center + chrono.ChVectorD(crank_rad+rod_length,0,0),
95  chrono.Q_ROTATE_Z_TO_X)
96  )
97 mysystem.Add(mjointC)
98 
99 
100 
101 # ---------------------------------------------------------------------
102 #
103 # Create an Irrlicht application to visualize the system
104 #
105 
106 myapplication = chronoirr.ChIrrApp(mysystem, 'PyChrono example', chronoirr.dimension2du(1024,768))
107 
108 myapplication.AddTypicalSky()
109 myapplication.AddTypicalLogo(chrono.GetChronoDataPath() + 'logo_pychrono_alpha.png')
110 myapplication.AddTypicalCamera(chronoirr.vector3df(1,1,3), chronoirr.vector3df(0,1,0))
111 myapplication.AddTypicalLights()
112 
113  # ==IMPORTANT!== Use this function for adding a ChIrrNodeAsset to all items
114  # in the system. These ChIrrNodeAsset assets are 'proxies' to the Irrlicht meshes.
115  # If you need a finer control on which item really needs a visualization proxy in
116  # Irrlicht, just use application.AssetBind(myitem); on a per-item basis.
117 
118 myapplication.AssetBindAll();
119 
120  # ==IMPORTANT!== Use this function for 'converting' into Irrlicht meshes the assets
121  # that you added to the bodies into 3D shapes, they can be visualized by Irrlicht!
122 
123 myapplication.AssetUpdateAll();
124 
125 
126 # ---------------------------------------------------------------------
127 #
128 # Run the simulation
129 #
130 
131 # Initialize these lists to store values to plot.
132 array_time = []
133 array_angle = []
134 array_pos = []
135 array_speed = []
136 
137 myapplication.SetTimestep(0.005)
138 
139 # Run the interactive simulation loop
140 while(myapplication.GetDevice().run()):
141 
142  # for plotting, append instantaneous values:
143  array_time.append(mysystem.GetChTime())
144  array_angle.append(my_motor.GetMotorRot())
145  array_pos.append(mpiston.GetPos().x)
146  array_speed.append(mpiston.GetPos_dt().x)
147 
148  # here happens the visualization and step time integration
149  myapplication.BeginScene()
150  myapplication.DrawAll()
151  myapplication.DoStep()
152  myapplication.EndScene()
153 
154  # stop simulation after 2 seconds
155  if mysystem.GetChTime() > 2:
156  myapplication.GetDevice().closeDevice()
157 
158 
159 # Use matplotlib to make two plots when simulation ended:
160 fig, (ax1, ax2) = plt.subplots(2, sharex = True)
161 
162 ax1.plot(array_angle, array_pos)
163 ax1.set(ylabel='position [m]')
164 ax1.grid()
165 
166 ax2.plot(array_angle, array_speed, 'r--')
167 ax2.set(ylabel='speed [m]',xlabel='angle [rad]')
168 ax2.grid()
169 
170 # trick to plot \pi on x axis of plots instead of 1 2 3 4 etc.
171 plt.xticks(np.linspace(0, 2*np.pi, 5),['0','$\pi/2$','$\pi$','$3\pi/2$','$2\pi$'])
172 
173 
174 
175 
176