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