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