ANCFbeltDrive.py
You can view and download this file on Github: ANCFbeltDrive.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: ANCF cable elements in contact with circles defined by GeneralContact;
5# this is an advanced example for a belt drive according to literature;
6# note that the belt is very thin in order to avoid bending effects, to make this comparible to analytical results
7#
8# Author: Johannes Gerstmayr
9# Date: 2022-01-31
10#
11# 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.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14import exudyn as exu
15from exudyn.itemInterface import *
16import numpy as np
17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
18import exudyn.graphics as graphics
19from exudyn.beams import *
20from math import sin, cos, sqrt, pi
21
22testIsActive = exu.sys.get('testIsActive', False)
23
24SC = exu.SystemContainer()
25mbs = SC.AddSystem()
26
27#background
28rect = [-1,-1.5,3,1.5] #xmin,ymin,xmax,ymax
29background0 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
30oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
31 visualization=VObjectGround(graphicsData= [background0])))
32nGround = mbs.AddNode(NodePointGround())
33mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
34
35
36#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37#contact
38doImplicit = True
39useContact = True
40useFriction = True
41dryFriction = 0.5
42contactStiffness = 1e5
43contactDamping = 1e-3*contactStiffness
44
45if useContact:
46 gContact = mbs.AddGeneralContact()
47 gContact.verboseMode = 1*0
48 gContact.frictionProportionalZone = 1
49 gContact.ancfCableUseExactMethod = False
50 gContact.ancfCableNumberOfContactSegments = 8
51 ssx = 16#32 #search tree size
52 ssy = 8#32 #search tree size
53 ssz = 1 #search tree size
54 gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
55 #gContact.SetSearchTreeBox(pMin=np.array([-1,-1,-1]), pMax=np.array([4,1,1]))
56torque=10*2
57#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
58#wheel:
59dampingWheel1 = 1 #5 #add breaking torque, to limit velocity
60#cable:
61
62numberOfElements = 16 # per section; 32 recommended
63curvedRefConf=False # this flag could initialize the elements to be produced curved -> not suitable for belt drive!
64L=2 # length of ANCF element in m
65E=1e10 # Young's modulus of ANCF element in N/m^2
66rhoBeam=1000 # density of ANCF element in kg/m^3
67b=0.002 # width of rectangular ANCF element in m
68h=0.002 # height of rectangular ANCF element in m
69A=b*h # cross sectional area of ANCF element in m^2
70I=b*h**3/12 # second moment of area of ANCF element in m^4
71dEI = 0*1e-3*E*I
72dEA = 1e-2*E*A
73# f=3*E*I/L**2 # tip load applied to ANCF element in N
74g=-9.81
75dimZ = b #z.dimension
76preStretch=-0.002
77# exu.Print("load f="+str(f))
78# exu.Print("EI="+str(E*I))
79
80# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
81# mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
82
83cableTemplate = Cable2D(#length = L / nElements, #set in GenerateBeamElementsAlongLine(...)
84 massPerLength = rhoBeam*A,
85 bendingStiffness = E*I,
86 axialStiffness = E*A,
87 bendingDamping = dEI,
88 axialDamping = dEA,
89 referenceAxialStrain = preStretch, #prestretch
90 #nodeNumbers = [0, 0], #will be filled in GenerateBeamElementsAlongLine(...)
91 visualization=VCable2D(drawHeight=2*h),
92 )
93exu.Print("pre-stretch force=", preStretch*E*A)
94exu.Print("beam mass per length=", rhoBeam*A)
95#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96#create belt drive:
97distanceWheels = 2 #distance of wheel centers
98wheelCenter0 = np.array([0,0,0])
99wheelCenter1 = np.array([distanceWheels,0,0])
100
101rWheel0 = 0.5
102rWheel1 = 0.5
103rRoll = 0.1
104rollCenter1 = np.array([0.5*distanceWheels,-rWheel0+rRoll,0])
105
106mWheel = 4
107mRoll = 1e-6 #small as compared to beam with rhoA=4e-3
108kRoll = 1e-2
109dRelRoll = 0.005
110fRoll = 0.01
111
112velocityControl = True
113
114yAxis = np.array([0,1.,0])
115ancfList=[]
116
117if True:
118 startAngle = -pi
119 arcAngle = -pi
120 positionOfNode0 = wheelCenter0-rWheel0*yAxis # starting point of line
121 ancf=GenerateCircularArcANCFCable2D(mbs, positionOfNode0,
122 rWheel0, startAngle, arcAngle, numberOfElements,
123 cableTemplate,
124 massProportionalLoad = [0,g,0], #optionally add gravity
125 #fixedConstraintsNode0 = [1,1,1,1], #add constraints for pos and rot (r'_y)
126 #fixedConstraintsNode1 = [1,1,1,1],
127 setCurvedReferenceConfiguration=curvedRefConf,
128 #verboseMode=True
129 )
130 ancfList+=[ancf]
131 beamInfo=GenerateBeamElementsAlongLine(mbs,
132 ancf[3][-1], wheelCenter1+rWheel1*yAxis,
133 numberOfElements,
134 beamTemplate=cableTemplate, #this defines the beam element properties
135 gravity = [0,g,0], #optionally add gravity
136 # groundConstraintsStart = [1,1,1,1], #add constraints for pos and rot (r'_y)
137 # groundConstraintsEnd = [1,1,1,1],
138 nodeNumberStart=ancf[0][-1]
139 )
140 ancf = [beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']] #as the list of GenerateCircularArcANCFCable2D
141 ancfList+=[ancf]
142
143 if True:
144 startAngle = 0
145 arcAngle = -pi
146 ancf=GenerateCircularArcANCFCable2D(mbs, ancf[3][-1],
147 rWheel1, startAngle, arcAngle, numberOfElements,
148 cableTemplate,
149 massProportionalLoad = [0,g,0], #optionally add gravity
150 #fixedConstraintsNode0 = [1,1,1,1], #add constraints for pos and rot (r'_y)
151 #fixedConstraintsNode1 = [1,1,1,1],
152 setCurvedReferenceConfiguration=curvedRefConf,
153 nodeNumber0=ancf[0][-1]
154 )
155 ancfList+=[ancf]
156 beamInfo=GenerateBeamElementsAlongLine(mbs,
157 ancf[3][-1], ancfList[0][3][0],
158 numberOfElements,
159 beamTemplate=cableTemplate, #this defines the beam element properties
160 gravity = [0,g,0], #optionally add gravity
161 nodeNumberStart=ancf[0][-1],
162 nodeNumberEnd=ancfList[0][0][0]
163 )
164 ancf = [beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']] #as the list of GenerateCircularArcANCFCable2D
165 ancfList+=[ancf]
166
167
168
169def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
170 vMax = 10 #5m/s
171 if 2*t < vMax:
172 return 2*t
173 return vMax
174
175#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
176sAngVel=[]
177#add contact:
178if useContact:
179
180
181 halfHeight = 0.5*h*0
182 wheels = [{'center':wheelCenter0, 'radius':rWheel0-halfHeight, 'mass':mWheel},
183 {'center':wheelCenter1, 'radius':rWheel1-halfHeight, 'mass':mWheel},
184 {'center':rollCenter1, 'radius':rRoll-halfHeight, 'mass':mRoll}, #small mass for roll, not to influence beam
185 ]
186
187 for i, wheel in enumerate(wheels):
188 r = wheel['radius']
189 p = wheel['center']
190 mass = wheel['mass']
191 rot0 = 0 #initial rotation
192 pRef = [p[0], p[1], rot0]
193 gList = [graphics.Cylinder(vAxis=[0,0,dimZ], radius=r*0.99, #draw smaller to see cable element
194 color= graphics.color.dodgerblue, nTiles=32),
195 graphics.Arrow(pAxis=[0,0,0.02*r], vAxis=[r,0,0], radius=0.02*r, color=graphics.color.orange)]
196
197 omega0 = 0 #initial angular velocity
198 v0 = np.array([0,0,omega0])
199
200 RBinertia = InertiaCylinder(mass/(r**2*np.pi*b), b, r, axis=2)
201
202 nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
203 visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
204 oMass = mbs.AddObject(ObjectRigidBody2D(mass=RBinertia.mass, inertia=RBinertia.GetInertia6D()[2],
205 nodeNumber=nMass, visualization=
206 VObjectRigidBody2D(graphicsData=gList)))
207 mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
208 mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p))
209 frictionMaterialIndex=0
210 if i < 2:
211 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode]))
212 else:
213 #add measurement roll:
214 dRoll = dRelRoll * (2*sqrt(kRoll*mRoll))
215 exu.Print('measurement roll resonance=',sqrt(kRoll/mRoll)/(2*pi))
216 frictionMaterialIndex=1 #no friction
217 mCWheel0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0))
218 mCWheel1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1))
219 mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCWheel0]))
220 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCWheel1],
221 stiffness = 1e-2, damping=dRoll,
222 offset=-fRoll/kRoll,
223 visualization=VCoordinateSpringDamper(show=False)))
224 sMeasureRoll = mbs.AddSensor(SensorNode(nodeNumber=nMass,
225 storeInternal=True, #read by PlotSensor, not from a file (#2492)
226 outputVariableType=exu.OutputVariableType.Displacement))
227
228 if i == 0:
229 if velocityControl:
230 mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
231 mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
232 velocityLevel=True, offsetUserFunction_t=UFvelocityDrive))
233
234 else:
235 mbs.AddLoad(LoadTorqueVector(markerNumber=mNode, loadVector=[0,0,torque]))
236 if i == 1:
237 mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
238 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel],
239 damping=dampingWheel1,
240 visualization=VCoordinateSpringDamper(show=False)))
241
242 gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=contactStiffness,
243 contactDamping=contactDamping, frictionMaterialIndex=frictionMaterialIndex)
244
245 sAngVel += [mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True, #(#2492)
246 outputVariableType=exu.OutputVariableType.AngularVelocity))]
247
248 allCables = []
249 for ancf in ancfList:
250 allCables += ancf[1]
251
252 for oIndex in allCables:
253 gContact.AddANCFCable(objectIndex=oIndex, halfHeight=halfHeight, #halfHeight should be h/2, but then cylinders should be smaller
254 contactStiffness=contactStiffness, contactDamping=contactDamping,
255 frictionMaterialIndex=0)
256
257 frictionMatrix = np.zeros((2,2))
258 frictionMatrix[0,0]=int(useFriction)*dryFriction
259 frictionMatrix[0,1]=0 #no friction between some rolls and cable
260 frictionMatrix[1,0]=0 #no friction between some rolls and cable
261 gContact.SetFrictionPairings(frictionMatrix)
262
263
264
265
266mbs.Assemble()
267#exu.Print(mbs)
268
269simulationSettings = exu.SimulationSettings() #takes currently set values or default values
270
271simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
272simulationSettings.solution.file.name = 'solution/coordinatesSolution.txt'
273simulationSettings.solution.file.write = (not testIsActive) #only the SolutionViewer reads it (#2492)
274simulationSettings.solution.file.writePeriod = 0.005
275simulationSettings.solution.sensors.writePeriod = 0.001
276#simulationSettings.show.computationTime = True
277#1 thread, because the multithreaded assembly sums in a non-deterministic order: with 4
278#threads the result of this model moved by up to 5.4e-14 between runs, past the 5e-14 test
279#tolerance, and at 16 elements per section 1 thread is also FASTER (0.36s vs 0.69s). Use more
280#threads for > 100 ANCF elements, but then the value is not a reference any more (#2368)
281simulationSettings.parallel.numberOfThreads = 1
282simulationSettings.show.statistics = True
283
284doDynamic = True
285tEnd = 0.1
286h = 0.5e-3
287simulationSettings.timeIntegration.endTime = tEnd
288simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
289simulationSettings.timeIntegration.stepInformation= 3+128+256
290
291simulationSettings.timeIntegration.verboseMode = 1 #otherwise, load steps are shown ...
292
293simulationSettings.timeIntegration.newton.useModifiedNewton = True
294
295#+++++++++++++++++++++++++++++++++++++++++++++
296#static:
297simulationSettings.staticSolver.verboseMode = 1 #otherwise, load steps are shown ...
298simulationSettings.staticSolver.numberOfLoadSteps = 40
299simulationSettings.staticSolver.discontinuous.iterationTolerance = 1
300simulationSettings.staticSolver.stabilizerODE2term = 2 #may only act on position degrees of freedom
301#+++++++++++++++++++++++++++++++++++++++++++++
302
303SC.visualizationSettings.view0.scene.drawWorldBasis=True
304SC.visualizationSettings.nodes.show = True
305SC.visualizationSettings.nodes.defaultSize = h*20
306SC.visualizationSettings.loads.show = False
307
308SC.visualizationSettings.contour.outputVariableComponent=0
309SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
310
311#visualize contact:
312if False:
313 SC.visualizationSettings.contact.showSearchTree =True
314 SC.visualizationSettings.contact.showSearchTreeCells =True
315 SC.visualizationSettings.contact.showBoundingBoxes = True
316
317if not testIsActive:
318 SC.renderer.Start()
319 SC.renderer.DoIdleTasks()
320
321if doDynamic :
322 mbs.SolveDynamic(simulationSettings) #183 Newton iterations, 0.114 seconds
323else:
324 mbs.SolveStatic(simulationSettings) #183 Newton iterations, 0.114 seconds
325
326
327if not testIsActive:
328 SC.visualizationSettings.general.autoFitScene = False
329 SC.visualizationSettings.general.graphicsUpdateInterval=0.02
330
331 sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
332 exu.Print('start SolutionViewer')
333 mbs.SolutionViewer(sol)
334
335
336if not testIsActive:
337 SC.renderer.DoIdleTasks()
338 SC.renderer.Stop() #safely close rendering window!
339
340 mbs.PlotSensor(sensorNumbers=[sAngVel[0],sAngVel[1]], components=2, closeAll=True)
341 mbs.PlotSensor(sensorNumbers=sMeasureRoll, components=1)
342
343
344#print representative result:
345#sol = mbs.systemData.GetODE2Coordinates()
346#n = len(sol)
347measurePos = mbs.GetSensorValues(sMeasureRoll)
348
349exu.Print('solution of ANCFbeltDrive (sensor y-disp): ',measurePos[1])
350#updated 2026-09-17; the old value was not reproducible - missing simulation parameters.
351#Re-measured single-threaded (see numberOfThreads above), which IS reproducible bit for bit;
352#the 4-thread value was -0.0011715990134786858 (#2368)
353exu.sys['testResult'] = measurePos[1] # use y-coordinate