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