sliderCrank3DwithANCFbeltDrive2.py
You can view and download this file on Github: sliderCrank3DwithANCFbeltDrive2.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Create slider-crank mechanism with belt drive modeled with ANCF cable elements
5#
6# Authors: Martin Knapp and Lukas March
7# Date: Created on Thu May 19 12:22:52 2020
8#
9# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
10#
11# Notes: PROJECT Exercise: Drive System + Crank System; VU Industrielle Mechatronik 2 - Robotics and Simulation
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15#
16#
17#AUTHORS: Martin Knapp & Lukas March
18#
19#
20#Copyright: This file is part of Exudyn. Exudyn is free software.
21#You can redistribute it and/or modify it under the terms of the Exudyn license.
22#See 'LICENSE.txt' for more details.
23#"""
24
25#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27#Tested with EXUDYN Version 0.1.342.
28
29import sys
30import os
31
32
33#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35import sys
36import os
37import numpy as np
38import matplotlib.pyplot as plt
39import matplotlib.ticker as ticker
40
41
42import exudyn as exu
43import numpy as np
44from exudyn.itemInterface import (MarkerNodeRotationCoordinate, ObjectConnectorCartesianSpringDamper,
45 LoadTorqueVector, VObjectJointPrismatic2D, ObjectJointPrismatic2D, Torque,
46 MassPoint2D, RigidBody2D, NodePoint2D, RevoluteJoint2D, CoordinateConstraint,
47 ObjectGround, ObjectContactFrictionCircleCable2D, NodeGenericData,
48 MarkerBodyCable2DShape, NodePoint2DSlope1, ObjectANCFCable2D, MarkerBodyPosition,
49 VObjectJointRevolute2D, VObjectRigidBody2D, NodePointGround, MarkerNodePosition,
50 MarkerNodeCoordinate, Force, SensorBody, NodeRigidBody2D, ObjectRigidBody2D,
51 MarkerBodyRigid, ObjectJointRevolute2D, SensorLoad)
52from exudyn.utilities import AddRigidBody, RigidBodyInertia, ObjectConnectorCoordinate, InertiaCuboid
53import exudyn.graphics as graphics #only import if it does not conflict
54#import visHelper
55#visHelper.init()
56
57#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
58
59PLOTS_PATH = "plots/"
60fontSize = 20
61
62#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
63# Change plot axis
64def set_axis(plt, equal = False):
65 ax=plt.gca() #get current axes
66 ax.grid(True , 'major', 'both')
67 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
68 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
69 if equal:
70 ax.set_aspect('equal')
71 plt.legend() #show labels as legend
72
73def plotOmegaDisk0():
74 ang_vel = np.loadtxt("angular_velocity_disk0.txt", comments='#', delimiter=',')
75
76 fig = plt.figure(figsize=[13,5])
77 plt.plot(ang_vel[:,0], ang_vel[:,3], 'r-', label='$\omega_{disk0}$')
78 set_axis(plt)
79 ax=plt.gca()
80 ax.set_xlabel('$t [s]$', fontsize=fontSize)
81 ax.set_ylabel('$\omega_{disk0} [rad/s]$', fontsize=fontSize)
82 ax.grid(True, 'major', 'both')
83 plt.legend()
84 plt.tight_layout()
85 plt.show()
86 fig.savefig(PLOTS_PATH + 'angular_velocity_disk0.pdf', format='pdf')
87
88def plotOmegaDisk1():
89 ang_vel = np.loadtxt("angular_velocity_disk1.txt", comments='#', delimiter=',')
90
91 fig = plt.figure(figsize=[13,5])
92 plt.plot(ang_vel[:,0], ang_vel[:,3], 'r-', label='$\omega_{disk1}$')
93 set_axis(plt)
94 ax=plt.gca()
95 ax.set_xlabel('$t [s]$', fontsize=fontSize)
96 ax.set_ylabel('$\omega_{disk1} [rad/s]$', fontsize=fontSize)
97 ax.grid(True, 'major', 'both')
98 plt.legend()
99 plt.tight_layout()
100 plt.show()
101 fig.savefig(PLOTS_PATH + 'angular_velocity_disk1.pdf', format='pdf')
102
103def plotTorque():
104 ang_vel = np.loadtxt("torque.txt", comments='#', delimiter=',')
105
106 fig = plt.figure(figsize=[13,5])
107 plt.plot(ang_vel[:,0], ang_vel[:,3], 'r-', label='$\\tau$')
108 set_axis(plt)
109 ax=plt.gca()
110 ax.set_xlabel('$t [s]$', fontsize=fontSize)
111 ax.set_ylabel('$\\tau [Nm]$', fontsize=fontSize)
112 ax.grid(True, 'major', 'both')
113 plt.legend()
114 plt.tight_layout()
115 plt.show()
116 fig.savefig(PLOTS_PATH + 'torque.pdf', format='pdf')
117
118def plotCrankPos():
119 ang_vel = np.loadtxt("crank_pos.txt", comments='#', delimiter=',')
120
121 fig = plt.figure(figsize=[13,5])
122 plt.plot(ang_vel[:,0], ang_vel[:,1], 'r-', label='$x_{Pos}$')
123 plt.plot(ang_vel[:,0], ang_vel[:,2], 'b-', label='$y_{Pos}$')
124 set_axis(plt)
125 ax=plt.gca()
126 ax.set_xlabel('$t [s]$', fontsize=fontSize)
127 ax.set_ylabel('$x,y [m]$', fontsize=fontSize)
128 ax.grid(True, 'major', 'both')
129 plt.legend()
130 plt.tight_layout()
131 plt.show()
132 fig.savefig(PLOTS_PATH + 'crank_pos.pdf', format='pdf')
133
134def plotBelt():
135 belt = np.loadtxt("belt.txt", comments='#', delimiter=',')
136 belt_slope = np.loadtxt("belt_slope.txt", comments='#', delimiter=',')
137
138 fig = plt.figure(figsize=[13,5])
139 plt.plot(belt[:,0], belt[:,1], 'r-', label='Belt')
140 scale = 100 # scale arrow length
141 for i in range(belt_slope.shape[0]-1):
142 plt.arrow(belt_slope[i,0], belt_slope[i,2],
143 belt_slope[i,1] / scale, belt_slope[i,3] / scale, zorder=10)
144 set_axis(plt, True)
145 ax=plt.gca()
146 ax.set_xlabel('$x [m]$', fontsize=fontSize)
147 ax.set_ylabel('$y [m]$', fontsize=fontSize)
148 ax.grid(True, 'major', 'both')
149 plt.legend()
150 plt.tight_layout()
151 plt.show()
152 fig.savefig(PLOTS_PATH + 'belt.pdf', format='pdf')
153
154def vishelperInit():
155 plt.close('all')
156 if not os.path.isdir(PLOTS_PATH):
157 os.mkdir(PLOTS_PATH)
158
159def visHelperPlot_all():
160 plotBelt()
161 plotOmegaDisk0()
162 plotOmegaDisk1()
163 plotTorque()
164 plotCrankPos()
165
166#if __name__ is "__main__":
167# init()
168# plot_all()
169
170
171vishelperInit()
172
173
174
175print("Exudyn used:", exu.GetVersionString())
176
177#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
178#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
179#SYSTEM SETTINGS
180
181# 0 - belt-drive-system
182# 1 - crank system 2D
183# 2 - crank system 3D
184# 3 - belt-drive-system + crank system 2D
185# 4 - belt-drive-system + crank system 3D
186
187sys_set = 4
188
189useGraphics = True
190export_images = useGraphics
191PLOTS_PATH = "plots/"
192
193if export_images:
194 if not os.path.isdir(PLOTS_PATH):
195 os.mkdir(PLOTS_PATH)
196
197# belt-drive-system
198test_belt_function = True # Draws the belt curve and the tangential vectors
199enable_force = True # Enable Preload Force in the belt
200enable_disk_friction = True # When True the disks will have contact to the belt else the belt is free floating
201enable_controller = True # If True the disk0 will get a torque regulated by a P-controller else a fixed torque
202n = 50 # Belt element count
203
204#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
205#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
206
207#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
208#geometry parameters
209L_0 = 0.5 #distance between the disks [m]
210L_A = 0.15 #length of crank [m]
211h_A = 0.025 #cross-section height and width of crank [m]
212L_B = 0.3 #length of connecting rod [m]
213h_B = 0.025 #cross-section height and width of connecting rod [m]
214r_0 = 0.05 #radius of disk0 [m]
215r_1 = 0.1 #radius of disk1 [m]
216b_a0 = 0.1 #section-length of crank [m]
217b_a1 = 0.05 #section-length of crank [m]
218b_a2 = 0.05 #section-length of crank [m]
219h_S = 0.1 #cross-section height and width of slider [m]
220#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
221#belt parameters
222d_belt = 0.01 #cross-section height and width [m]
223csa_belt = d_belt**2 #cross-section area [m^2]
224E_belt = 2e9 #E modulus [N/m^2]
225rho_belt = 1e3 #density [kg/m^3]
226F_p = 1e4 #preload force[N]
227
228contactStiffness=1e5 #contactStiffness between the belt and the disks ->
229 #holds the belt around disks
230contactDamping=200 #damping between the belt and the disks
231mu = 0.9 #friction coefficent between the belt and the disks
232velPenalty = 5e2 #frictionVelocityPenalty between tangential velocities
233 #of the belt against the disks
234controller_set_vel = 100 # Desired angular velocity in rad/s of the disk0: 100 rad/s = ~ 955 U/min
235controller_p = 30.0 # P-Factor of controller
236
237MassPerLength_belt = rho_belt * csa_belt
238bendStiffness_belt = E_belt * (d_belt**4 / 12)
239#Hook: F = E * A * epsilon
240axialStiffness_belt = E_belt * csa_belt
241
242#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243#mass parameters
244m_disk0 = 0.5 #mass of disk0 [kg]
245m_disk1 = 1 #mass of disk1 [kg]
246m_crank = 0.5 #mass of crank [kg]
247density_conrod = 1000 #Density of the conrod [kg/m^3]
248m_slider = 0.2 #mass of slider [kg]
249
250#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
251#inertia parameters
252#inertia disks
253J_disk0 = 0.5 * m_disk0 * r_0**2 #inertia disk0 [Nm^2]
254J_disk1 = 0.5 * m_disk1 * r_1**2 #inertia disk1 [Nm^2]
255
256#inertia crank
257J_xx_crank = 5e-3 #inertia xx [Nm^2]
258J_yy_crank = 7e-3 #inertia yy [Nm^2]
259J_zz_crank = 2.5e-3 #inertia zz [Nm^2]
260inertia_crank = RigidBodyInertia(mass=m_crank,
261 inertiaTensor=np.diag([J_xx_crank,
262 J_yy_crank,
263 J_zz_crank]))
264#inertia conrod
265inertia_conrod = InertiaCuboid(density=density_conrod,
266 sideLengths = [L_B,h_B,h_B])
267
268#inertia slider
269#Dummy inertia because MassPoint will result in System Jacobian not invertible!
270inertia_slider = InertiaCuboid(density=(m_slider/h_S**3),
271 sideLengths = [h_S]*3)
272
273#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
274#bearing parameters
275k = 4e4 #stiffness of bearing [N/m]
276d = 8e2 #damping of bearing [N/ms]
277#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
278#load parameters
279M = 0.1 #torque [Nm]
280g = [0,-9.81,0] #acceleration of earth [m/s^2]
281load_crank = [0,-9.81*m_crank,0] #gravity [N]
282load_conrod = [0,-9.81*inertia_conrod.mass,0] #gravity [N]
283load_slider = [0,-9.81*m_slider,0] #gravity [N]
284
285#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
286#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
287
288SC = exu.SystemContainer()
289mbs = SC.AddSystem()
290
291#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
292#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
293#BELT DRIVE SYSTEM
294if sys_set == 0 or sys_set > 2:
295
296 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
297 #grounds
298
299 #ground-node and ground-marker in center of disk0
300 nG_disk0 = mbs.AddNode(NodePointGround(referenceCoordinates = [L_0,0,0]))
301 mG_disk0 = mbs.AddMarker(MarkerNodePosition(nodeNumber = nG_disk0))
302
303 #ground-node and ground-marker in center of disk1
304 nG_disk1 = mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
305 mG_disk1 = mbs.AddMarker(MarkerNodePosition(nodeNumber = nG_disk1))
306
307 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
308 #disk0
309
310 #initial coordinates disk0
311 u0_disk0 = [L_0,0,0]
312 #NodeRigidBody2D on center of gravity disk0
313 nRB2D_disk0 = mbs.AddNode(NodeRigidBody2D(referenceCoordinates = [0,0,0],
314 initialCoordinates = u0_disk0,
315 initialVelocities = [0,0,0]))
316 #visualization of disk0
317 bodyVis_disk0 = VObjectRigidBody2D(graphicsData=[{'type':'Circle',
318 'color':[0,0,0,1],
319 'radius':r_0,
320 'position':[0,0,0],
321 'normal':[0,0,1]},
322 {'type':'Line',
323 'color':[1,0,0,1],
324 'data':[0,0,0,r_0,0,0]}])
325 #ObjectRigidBody2D disk0
326 oRB2D_disk0 = mbs.AddObject(ObjectRigidBody2D(nodeNumber = nRB2D_disk0,
327 physicsMass=m_disk0,
328 physicsInertia=J_disk0,
329 visualization = bodyVis_disk0))
330 #MarkerBodyRigid on center of disk0
331 mNP_disk0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB2D_disk0,
332 localPosition=[0,0,0]))
333
334 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
335 #disk1
336
337 #initial coordinates disk1
338 u0_disk1 = [0,0,0]
339 #NodeRigidBody2D on center of disk1
340 nRB2D_disk1 = mbs.AddNode(NodeRigidBody2D(referenceCoordinates = [0,0,0],
341 initialCoordinates = u0_disk1,
342 initialVelocities = [0,0,0]))
343 #visualization of disk1
344 bodyVis_disk1 = VObjectRigidBody2D(graphicsData=[{'type':'Circle',
345 'color':[0,0,0,1],
346 'radius':r_1,
347 'position':[0,0,0],
348 'normal':[0,0,1]},
349 {'type':'Line',
350 'color':[1,0,0,1],
351 'data':[0,0,0,r_1,0,0]}])
352 #ObjectRigidBody2D disk1
353 oRB2D_disk1 = mbs.AddObject(ObjectRigidBody2D(nodeNumber = nRB2D_disk1,
354 physicsMass=m_disk1,
355 physicsInertia=J_disk1,
356 visualization = bodyVis_disk1))
357 #MarkerBodyRigid on center of disk1
358 mNP_disk1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB2D_disk1,
359 localPosition=[0,0,0]))
360
361 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
362 #joints
363
364 #visualization of joint0 in the center of disk0
365 jointVis_disk0 = VObjectJointRevolute2D(show=True, drawSize=0.02,
366 color=[0,0,0,1])
367 #ObjectJointRevolute2D between ground and disk0
368 oJR2D_disk0 = mbs.AddObject(ObjectJointRevolute2D(markerNumbers = [mG_disk0,mNP_disk0],
369 visualization = jointVis_disk0))
370 #visualization of joint1 in the center of disk1
371 jointVis_disk1 = VObjectJointRevolute2D(show=True, drawSize=0.02, color=[0,0,0,1])
372 #ObjectJointRevolute2D between ground and disk1
373 oJR2D_disk1 = mbs.AddObject(ObjectJointRevolute2D(markerNumbers = [mG_disk1,mNP_disk1],
374 visualization = jointVis_disk1))
375
376 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
377 #function to get length of the belt
378 def belt_get_lengths(L_0, r_l, r_r):
379 alpha = np.arcsin((r_l-r_r)/L_0) #angle of the arc length
380 b_belt = L_0*np.cos(alpha) #branch between the disks
381 al_dr_belt = r_r*(np.pi-2*alpha) #arc length disk0
382 al_dl_belt = r_l*(np.pi+2*alpha) #arc length disk1
383 len_belt = 2*b_belt+al_dr_belt+al_dl_belt #belt length
384 return [alpha, b_belt, al_dl_belt, al_dr_belt, len_belt]
385
386 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
387 #curve parameterization of belt curve
388 def belt_function(p, L_0, r_l, r_r):
389 """
390 Calcultes the x,y-Position and the tangential Vectors of the belt
391 curve at the length parameter p [0-1] given.
392 p = 0 is on up on the left disk where the belt leaves the disk.
393 """
394 [alpha, b_belt, al_dl_belt, al_dr_belt, len_belt] = belt_get_lengths(L_0, r_l, r_r)
395
396 # scale [0-1] to [0 - len_belt]
397 p = p * len_belt
398
399 if p < 0.0 or p > len_belt:
400 return False;
401 elif p < b_belt:
402 element_type = "s_u" #straight up
403 # straight branch up
404 p_element = p
405
406 # calculate start point (p = 0):
407 x_offset = r_l*np.sin(alpha)
408 y_offset = r_l*np.cos(alpha)
409
410 x_slopex = np.cos(alpha)
411 y_slopex = -np.sin(alpha)
412
413 x_pos = x_offset + p_element * x_slopex
414 y_pos = y_offset + p_element * y_slopex
415 elif p < (b_belt + al_dr_belt):
416 element_type = "d_r" #disk right
417 # arc at right disk:
418 p_element = p - b_belt
419
420 alpha_element = p_element / r_r
421
422 # calculate start point at arc:
423 alpha_offset = alpha
424
425 alpha_i = alpha_offset + alpha_element
426 x_pos = L_0 + r_r*np.sin(alpha_i)
427 y_pos = r_r*np.cos(alpha_i)
428 x_slopex = np.cos(alpha_i)
429 y_slopex = -np.sin(alpha_i)
430
431 elif p <= (2 * b_belt + al_dr_belt):
432 element_type = "s_d" #straight down
433 # straight branch down
434 p_element = p - (b_belt + al_dr_belt)
435
436 # calculate start point (p = 0):
437 x_offset = L_0 + r_r*np.sin(alpha)
438 y_offset = -r_r*np.cos(alpha)
439
440 x_slopex = -np.cos(alpha)
441 y_slopex = -np.sin(alpha)
442
443 x_pos = x_offset + p_element * x_slopex
444 y_pos = y_offset + p_element * y_slopex
445
446 elif p <= len_belt:
447 element_type = "d_l" #disk left
448 # arc at left disk:
449 p_element = p - (2 * b_belt + al_dr_belt)
450
451 alpha_element = p_element / r_l
452
453 # calculate start point at arc:q
454 alpha_offset = 2*np.pi - alpha
455
456 alpha_i = alpha_offset + alpha_element
457 x_pos = -r_l*np.sin(alpha_i)
458 y_pos = -r_l*np.cos(alpha_i)
459 x_slopex = -np.cos(alpha_i)
460 y_slopex = np.sin(alpha_i)
461
462 return [x_pos, y_pos, x_slopex, y_slopex, element_type]
463
464 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
465 #Test belt_function
466 if test_belt_function:
467 x = []
468 y = []
469 for pi in range(0,1001):
470 # calculate position
471 p = pi / 1000.0
472 [xi, yi, sxi, syi, el_typei] = belt_function(p, L_0, r_1, r_0)
473 x.append(xi)
474 y.append(yi)
475 np.savetxt("belt.txt", np.array([x,y]).T, delimiter=',')
476
477 slope_x = []
478 slope_y = []
479 for pi in range(n):
480 # calculate position
481 p = pi / n
482 [xi, yi, sxi, syi, el_typei] = belt_function(p, L_0, r_1, r_0)
483 slope_x.append([xi, sxi])
484 slope_y.append([yi, syi])
485 np.savetxt("belt_slope.txt", np.hstack((slope_x,slope_y)),delimiter=',')
486
487 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
488 #belt nodes
489 nNP2DS = []
490 for i in range(n):
491 p = i / n
492 [xi, yi, sxi, syi, el_typei] = belt_function(p, L_0, r_1, r_0)
493 nNP2DS.append(mbs.AddNode(NodePoint2DSlope1(name=el_typei+"_"+str(i),
494 referenceCoordinates = [xi,yi,sxi,syi])))
495
496 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
497 #belt cable elements
498 [alpha, b_belt, al_dl_belt, al_dr_belt, len_belt] = belt_get_lengths(L_0, r_1, r_0)
499 #calculate belt length with preload force:
500 epsilon_belt = -F_p / axialStiffness_belt
501 len_belt_p = len_belt * (1.0 + epsilon_belt)
502
503 print("Belt length: ", len_belt, " and after preload force: ", len_belt_p)
504
505 el_len = len_belt / len(nNP2DS)
506 oANCFC2D_b = []
507
508 if not enable_force:
509 epsilon_belt = 0.0 # Without preload force
510
511 print("epsilon: ", epsilon_belt)
512
513 #ANCF cable objects
514 b_len = 0.0
515 for i in range(len(nNP2DS)):
516 b_len += el_len # Sum all elements -> after loop len_belt must be equal to this sum
517 oANCFC2D_b.append(mbs.AddObject(ObjectANCFCable2D(nodeNumbers=[nNP2DS[i], nNP2DS[(i+1)%len(nNP2DS)]],
518 physicsLength=el_len,
519 physicsBendingStiffness=bendStiffness_belt,
520 physicsMassPerLength=MassPerLength_belt,
521 physicsAxialStiffness=axialStiffness_belt,
522 physicsReferenceAxialStrain=epsilon_belt)))
523
524 print("Belt length from element length sum: ", b_len)
525
526 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
527 #contact friction belt
528
529 if enable_disk_friction:
530
531 nCableGenData_disk0 = []
532 contact_disk0 = []
533 nCableGenData_disk1 = []
534 contact_disk1 = []
535 mCableShape = []
536
537 for i in range(0,len(oANCFC2D_b)):
538 nCS = 8
539 #NodeGenericData disk0
540 nCableGenData_disk0.append(mbs.AddNode(NodeGenericData(initialCoordinates = [0,0,0]*nCS,
541 numberOfDataCoordinates=3*nCS)))
542 #NodeGenericData disk1
543 nCableGenData_disk1.append(mbs.AddNode(NodeGenericData(initialCoordinates = [0,0,0]*nCS,
544 numberOfDataCoordinates=3*nCS)))
545 #MarkerBodyCable2DShape on cable
546 mCableShape.append(mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=oANCFC2D_b[i],
547 numberOfSegments=nCS)))
548 #contact friction to disk0
549 contact_disk0.append(mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mNP_disk0,mCableShape[-1]],
550 nodeNumber=nCableGenData_disk0[-1],
551 circleRadius=r_0,
552 contactStiffness=contactStiffness,
553 contactDamping=contactDamping,
554 numberOfContactSegments=nCS,
555 frictionCoefficient=mu,
556 frictionVelocityPenalty=velPenalty)))
557 #contact friction to disk1
558 contact_disk1.append(mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mNP_disk1, mCableShape[-1]],
559 nodeNumber=nCableGenData_disk1[-1],
560 circleRadius=r_1,
561 contactStiffness=contactStiffness,
562 contactDamping=contactDamping,
563 numberOfContactSegments=nCS,
564 frictionCoefficient=mu,
565 frictionVelocityPenalty=velPenalty)))
566
567 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
568 # Velocity controller
569
570 s_disk0 = mbs.AddSensor(SensorBody(bodyNumber=oRB2D_disk0, writeToFile=True,
571 fileName="angular_velocity_disk0.txt",
572 outputVariableType=exu.OutputVariableType.AngularVelocity))
573
574 def p_controller(mbs, t, loadVector):
575 vel = mbs.GetSensorValues(s_disk0)[2]
576 torque = controller_p * (controller_set_vel - vel)
577 return [0,0,torque]
578
579 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
580 #torque on disk0
581 if enable_controller:
582 l_Torquedisk0 = mbs.AddLoad(LoadTorqueVector(markerNumber=mNP_disk0,
583 loadVectorUserFunction=p_controller))
584
585 else:
586 l_Torquedisk0 = mbs.AddLoad(LoadTorqueVector(markerNumber=mNP_disk0,loadVector=[0,0,M]))
587
588 s_disk1 = mbs.AddSensor(SensorBody(bodyNumber=oRB2D_disk1, writeToFile=True,
589 fileName="angular_velocity_disk1.txt",
590 outputVariableType=exu.OutputVariableType.AngularVelocity))
591 s_load = mbs.AddSensor(SensorLoad(loadNumber=l_Torquedisk0, writeToFile=True,
592 fileName="torque.txt"))
593
594
595#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
596#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
597#CRANK DRIVE SYSTEM 2D
598if sys_set == 1 or sys_set == 3:
599
600 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
601 #ground
602 nPG = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
603 oG = mbs.AddObject(ObjectGround(referencePosition=[0,0,0]))
604 mBPG = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oG, localPosition=[0,0,0]))
605 mNCG = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nPG, coordinate=0))
606
607 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
608 #crank
609
610 #visualization of crank
611 graphics_crank = graphics.RigidLink(p0=[-0.5*L_A,0,-h_A/2],
612 p1=[0.5*L_A ,0,-h_A/2],
613 axis0=[0,0,1], axis1=[0,0,1],
614 radius=[0.5*h_A,0.5*h_A],
615 thickness=h_A, width=[h_A,h_A],
616 color=[0.8,0.8,0.8,1.],nTiles=16)
617 #node on center of gravity of crank
618 nRB2D_crank = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[-L_A*0.5,0,0],
619 initialVelocities=[0,0,0]))
620 #RigidBody2D crank
621 oRB2D_crank = mbs.AddObject(RigidBody2D(physicsMass=m_crank,
622 physicsInertia=J_zz_crank,
623 nodeNumber=nRB2D_crank,
624 visualization=VObjectRigidBody2D(graphicsData=[graphics_crank])))
625
626 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
627 #connecting rod
628
629 #visualization of connecting rod
630 graphics_conrod = graphics.RigidLink(p0=[-0.5*L_B,0,h_B/2],
631 p1=[0.5*L_B,0,h_B/2],
632 axis0=[0,0,1], axis1=[0,0,1],
633 radius=[0.5*h_B,0.5*h_B],
634 thickness=h_B, width=[h_B,h_B],
635 color=[0.7,0.7,0.7,1],nTiles=36)
636 #node on center of gravity of connecting rod
637 nRB2D_conrod = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[-(L_A+L_B*0.5),0,0],
638 initialVelocities=[0,0,0]))
639 #RigidBody2D connecting rod
640 oRB2D_conrod = mbs.AddObject(RigidBody2D(physicsMass=inertia_conrod.mass,
641 physicsInertia=inertia_conrod.inertiaTensor[1,1],
642 nodeNumber=nRB2D_conrod,
643 visualization=VObjectRigidBody2D(graphicsData= [graphics_conrod])))
644
645 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
646 #slider
647
648 #visualization of slider
649 c=0.025 #dimension of slider
650 graphics_slider = graphics.BrickXYZ(-c,-c,-c*2,c,c,0,
651 color=[0.2,0.2,0.2,0.9])
652 #node on center of gravity of slider
653 nP2D_slider = mbs.AddNode(NodePoint2D(referenceCoordinates=[-(L_A+L_B),0]))
654 #MassPoint2D slider
655 oMP2D_slider = mbs.AddObject(MassPoint2D(physicsMass=m_slider,
656 nodeNumber=nP2D_slider,
657 visualization=VObjectRigidBody2D(graphicsData= [graphics_slider])))
658
659 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
660 #markers for joints
661 mBPLeft_crank = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRB2D_crank,
662 localPosition=[-L_A*0.5,0.,0.]))
663 mBPRight_crank = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRB2D_crank,
664 localPosition=[L_A*0.5,0.,0.]))
665 mBPLeft_conrod = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRB2D_conrod,
666 localPosition=[-L_B*0.5,0.,0.]))
667 mBPRight_conrod = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRB2D_conrod,
668 localPosition=[L_B*0.5,0.,0.]))
669 mBP_slider = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMP2D_slider,
670 localPosition=[ 0.,0.,0.]))
671
672 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
673 #joints
674 oRJ2D_ground_crank = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mBPG,mBPRight_crank]))
675 oRJ2D_crank_conrod = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mBPLeft_crank,mBPRight_conrod]))
676 oRJ2D_slider = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mBP_slider,mBPLeft_conrod]))
677
678 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
679 #markers for node constraints
680 mNC_Y_slider = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nP2D_slider, coordinate=1)) #y-coordinate is constrained
681
682 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
683 #coordinate constraints for slider (free motion in x-direction)
684 oCC_ground_slider = mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCG, mNC_Y_slider]))
685
686 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
687 #markers for load
688 mBR_crank_torque = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB2D_crank, localPosition=[L_A/2,0.,0.]))
689 mBR_crank_gravity = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB2D_crank, localPosition=[0.,0.,0.]))
690 mBR_conrod_gravity = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB2D_conrod, localPosition=[0.,0.,0.]))
691
692 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
693 #loads and driving forces
694
695 #gravity of crank
696 lC_crank_gravity = mbs.AddLoad(Force(markerNumber = mBR_crank_gravity,
697 loadVector = load_crank))
698 #gravity of conrod
699 lC_conrod_gravity = mbs.AddLoad(Force(markerNumber = mBR_conrod_gravity,
700 loadVector = load_conrod))
701 #gravity of slider
702 lC_slider_gravity = mbs.AddLoad(Force(markerNumber = mBP_slider,
703 loadVector = load_slider))
704
705 if sys_set == 1:
706 #torque at crank
707 lC_crank_torque = mbs.AddLoad(Torque(markerNumber = mBR_crank_torque,
708 loadVector = [0, 0, M]))
709
710 if sys_set == 3:
711 #ConnectorCoordinate - crank gets torque of disk1
712 mNC_disk1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRB2D_disk1,coordinate=2))
713 mNC_crank = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRB2D_crank,coordinate=2))
714 oCC_disk1_crank = mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mNC_disk1,mNC_crank],
715 velocityLevel=True))
716
717#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
718#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
719#CRANK DRIVE SYSTEM 3D
720if sys_set == 2 or sys_set == 4:
721
722 nodeType=exu.NodeType.RotationEulerParameters
723
724 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
725 #ground
726 oG_Left = mbs.AddObject(ObjectGround())
727 oG_Right = mbs.AddObject(ObjectGround())
728 oG_slider = mbs.AddObject(ObjectGround())
729
730 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
731 #crank
732 graphics_crank_1 = graphics.RigidLink(p0=[0,0,0],p1=[0,0,-b_a0], axis1=[0,0,1], radius=[h_A/2,h_A/1.3],
733 thickness = h_A, width=[0,h_A/2], color=[0.8,0.8,0.8,1])
734 graphics_crank_2 = graphics.RigidLink(p0=[0,0,0],p1=[-L_A,0,0], radius=[h_A/2,h_A/2],
735 thickness = h_A, color=[0.8,0.8,0.8,1])
736 graphics_crank_3 = graphics.RigidLink(p0=[-L_A,0,0],p1=[-L_A,0,b_a1], radius=[h_A/2,h_A/2],
737 thickness = h_A, color=[0.8,0.8,0.8,1])
738 graphics_crank_4 = graphics.RigidLink(p0=[-L_A,0,b_a1],p1=[0,0,b_a1], radius=[h_A/2,h_A/2],
739 thickness = h_A, color=[0.8,0.8,0.8,1])
740 graphics_crank_5 = graphics.RigidLink(p0=[0,0,b_a1],p1=[0,0,b_a1+b_a2],axis1=[0,0,1], radius=[h_A/2,h_A/1.3],
741 thickness = h_A, width=[0,h_A/2], color=[0.8,0.8,0.8,1])
742 [nRG_crank,oRB_crank]=AddRigidBody(mainSys = mbs, inertia=inertia_crank, nodeType=str(nodeType),
743 position=[0,0,b_a0], angularVelocity=[0,0,0],
744 gravity=g, graphicsDataList=[graphics_crank_1,
745 graphics_crank_2,graphics_crank_3, graphics_crank_4,graphics_crank_5])
746
747 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
748 #connecting rod
749 graphics_conrod = graphics.RigidLink(p0=[L_B/2,0,0],p1=[-L_B/2,0,0], axis0=[0,0,1], axis1=[0,0,1], radius=[h_B/1.5,h_B/2],
750 thickness = h_B, width=[h_B,h_B], color=[0.5,0.5,0.5,1])
751 [nRG_conrod,oRB_conrod]=AddRigidBody(mainSys = mbs, inertia=inertia_conrod, nodeType=str(nodeType), angularVelocity=[0,0,0],
752 position=[-L_A-L_B/2,0,b_a0+b_a1/2], gravity=g, graphicsDataList=[graphics_conrod])
753
754 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
755 #slider
756 d=0.07
757 graphics_slider = graphics.BrickXYZ(-d/2,-d/2,-d-h_B/2,d/2,d/2,-h_B/2,
758 color=[0.2,0.2,0.2,0.9])
759 [nRB_slider,oRB_slider]=AddRigidBody(mainSys = mbs,
760 inertia=inertia_slider,
761 nodeType=str(nodeType),
762 angularVelocity=[0,0,0],
763 position=[-(L_A+L_B),0,b_a0+b_a1/2],
764 gravity=g,
765 graphicsDataList=[graphics_slider])
766
767 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
768 #marker for joints
769 mG_Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oG_Left,
770 localPosition=[0,0,0]))
771 mG_Right = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oG_Right,
772 localPosition=[0,0,b_a0+b_a1+b_a2]))
773 mG_slider = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oG_slider,
774 localPosition=[-(L_A+L_B),0,b_a0+b_a1/2]))
775 mBR_crank_Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_crank,
776 localPosition=[0,0,-b_a0]))
777 mBR_crank_Right = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_crank,
778 localPosition=[0,0,b_a1+b_a2]))
779 mBR_crank_conrod = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_crank,
780 localPosition=[-L_A,0,b_a1/2]))
781 mBR_conrod_crank = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_conrod,
782 localPosition=[L_B/2,0,0]))
783 mBR_conrod_slider = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_conrod,
784 localPosition=[-L_B/2,0,0]))
785 mBR_slider_conrod = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_slider,
786 localPosition=[0,0,0]))
787 mBR_slider = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_slider,
788 localPosition=[0,0,0]))
789
790 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++q
791 #joints
792 oGJ_crank_Left = mbs.AddObject(ObjectConnectorCartesianSpringDamper(markerNumbers=[mG_Left, mBR_crank_Left],
793 stiffness=[k]*3, damping=[d]*3))
794 oGJ_crank_Right = mbs.AddObject(ObjectConnectorCartesianSpringDamper(markerNumbers=[mG_Right, mBR_crank_Right],
795 stiffness=[k]*3, damping=[d]*3))
796 oRJ2D_crank_conrod = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mBR_crank_conrod,mBR_conrod_crank],
797 visualization=VObjectJointRevolute2D(drawSize=0.05)))
798 oRJ2D_conrod_slider = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mBR_conrod_slider,mBR_slider_conrod],
799 visualization=VObjectJointRevolute2D(drawSize=0.05)))
800 oJP2D_slider = mbs.AddObject(ObjectJointPrismatic2D(markerNumbers=[mG_slider,mBR_slider],
801 axisMarker0 = [1.,0.,0.],
802 normalMarker1 = [0.,1.,0.],
803 visualization=VObjectJointPrismatic2D(drawSize=0.01)))
804 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
805 #load and driving forces
806
807 if sys_set == 2:
808 #markers for load
809 mBR_crank_torque = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRB_crank,
810 localPosition=[0.,0.,b_a1+b_a2]))
811 #driving forces
812 lC_crank_torque = mbs.AddLoad(Torque(markerNumber = mBR_crank_torque,
813 loadVector = [0, 0, M/2]))
814
815 if sys_set == 4:
816 #ConnectorCoordinate - crank gets torque of disk1
817 #markers for Connector
818 mNC_disk1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRB2D_disk1,coordinate=2))
819 # disk1 is a NodeRigidBody2D -> coordinates = [q_0,q_1,psi] -> Rotation psi
820 mNC_crank = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nRG_crank,rotationCoordinate=2))
821 # crank is a 3D Node -> Rotation around z-axis
822 #Connector Coordinate
823 oCC_disk1_crank = mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mNC_disk1,mNC_crank],
824 velocityLevel=True))
825
826 s_crank = mbs.AddSensor(SensorBody(bodyNumber=oRB_crank, writeToFile=True,
827 fileName="crank_pos.txt", localPosition=[0,0,-b_a0],
828 outputVariableType=exu.OutputVariableType.Position))
829
830
831#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
832#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
833
834mbs.Assemble()
835print(mbs)
836
837#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
838#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
839#simulation
840
841if sys_set == 0:
842 tEnd = 1 #end time of the simulation
843 steps = 1000 #number of steps
844elif sys_set > 0:
845 tEnd = 2 #end time of the simulation
846 steps = 1000 #number of steps
847
848sims = exu.SimulationSettings()
849sims.timeIntegration.generalizedAlpha.spectralRadius=1
850if export_images:
851 sims.solutionSettings.recordImagesInterval = 0.001
852 if not os.path.isdir("images"):
853 os.mkdir("images")
854else:
855 sims.solutionSettings.recordImagesInterval = -1
856sims.pauseAfterEachStep = False
857sims.displayStatistics = True
858sims.displayComputationTime = True
859sims.timeIntegration.numberOfSteps = steps
860sims.timeIntegration.endTime = tEnd
861sims.solutionSettings.sensorsWritePeriod = 0.001
862
863#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
864#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
865#visualizaion
866SC.visualizationSettings.general.circleTiling=128
867dSize = 0.02
868SC.visualizationSettings.nodes.defaultSize = dSize
869SC.visualizationSettings.markers.defaultSize = dSize
870SC.visualizationSettings.bodies.defaultSize = [dSize, dSize, dSize]
871SC.visualizationSettings.connectors.defaultSize = dSize
872SC.visualizationSettings.nodes.show=True
873SC.visualizationSettings.loads.show=True
874SC.visualizationSettings.markers.show=True
875SC.visualizationSettings.sensors.show=True
876SC.visualizationSettings.general.autoFitScene=False
877if sys_set == 0:
878 SC.visualizationSettings.openGL.initialCenterPoint = [L_0/2,0,0]
879 SC.visualizationSettings.openGL.initialZoom = 0.5
880elif sys_set == 1:
881 SC.visualizationSettings.openGL.initialCenterPoint = [-L_A,0,0]
882 SC.visualizationSettings.openGL.initialZoom = 0.4
883elif sys_set == 2:
884 SC.visualizationSettings.openGL.initialCenterPoint = [-0.1,0,0]
885 SC.visualizationSettings.openGL.initialZoom = 0.3
886elif sys_set == 3:
887 SC.visualizationSettings.openGL.initialCenterPoint = [0,0,0]
888 SC.visualizationSettings.openGL.initialZoom = 0.5
889elif sys_set == 4:
890 SC.visualizationSettings.openGL.initialCenterPoint = [0,0,0]
891 SC.visualizationSettings.openGL.initialZoom = 0.5
892SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Force
893
894if sys_set == 0:
895 rot_alpha=0
896 rot_beta=0
897 rot_gamma=0
898elif sys_set > 0:
899 rot_alpha=-0.6
900 rot_beta=0.6
901 rot_gamma=0.37
902
903R_x = np.array([[ 1, 0, 0],
904 [ 0, np.cos(rot_alpha), -np.sin(rot_alpha)],
905 [ 0, np.sin(rot_alpha), np.cos(rot_alpha)]])
906R_y = np.array([[ np.cos(rot_beta), 0, np.sin(rot_beta)],
907 [ 0, 1, 0],
908 [ -np.sin(rot_beta), 0, np.cos(rot_beta)]])
909R_z = np.array([[ np.cos(rot_gamma), -np.sin(rot_gamma), 0],
910 [ np.sin(rot_gamma), np.cos(rot_gamma), 0],
911 [ 0, 0, 1]])
912IMR = np.dot(R_x,R_y)
913IMR = np.dot(IMR,R_z)
914SC.visualizationSettings.openGL.initialModelRotation = [[IMR[0,0],IMR[0,1],IMR[0,2]],
915 [IMR[1,0],IMR[1,1],IMR[1,2]],
916 [IMR[2,0],IMR[2,1],IMR[2,2]]]
917
918#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
919#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
920#Rendering
921exu.StartRenderer() #start graphics visualization
922mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
923mbs.SolveDynamic(sims)
924SC.WaitForRenderEngineStopFlag() #wait for pressing 'Q' to quit
925exu.StopRenderer() #safely close rendering window!
926
927#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
928#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
929
930if export_images:
931 visHelperPlot_all()