sliderCrank3DwithANCFbeltDrive.py
You can view and download this file on Github: sliderCrank3DwithANCFbeltDrive.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: David Wibmer and Dominik Sponring
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#note: tested with PYTHON version = 3.6.10 and EXUDYN version = 0.1.342
15
16import exudyn as exu
17#a local run must not open a plot window when windows are suppressed (#2477)
18if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
19from exudyn.itemInterface import*
20from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
21import exudyn.graphics as graphics
22
23import numpy as np
24
25
26#Print EXUDYN version
27print('EXUDYN version='+exu.config.Version())
28
29
30#Paramters
31#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33#General
34tSim = 2 #Simulation time + time to fasten belt
35duration = 0.25 #Time for fasten up belt
36
37omega = 80 #Angular velocity to controll
38sweep = True #Controll a sweep of omega
39
40gravitation = False
41
42#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
43##Crank-System##
44
45#Parameters crank
46L_A = 0.15 #[m]
47
48ba_0 = 0.1 #[m]
49ba_1 = 0.05 #[m]
50ba_2 = 0.05 #[m]
51
52Jxx = 0.005 #[kg m^2]
53Jyy = 0.007 #[kg m^2]
54Jzz = 0.0025 #[kg m^2]
55
56massCrank = 0.5 #[Kg]
57
58#Parameters Slider
59massSlider = 0.2 #[Kg]
60
61#Parameters connection rod
62density = 7850 #[Kg/m^3]
63h_B = 0.025 #[m]
64L_B = 0.3 #[m]
65V_CR = L_B*h_B**2 #[m^3]
66massCR = V_CR*density #[Kg]
67Jxx_CR = (1/12)*massCR*(h_B**2 + h_B**2) #[kg m^2]
68Jyy_CR = (1/12)*massCR*(h_B**2 + L_B**2) #[kg m^2]
69Jzz_CR = (1/12)*massCR*(h_B**2 + L_B**2) #[kg m^2]
70
71#Parameters bearing
72 ##rigidity has been increased##
73stiffnes = 4e5 #[N/m]
74damping = 8e2 #[N/(ms)]
75#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
76##Belt-System##
77 ##wheelbase has been inreased##
78L0 = 0.51 #wheelbase [m]
79r0 = 0.05 #radius disk 0 [m]
80r1 = 0.1 #radius disk 1 [m]
81a = 0.01 #belt width [m]
82A = a*a #belt cross section [m^2]
83L_eff = 1.5 #effective belt length [m]
84
85#Material parameters
86E_B = 2e9 #modulus of elasticity [N/m^2]
87rho = 1000 #density [kg/m^3]
88I = a*a*a*a/12 #second moment of area of ANCF element in m^4
89mu = 0.7 #friction coefficient []
90
91m_disk0 = 0.5 #mass disk0[kg]
92m_disk1 = 1 #mass disk1[kg]
93
94J_disk0 = m_disk0*r0*r0/2 #Moment of inertia disk0 [kg m^2]
95J_disk1 = m_disk1*r1*r1/2 #Moment of inertia disk1 [kg m^2]
96
97zOff = 0.12 #additional offset, mainly for drawing reasons
98#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
99#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
100#Creat system container
101SC = exu.SystemContainer()
102
103mbs = SC.AddSystem()
104
105#Generate ground
106#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
107nGround = mbs.AddNode(NodePointGround())
108oGround = mbs.AddObject(ObjectGround(referencePosition=[0,0,zOff]))
109mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
110
111
112
113#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
114#Create crank system
115#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
116
117#Generate Visualisation-Objects
118#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
119
120vSlider = graphics.Cylinder([0.05,0,0], [-0.1,0,0],
121 0.05, [1,0,0,1], nTiles=32)
122vRod = graphics.Brick(centerPoint=[0, 0, 0], size=[L_B, h_B, h_B], color=[0,1,0,1])
123
124vCrank0 = graphics.Cylinder([0,0,-2*ba_1], [0,0,0.01],
125 r1+a/2,color=[0.3,0.3,0.9,1], nTiles=64)
126vCrank1 = graphics.Cylinder([0,0,0.01], [0,0,-ba_0-0.01],
127 0.01,color=[0.3,0.3,0.9,1])
128vCrank2 = graphics.Cylinder([0,0,ba_1-0.01], [0,0,ba_2+0.01],
129 0.01, color=[0.3,0.3,0.9,1])
130#vCrank3 = graphics.Cylinder([-L_A/2,0,-0.0125], [0,0,0.025],
131# 0.1, color=[0.3,0.3,0.9,0.9])
132#vCrank4 = graphics.Cylinder([-L_A/2,0,0.0375], [0,0,0.025],
133# 0.1, color=[0.3,0.3,0.9,0.9])
134
135vCrank3 = graphics.Brick([-L_A/2,0,+0.005], [L_A+0.01,0.01,0.008],
136 color=[0.3,0.3,0.9,0.9])
137vCrank4 = graphics.Brick([-L_A/2,0,0.05-0.005], [L_A+0.01,0.01,0.008],
138 color=[0.3,0.3,0.9,0.9])
139
140vCrank5 = graphics.Cylinder([-L_A,0,0.0], [0,0,0.05],
141 0.01,color=[0.3,0.3,0.9,1])
142
143vDisk_line0 = graphics.Lines([[0,-0.001,0], [r0,-0.001,0], [r0,0.001,0], [0,0.001,0], [0,-0.001,0]])
144vDisk_line1 = graphics.Lines([[0,-0.001,0], [r1,-0.001,0], [r1,0.001,0], [0,0.001,0], [0,-0.001,0]])
145
146cylDisc0 = graphics.Cylinder([0,0,-0.005], [0,0,0.01],
147 r0+a/2,color=[0.3,0.3,0.9,1], nTiles=32)
148
149#Generate Nodes and Objects
150#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
151#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
152
153##Crank##
154#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155ep0 = [1,0,0,0] #no rotation
156nCrank_3D = mbs.AddNode(RigidEP(referenceCoordinates=[0,0,-ba_1/2+zOff]+ep0))
157oCrank_3D = mbs.AddObject(RigidBody(mass=massCrank,
158 inertia=[Jxx,Jyy,Jzz,0,0,0],
159 nodeNumber=nCrank_3D,
160 visualization=VObjectRigidBody2D(graphicsData=[vCrank0,
161 vCrank1,
162 vCrank2,
163 vCrank3,
164 vCrank4,vCrank5])))
165
166##Rod##
167#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
168nRod = mbs.AddNode(RigidEP(referenceCoordinates=[-(L_A+L_B/2),0,0+zOff]+ep0));
169oRod = mbs.AddObject(RigidBody(mass=massCR,
170 inertia=[Jxx_CR,Jyy_CR,Jzz_CR,0,0,0],
171 nodeNumber=nRod,
172 visualization=VObjectRigidBody2D(graphicsData=[vRod])))
173
174
175##Slider##
176#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
177nSlider = mbs.AddNode(Point(referenceCoordinates=[-(L_A+L_B), 0,0+zOff]))
178oSlider = mbs.AddObject(MassPoint(mass = massSlider,
179 nodeNumber = nSlider,
180 visualization=VObjectMassPoint(graphicsData= [vSlider])))
181
182
183#Generate Joints
184#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
185#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
186
187##Markers##
188#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
189mCrank_Rod = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oCrank_3D,
190 localPosition=[-L_A,0,ba_1/2]))
191mjoint2leftC = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oCrank_3D,
192 localPosition = [0,0,-ba_0]))
193mjoint2rightC = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oCrank_3D,
194 localPosition = [0,0,ba_1+ba_2]))
195
196mRodLeft = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRod,
197 localPosition=[-L_B/2,0,0]))
198mRodRight = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRod,
199 localPosition=[ L_B/2,0,0]))
200
201mSlider = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oSlider,
202 localPosition=[0,0,0]))
203
204mjoint2leftG = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround,
205 localPosition=[0,0,-(ba_0+ba_1/2)]))
206mjoint2rightG = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround,
207 localPosition=[0,0,ba_2+ba_1/2]))
208
209##Joints##
210#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
211mbs.AddObject(SphericalJoint(markerNumbers = [mjoint2leftG,mjoint2leftC],
212 constrainedAxes = [0,0,0],
213 visualization = VObjectJointSpherical(jointRadius= 0.01)))
214mbs.AddObject(SphericalJoint(markerNumbers = [mjoint2rightG,mjoint2rightC],
215 constrainedAxes = [0,0,0],
216 visualization = VObjectJointSpherical(jointRadius= 0.01)))
217mbs.AddObject(SphericalJoint(markerNumbers = [mCrank_Rod,mRodRight],
218 constrainedAxes = [1,1,1],
219 visualization = VObjectJointSpherical(jointRadius= 0.01)))
220mbs.AddObject(SphericalJoint(markerNumbers = [mRodLeft,mSlider],
221 constrainedAxes = [1,1,1],
222 visualization = VObjectJointSpherical(jointRadius= 0.01)))
223
224
225##Joints for bearing##
226#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
227mbs.AddObject(CartesianSpringDamper(markerNumbers=[mjoint2leftG,mjoint2leftC],
228 stiffness=[stiffnes,stiffnes,stiffnes],
229 damping=[damping,damping,damping]))
230mbs.AddObject(CartesianSpringDamper(markerNumbers=[mjoint2rightG,mjoint2rightC],
231 stiffness=[stiffnes,stiffnes,stiffnes],
232 damping=[damping,damping,damping]))
233
234
235##Constrains##
236#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
237#Y-coordinate is constrained
238mSliderY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlider,
239 coordinate=1))
240#Z-coordinate is constrained
241mSliderZ = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlider,
242 coordinate=2))
243
244mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mSliderY]))
245mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mSliderZ]))
246
247
248
249
250#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
251#Create Belt system
252#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
253
254
255#Generate disks
256#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
257#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
258#Disk 1 (Driven Disk)
259nDisk1 = mbs.AddNode(Rigid2D(referenceCoordinates=[0, 0, 0]))
260oDisk1 = mbs.AddObject(RigidBody2D(mass=m_disk1, inertia=J_disk1,
261 nodeNumber=nDisk1,
262 visualization=VObjectRigidBody2D(graphicsData=[vDisk_line1])))
263mDisk1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oDisk1))
264
265#Disk 0 (Driver Disk)
266#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
267nDisk0 = mbs.AddNode(Rigid2D(referenceCoordinates=[0, 0, 0]))
268oDisk0 = mbs.AddObject(RigidBody2D(mass=m_disk0, inertia=J_disk0,
269 nodeNumber=nDisk0,
270 visualization=VObjectRigidBody2D(graphicsData=[vDisk_line0,cylDisc0])))
271mDisk0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oDisk0))
272#Marker for contactCable
273#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
274mBDisk0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oDisk0))
275#Marker for coordinateConstraint
276#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
277mNDisk00 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDisk0, coordinate=0))
278mNDisk01 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDisk0, coordinate=1))
279
280
281#Generate disc joints
282#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
283def OffsetUF(mbs, t, itemIndex, lOffset):
284 if t < duration: return L0*(1-np.cos(t*np.pi/duration))/2
285 else: return L0
286
287#Joint for Disk1 ->Rigid
288#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
289mpoint1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround))
290mbs.AddObject(RevoluteJoint2D(markerNumbers=[mpoint1, mDisk1]))
291#Joint for Disk0 ->tighten the belt
292#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
293mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround, mNDisk01]))
294oTighten = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround, mNDisk00],
295 offsetUserFunction=OffsetUF))
296
297
298#Generate belt
299#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
300#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
301cableList = []
302nodeList = []
303markerList = []
304
305#Calculate element parameters
306#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
307beltRadius = L_eff/(2*np.pi)
308nElements = 30*2
309angleSegments = 2*np.pi/nElements
310arcLenght = angleSegments*beltRadius
311
312#Create belt-nodes
313#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
314for i in range(nElements):
315 s = np.sin(angleSegments*i)
316 c = np.cos(angleSegments*i)
317 nodeList += [mbs.AddNode(Point2DS1(referenceCoordinates = [c*beltRadius,s*beltRadius,-s,c]))]
318
319#Create belt-objects
320#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
321for i in range(1, nElements):
322 cableList += [mbs.AddObject(Cable2D(length=arcLenght, massPerLength=rho*A,
323 bendingStiffness=E_B*I, axialStiffness=E_B*A,
324 referenceCurvature=1/beltRadius,
325 nodeNumbers=[nodeList[i-1],nodeList[i]]))]
326
327#Connect first and last belt-object
328#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
329cableList += [mbs.AddObject(Cable2D(length=arcLenght, massPerLength=rho*A,
330 bendingStiffness=E_B*I, axialStiffness=E_B*A,
331 referenceCurvature=1/beltRadius,
332 nodeNumbers=[nodeList[-1],nodeList[0]]))]
333
334
335
336#Add gravity
337#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
338if gravitation:
339 for i in range(len(nodeList)):
340 markerList += [mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeList[i]))]
341 mbs.AddLoad(Force(markerNumber=markerList[i], loadVector=[0, -9.81*rho*A*arcLenght, 0]))
342
343
344#Add contact
345#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
346#Contact parameters
347cStiffness = 1e6
348cDamping = 1000
349
350
351nSegments = 3
352initialGapList = [0.1, 0.1, 0.1]*nSegments
353
354for i in range(len(cableList)):
355 #Generate markers for the ANCF-Elements
356 mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[i], numberOfSegments=nSegments))
357 #Generate contact for disk1
358 nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,
359 numberOfDataCoordinates=3*nSegments))
360 mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mDisk1, mCable],
361 nodeNumber = nodeDataContactCable,
362 contactStiffness = cStiffness,
363 contactDamping = cDamping,
364 frictionVelocityPenalty = 1000,
365 frictionCoefficient = mu,
366 circleRadius = r1+a/2))
367 #Generate contact for disk0
368 nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,
369 numberOfDataCoordinates=3*nSegments))
370 mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mBDisk0, mCable],
371 nodeNumber = nodeDataContactCable,
372 contactStiffness = cStiffness,
373 contactDamping = cDamping,
374 frictionVelocityPenalty = 1000,
375 frictionCoefficient = mu,
376 circleRadius = r0+a/2))
377
378
379#Generate Load
380#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
381def loadUF(mbs, t,loadVector):
382 if t < 2*duration : return [0,0,0]
383 if sweep == True:
384 omega_soll = 10 + omega*(t-2*duration)/tSim
385 else:
386 omega_soll = omega
387
388 omega_ist = mbs.GetNodeOutput(nCrank_3D, exu.OutputVariableType.AngularVelocity)[2]
389
390 tor = (omega_soll - omega_ist)*5
391
392 return [0,0,tor]
393
394
395#apply torque at disk0
396#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
397mbs.AddLoad(Torque(markerNumber=mDisk0,
398 loadVectorUserFunction=loadUF))
399
400mbs.AddObject(GenericJoint(markerNumbers=[mDisk1,mjoint2rightC],
401 constrainedAxes=[0,0,0,0,0,1],
402 visualization=VObjectJointGeneric(show=False)))
403
404
405#Generate Sensors
406#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
407
408#Sensor for measuring disk0
409mbs.AddSensor(SensorObject(objectNumber=oTighten,
410 fileName='solution/Preload_overallS.txt',
411 outputVariableType=exu.OutputVariableType.Force))
412mbs.AddSensor(SensorBody(bodyNumber=oDisk0,
413 fileName='solution/Pos_Disk0_overallS.txt',
414 outputVariableType=exu.OutputVariableType.Position))
415
416
417mbs.AddSensor(SensorNode(nodeNumber=nCrank_3D,
418 fileName='solution/Angular_velocity_overallS.txt',
419 outputVariableType=exu.OutputVariableType.AngularVelocity))
420
421
422#Assamble the system
423#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
424mbs.Assemble()
425
426
427
428#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
429#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
430simulationSettings = exu.SimulationSettings()
431
432simulationSettings.timeIntegration.numberOfSteps = tSim*1000
433simulationSettings.timeIntegration.endTime = tSim
434simulationSettings.solution.file.write = True
435simulationSettings.timeIntegration.verboseMode = 1
436simulationSettings.timeIntegration.newton.relativeTolerance = 1e-10
437
438simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*10 #10000
439simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
440
441simulationSettings.timeIntegration.newton.useModifiedNewton = False
442simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
443simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
444simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 6.055454452393343e-06*0.1 #eps^(1/3)
445simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1e8
446
447simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
448simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6
449simulationSettings.show.statistics = True
450
451
452
453#Visualisation Settings
454#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
455SC.visualizationSettings.nodes.showNumbers = False
456SC.visualizationSettings.bodies.showNumbers = False
457SC.visualizationSettings.connectors.showNumbers = False
458SC.visualizationSettings.connectors.defaultSize = 0.005
459SC.visualizationSettings.contact.contactPointsDefaultSize = 0.01
460SC.visualizationSettings.connectors.showContact = True
461
462#create animation:
463if False:
464 simulationSettings.solution.recordImagesInterval = 0.002
465 SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
466 SC.visualizationSettings.view0.window.renderWindowSize = [1920,1080]
467 SC.visualizationSettings.openGL.multiSampling = 4
468
469
470SC.renderer.Start()
471if 'lastRenderState' in vars():
472 SC.renderer.SetState(lastRenderState) #load last model view
473SC.renderer.DoIdleTasks()
474
475mbs.SolveDynamic(simulationSettings)
476
477
478SC.renderer.DoIdleTasks()
479SC.renderer.Stop() #safely close rendering window!
480lastRenderState = SC.renderer.GetState() #store model view for next simulation
481
482
483
484
485#Show Sensor
486#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
487import matplotlib.pyplot as plt
488import matplotlib.ticker as ticker
489
490#Import sensor data
491data = np.loadtxt(OutputFilePath('solution/Angular_velocity_overallS.txt'), comments='#', delimiter=',')
492plt.figure(1)
493plt.plot(data[:,0], data[:,3], 'r-', label='controlled sweep')
494
495
496ax=plt.gca() # get current axes
497ax.grid(True, 'major', 'both')
498ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
499plt.xlabel('time (s)')
500ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
501plt.ylabel('Angular velocity')
502plt.legend() #show labels as legend
503plt.tight_layout()
504plt.show()
505
506
507#Import sensor data
508data = np.loadtxt(OutputFilePath('solution/Preload_overallS.txt'), comments='#', delimiter=',')
509plt.figure(2)
510# 1.column = time 2.column = load
511plt.plot(data[:,0], data[:,1], 'r-', label='Preload')
512
513
514ax=plt.gca() # get current axes
515ax.grid(True, 'major', 'both')
516ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
517plt.xlabel('time (s)')
518ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
519plt.ylabel('force (N)')
520plt.legend() #show labels as legend
521plt.tight_layout()
522plt.show()
523
524
525
526
527#Import sensor data
528data = np.loadtxt(OutputFilePath('solution/Pos_Disk0_overallS.txt'), comments='#', delimiter=',')
529plt.figure(3)
530# 1.column = time 2.column = x-axis
531plt.plot(data[:,0], data[:,1], 'r-', label='Position X-Axis')
532
533
534ax=plt.gca() # get current axes
535ax.grid(True, 'major', 'both')
536ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
537plt.xlabel('time (s)')
538ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
539plt.ylabel('position (m)')
540plt.legend() #show labels as legend
541plt.tight_layout()
542plt.show()