beltDrivesComparison.py
You can view and download this file on Github: beltDrivesComparison.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Comparison of functionality of different belt drives
5#
6# Author: Johannes Gerstmayr
7# Date: 2022-11-05
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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import sys
14sys.exudynFast = True
15
16import exudyn as exu
17from exudyn.itemInterface import *
18from exudyn.utilities import *
19from exudyn.beams import *
20
21import numpy as np
22from math import sin, cos, pi, sqrt , asin, acos, atan2
23import copy
24
25SC = exu.SystemContainer()
26mbs = SC.AddSystem()
27
28
29
30#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31#settings:
32useGraphics= True
33useGeneralContact = False
34useContact = True
35staticEqulibrium = not useGeneralContact
36
37tEnd = 20 #end time of dynamic simulation
38stepSize = 2*1e-3 #step size
39useFriction = True
40dryFriction = 0.35
41nContactSegments = 8
42dEAfact = 2
43
44#GeneralContact:
45contactStiffness = 4e5
46contactDamping = 1e-3*contactStiffness
47
48#CableFriction contact:
49
50if not useGeneralContact:
51 contactStiffness = 4e5*0.1
52 contactDamping = 1e-2*contactStiffness*0.1
53 dEAfact = 1
54 stepSize = 0.25*1e-3 #step size
55 dryFriction = 0.2 #lower because more accurate ...
56 contactFrictionStiffness = contactStiffness*0.2
57 frictionVelocityPenalty = 1e-4*contactFrictionStiffness
58 nContactSegments = 2
59
60wheelMass = 1
61wheelInertia = 0.01
62
63
64#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
65#belt:
66g = 9.81
67gVec = [0,-g,0] # gravity vector
68E=1e9 # Young's modulus of ANCF element in N/m^2
69rhoBeam=1000 # density of ANCF element in kg/m^3
70b=0.002 # width of rectangular ANCF element in m
71h=0.002 # height of rectangular ANCF element in m
72A=b*h # cross sectional area of ANCF element in m^2
73I=b*h**3/12 # second moment of area of ANCF element in m^4
74EI = E*I # bending stiffness
75EA = E*A # axial stiffness
76dEI = 0*1e-3*EI #bending proportional damping
77dEA = dEAfact*1e-2*EA #axial strain proportional damping
78dimZ = 0.02 #z.dimension
79
80nANCFnodes = 50*2
81preStretch=-0.005 #relative stretch in beltdrive
82
83#wheels
84angularVelocity = 2*pi*1
85velocityControl = True
86rWheel = 0.25
87dWheels = 0.6
88leverArm = 0.5
89massArmTip = 1
90
91def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
92 return lOffset*min(angularVelocity*t, angularVelocity)
93
94
95#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96#create reeving system for each belt drive
97#complicated shape:
98circleList0 = [
99 [[0,0],rWheel,'L'],
100 [[dWheels,0],rWheel,'L'],
101 [[0,0],rWheel,'L'],
102 [[dWheels,0],rWheel,'L'],
103 ]
104circleList1 = [
105 [[0,0],rWheel,'R'],
106 [[dWheels,0],rWheel,'L'],
107 [[0,0],rWheel,'R'],
108 [[dWheels,0],rWheel,'L'],
109 ]
110factRadius = 0.5
111circleList2 = [
112 [[0,0],rWheel*factRadius,'L'],
113 [[dWheels,0],rWheel,'L'],
114 [[0,0],rWheel*factRadius,'L'],
115 [[dWheels,0],rWheel,'L'],
116 ]
117tensionerY = rWheel*0.7
118circleList3 = [
119 [[0,0],rWheel*factRadius,'L'],
120 [[dWheels*0.4,tensionerY],rWheel*0.2,'R'],
121 [[dWheels,0],rWheel,'L'],
122 [[dWheels*0.4,-tensionerY],rWheel*0.2,'R'],
123 [[0,0],rWheel*factRadius,'L'],
124 [[dWheels*0.4,tensionerY],rWheel*0.2,'R'],
125 ]
126
127
128reevingSystems = [circleList0,circleList1,circleList2,circleList3]
129# reevingSystems = [circleList0,circleList1,circleList2,]
130# reevingSystems = [circleList3]
131
132contactObjects = [] #ContactFrictionCable objects to modify in static initialization
133velControlList = [] #for static initialization
134fixWheelsList = [] #for static initialization
135
136for cnt, circleList in enumerate(reevingSystems):
137 offX = (dWheels+rWheel*3)*cnt
138 hasTensioner = int(len(circleList) == 6) #this is a special case with tensioners
139 # print('cnt=',cnt, ', tensioner=', hasTensioner)
140
141 for item in circleList:
142 item[0][0]+=offX
143 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
144 #create geometry:
145 reevingDict = CreateReevingCurve(circleList, drawingLinesPerCircle = 64,
146 removeLastLine=True, #allows closed curve
147 numberOfANCFnodes=nANCFnodes, graphicsNodeSize= 0.01)
148
149 del circleList[-1] #remove circles not needed for contact/visualization
150 del circleList[-1] #remove circles not needed for contact/visualization
151
152 gList=[]
153 if False: #visualize reeving curve, without simulation
154 gList = reevingDict['graphicsDataLines'] + reevingDict['graphicsDataCircles']
155
156 # print('belt length=',reevingDict['totalLength'])
157
158 oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
159 visualization=VObjectGround(graphicsData= gList)))
160
161 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
162 #create ANCF elements:
163 cableTemplate = Cable2D(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
164 physicsMassPerLength = rhoBeam*A,
165 physicsBendingStiffness = E*I,
166 physicsAxialStiffness = E*A,
167 physicsBendingDamping = dEI,
168 physicsAxialDamping = dEA,
169 physicsReferenceAxialStrain = preStretch, #prestretch
170 useReducedOrderIntegration=2,
171 visualization=VCable2D(drawHeight=2*h),
172 )
173
174 ancf = PointsAndSlopes2ANCFCable2D(mbs, reevingDict['ancfPointsSlopes'], reevingDict['elementLengths'],
175 cableTemplate, massProportionalLoad=gVec,
176 #fixedConstraintsNode0=[1,1,1,1], fixedConstraintsNode1=[1,1,1,1],
177 firstNodeIsLastNode=True, graphicsSizeConstraints=0.01)
178
179
180 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
181 #add contact:
182 if useContact:
183
184 if useGeneralContact:
185 gContact = mbs.AddGeneralContact()
186 gContact.verboseMode = 1
187 gContact.frictionProportionalZone = 0.1
188 gContact.ancfCableUseExactMethod = False
189 gContact.ancfCableNumberOfContactSegments = nContactSegments
190 ssx = 64#32 #search tree size
191 ssy = 12#32 #search tree size
192 ssz = 1 #search tree size
193 gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
194
195 sWheelRot = [] #sensors for angular velocity
196
197 nGround = mbs.AddNode(NodePointGround())
198 mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
199
200 if hasTensioner:
201 #add tensioning system
202 nTensioner = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[offX+0.4*dWheels,0,0],
203 visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
204 #gTensioner = [GraphicsDataOrthoCubePoint([0,0,dimZ],size=[0.1*rWheel,2.2*tensionerY,dimZ],color=color4orange)]
205 cyl0 = GraphicsDataCylinder([0,tensionerY+0.05*rWheel,dimZ],vAxis=[0,-2*tensionerY-0.1*rWheel,0],
206 radius=0.05*rWheel, color=color4orange)
207 cyl1 = GraphicsDataCylinder([0,tensionerY,dimZ],vAxis=[0.6*dWheels,-tensionerY,0], radius=0.05*rWheel, color=color4orange)
208 cyl2 = GraphicsDataCylinder([0,-tensionerY,dimZ],vAxis=[0.6*dWheels,tensionerY,0], radius=0.05*rWheel, color=color4orange)
209 gTensioner = [cyl0,cyl1,cyl2]
210
211 oTensioner = mbs.AddObject(ObjectRigidBody2D(physicsMass=wheelMass*10, physicsInertia=wheelInertia*10,
212 nodeNumber=nTensioner, visualization=
213 VObjectRigidBody2D(graphicsData=gTensioner)))
214 mTensionerCenter = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nTensioner))
215 # mbs.AddLoad(Force(markerNumber=mTensionerCenter, loadVector=[1,-g*wheelMass,0]))
216
217 mTensionerSupport = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oTensioner, localPosition=[0.6*dWheels,0,0]))
218 mTensionerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[offX+1*dWheels,0,0]))
219 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mTensionerGround, mTensionerSupport]))
220
221 mTensionerTop = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oTensioner, localPosition=[0,tensionerY,0]))
222 mTensionerBottom = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oTensioner, localPosition=[0,-tensionerY,0]))
223
224 if staticEqulibrium: #fix all wheels
225 mCoordinateTensioner = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nTensioner, coordinate=2))
226 cc = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateTensioner]))
227 fixWheelsList += [cc]
228
229 # print(mbs.GetMarker(mTensionerTop))
230 # print(mbs.GetMarker(mTensionerBottom))
231
232 for i, wheel in enumerate(circleList):
233 p = [wheel[0][0], wheel[0][1], 0] #position of wheel center
234 r = wheel[1]
235
236 rot0 = 0 #initial rotation
237 pRef = [p[0], p[1], rot0]
238 gList = [GraphicsDataCylinder(pAxis=[0,0,-0.5*dimZ],vAxis=[0,0,dimZ], radius=r-h,
239 color= color4dodgerblue, nTiles=64),
240 GraphicsDataArrow(pAxis=[0,0,dimZ], vAxis=[0.9*r,0,0], radius=0.01*r, color=color4orange)]
241
242 if i == 1+hasTensioner:
243 gList += [GraphicsDataOrthoCubePoint([0,-0.5*leverArm,-dimZ],size=[leverArm*0.1,leverArm,dimZ],color=color4red)]
244
245
246 omega0 = 0 #initial angular velocity
247 v0 = np.array([0,0,omega0])
248
249 nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
250 visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
251 oMass = mbs.AddObject(ObjectRigidBody2D(physicsMass=wheelMass, physicsInertia=wheelInertia,
252 nodeNumber=nMass, visualization=
253 VObjectRigidBody2D(graphicsData=gList)))
254 mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
255 mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p))
256
257 if hasTensioner==0 or i==0 or i==2:
258 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode]))
259 elif hasTensioner==1:
260 if i==1:
261 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mTensionerTop, mNode]))
262 elif i==3:
263 # mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode]))
264 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mTensionerBottom, mNode]))
265
266
267 sWheelRot += [mbs.AddSensor(SensorNode(nodeNumber=nMass,
268 fileName='solution/beltDrive'+str(cnt)+'wheel'+str(i)+'angVel.txt',
269 outputVariableType=exu.OutputVariableType.AngularVelocity))]
270
271 mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
272 if staticEqulibrium: #fix all wheels
273 cc = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel]))
274 fixWheelsList += [cc]
275
276 if i == 0:
277 if velocityControl:
278 fact = 1 #drive direction
279 if circleList[0][2] == 'R':
280 fact = -1
281 cc = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
282 velocityLevel=True,offset=fact,
283 offsetUserFunction=UFvelocityDrive))
284 velControlList += [cc]
285
286 elif i == 1+hasTensioner:
287 mLever = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[0,-leverArm,0]))
288 mbs.AddLoad(Force(markerNumber=mLever, loadVector=[0,-g*massArmTip,0]))
289
290 if useGeneralContact:
291 frictionMaterialIndex=0
292 # if hasTensioner==0 or i==0 or i==2:
293 gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=contactStiffness,
294 contactDamping=contactDamping, frictionMaterialIndex=frictionMaterialIndex)
295 else: #conventional contact, one contact per element and pulley
296 cableList = ancf[1]
297 mCircleBody = mNode
298 for k in range(len(cableList)):
299 nseg = nContactSegments
300 initialGapList = [0.1]*nseg + [-2]*nseg + [0]*nseg #initial gap of 0., isStick (0=slip, +-1=stick, -2 undefined initial state), lastStickingPosition (0)
301
302 mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[k],
303 numberOfSegments = nseg, verticalOffset=-h/2))
304 nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,
305 numberOfDataCoordinates=nseg*3 ))
306
307 co = mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mCircleBody, mCable], nodeNumber = nodeDataContactCable,
308 numberOfContactSegments=nseg,
309 contactStiffness = contactStiffness,
310 contactDamping = contactDamping,
311 frictionVelocityPenalty = frictionVelocityPenalty,
312 frictionStiffness = contactFrictionStiffness,
313 frictionCoefficient = dryFriction,
314 circleRadius = r,
315 # useSegmentNormals=False,
316 visualization=VObjectContactFrictionCircleCable2D(showContactCircle=False)))
317 contactObjects+=[co]
318
319
320
321 if useGeneralContact:
322 halfHeight = 0.5*h*0
323 for oIndex in ancf[1]:
324 gContact.AddANCFCable(objectIndex=oIndex, halfHeight=halfHeight, #halfHeight should be h/2, but then cylinders should be smaller
325 contactStiffness=contactStiffness, contactDamping=contactDamping,
326 frictionMaterialIndex=0)
327
328 frictionMatrix = np.zeros((2,2))
329 frictionMatrix[0,0]=int(useFriction)*dryFriction
330 frictionMatrix[0,1]=0 #no friction between some rolls and cable
331 frictionMatrix[1,0]=0 #no friction between some rolls and cable
332 gContact.SetFrictionPairings(frictionMatrix)
333
334
335mbs.Assemble()
336
337simulationSettings = exu.SimulationSettings() #takes currently set values or default values
338
339simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
340simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
341simulationSettings.solutionSettings.writeSolutionToFile = True
342simulationSettings.solutionSettings.solutionWritePeriod = 0.005
343simulationSettings.solutionSettings.sensorsWritePeriod = 0.001
344# simulationSettings.displayComputationTime = True
345simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
346# simulationSettings.displayStatistics = True
347
348simulationSettings.timeIntegration.endTime = tEnd
349simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
350simulationSettings.timeIntegration.stepInformation= 3+8+32+128+256
351
352simulationSettings.timeIntegration.verboseMode = 1
353
354simulationSettings.timeIntegration.newton.useModifiedNewton = True
355simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
356simulationSettings.timeIntegration.discontinuous.maxIterations = 5+2
357# if not useGeneralContact:
358# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-4
359
360simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
361
362SC.visualizationSettings.general.circleTiling = 24
363SC.visualizationSettings.loads.show=False
364SC.visualizationSettings.nodes.defaultSize = 0.005
365SC.visualizationSettings.nodes.show=False
366SC.visualizationSettings.openGL.multiSampling = 4
367
368SC.visualizationSettings.general.textSize=14
369SC.visualizationSettings.general.showSolverInformation=False
370SC.visualizationSettings.general.renderWindowString = 'Comparision of belt drive configurations:\n - prescribed drive velocity at left pulley\n - lever arm under gravity attached to right pulley\n - showing axial forces as vertical bars along beams'
371SC.visualizationSettings.window.renderWindowSize=[1920,1080]
372SC.visualizationSettings.general.drawCoordinateSystem=False
373
374if True:
375 SC.visualizationSettings.bodies.beams.axialTiling = 1
376 SC.visualizationSettings.bodies.beams.drawVertical = True
377 SC.visualizationSettings.bodies.beams.drawVerticalLines = True
378
379 SC.visualizationSettings.contour.outputVariableComponent=0
380 SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
381
382 SC.visualizationSettings.bodies.beams.drawVerticalFactor = 0.005
383 SC.visualizationSettings.bodies.beams.drawVerticalOffset = -6*0
384
385 SC.visualizationSettings.bodies.beams.reducedAxialInterploation = True
386
387
388#visualize contact:
389SC.visualizationSettings.connectors.showContact = True
390SC.visualizationSettings.contact.showContactForces = True
391SC.visualizationSettings.contact.contactForcesFactor = 0.025
392
393if False:
394 SC.visualizationSettings.contact.showSearchTree =True
395 SC.visualizationSettings.contact.showSearchTreeCells =True
396 SC.visualizationSettings.contact.showBoundingBoxes = True
397
398if useGraphics:
399 exu.StartRenderer()
400 mbs.WaitForUserToContinue()
401
402#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
403simulationSettings.staticSolver.adaptiveStep = False
404simulationSettings.staticSolver.loadStepGeometric = True;
405simulationSettings.staticSolver.loadStepGeometricRange=1e3
406simulationSettings.staticSolver.numberOfLoadSteps = 10
407#simulationSettings.staticSolver.useLoadFactor = False
408simulationSettings.staticSolver.stabilizerODE2term = 1e5*10
409simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
410simulationSettings.staticSolver.newton.absoluteTolerance = 1e-6
411simulationSettings.staticSolver.newton.numericalDifferentiation.minimumCoordinateSize = 1 #needed for static solver to converge
412simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-5 #needed for static solver to converge
413if staticEqulibrium: #precompute static equilibrium
414 for velControl in velControlList:
415 mbs.SetObjectParameter(velControl, 'activeConnector', False)
416
417 for obj in contactObjects:
418 mbs.SetObjectParameter(obj, 'frictionCoefficient', 0.)
419 mbs.SetObjectParameter(obj, 'frictionStiffness', 1) #do not set to zero, as it needs to do some initialization...
420
421 mbs.SolveStatic(simulationSettings, updateInitialValues=True)
422
423
424 for obj in contactObjects:
425 mbs.SetObjectParameter(obj, 'frictionCoefficient', dryFriction)
426 mbs.SetObjectParameter(obj, 'frictionStiffness', contactFrictionStiffness)
427
428 for velControl in velControlList:
429 mbs.SetObjectParameter(velControl, 'activeConnector', True)
430
431 for obj in fixWheelsList:
432 mbs.SetObjectParameter(obj, 'activeConnector', False)
433
434
435mbs.SolveDynamic(simulationSettings,
436 # solverType=exu.DynamicSolverType.TrapezoidalIndex2
437 ) #183 Newton iterations, 0.114 seconds
438
439
440
441if useGraphics and True:
442 SC.visualizationSettings.general.autoFitScene = False
443 SC.visualizationSettings.general.graphicsUpdateInterval=0.02
444
445 sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
446 mbs.SolutionViewer(sol)
447
448
449if useGraphics:
450 SC.WaitForRenderEngineStopFlag()
451 exu.StopRenderer() #safely close rendering window!
452
453 # if True:
454 #
455 # mbs.PlotSensor(sensorNumbers=[sAngVel[0],sAngVel[1]], components=2, closeAll=True)
456 # mbs.PlotSensor(sensorNumbers=sMeasureRoll, components=1)