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