chainDriveExample.py
You can view and download this file on Github: chainDriveExample.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test of a chain drive for ContactCurveCircles
5#
6# Author: Johannes Gerstmayr
7# Date: 2024-11-01
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 exudyn as exu
14#a local run must not open a plot window when windows are suppressed (#2477)
15if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
16from exudyn.utilities import *
17import exudyn.graphics as graphics
18import numpy as np
19from numpy import sin, cos, arcsin #used below; not exported by exudyn.utilities
20
21#**function: unique drawing for chain plate; reference position is at first roller;
22#link is oriented into x-axis, thickness in z-axis, height in y-axis (height of rectangular part);
23#zOff is the center of the thickness of the plate
24#returns a single graphicsData object
25def GraphicsChainPlate(lengthLink, heightLink, thicknessLink, radiusLink,
26 pOff=[0,0,0], color=[-1,-1,-1,-1], nTiles=16,
27 addEdges = False, addFaces=True):
28 if heightLink > 2*radiusLink:
29 raise ValueError('GraphicsChainPlate: height must by <= 2*radius!')
30 from numpy import sin, cos, arcsin
31
32 gList = []
33 vertices = []
34 segments = []
35 lastIndex = 2*(nTiles+1)-1
36 phiStart = arcsin(0.5*heightLink/radiusLink)
37
38 #check if rectangular section can be drawn:
39 if lengthLink - 2*cos(phiStart)*radiusLink < 0:
40 phiStart = pi/2
41
42 phiEnd = 2*pi-phiStart
43 phiRange = phiEnd-phiStart
44
45 for i in range(nTiles+1):
46 phi = i/nTiles*phiRange+phiStart
47 vertices.append([radiusLink*cos(phi),radiusLink*sin(phi)])
48
49 for i in range(nTiles+1):
50 phi = i/nTiles*phiRange+pi+phiStart
51 vertices.append([radiusLink*cos(phi)+lengthLink,radiusLink*sin(phi)])
52
53 segments = np.vstack((np.arange(2*(nTiles+1)),
54 np.arange(2*(nTiles+1))+1)).T
55 segments[-1,-1] = 0
56
57 return graphics.SolidExtrusion(vertices, segments.tolist(), thicknessLink, pOff=np.array(pOff)-[0,0,0.5*thicknessLink],
58 color=color,addEdges=addEdges,addFaces=addFaces)
59
60
61#**function: unique drawing for chain plate; reference position is at first roller;
62#link is oriented into x-axis, thickness in z-axis, height in y-axis (height of rectangular part);
63#zOff is the center of the thickness of the plate
64#returns a single graphicsData object
65def GraphicsChainLink(lengthLink, heightLink, thicknessLink, radiusLink, innerWidthLink,
66 radiusRoller, radiusPin, withPin, pOff=[0,0,0],
67 colorLink=graphics.color.grey, colorRoller=[0.7,0.7,0.8,1], colorPin=graphics.color.darkgrey,
68 nTiles=16, addEdges = False, addFaces=True):
69
70 gcl = []
71 gcl += [graphics.Cylinder([0,0,-0.5*(withPin)],[0,0,withPin], radius=radiusPin,
72 color=colorRoller, nTiles=nTiles, addEdges=addEdges, addFaces=addFaces)]
73 gcl += [GraphicsChainPlate(lengthLink, heightLink, thicknessLink, radiusLink,
74 pOff=[0,0,0.5*(innerWidthLink+thicknessLink)], color=colorLink,
75 nTiles=nTiles, addEdges=addEdges, addFaces=addFaces)]
76 gcl += [GraphicsChainPlate(lengthLink, heightLink, thicknessLink, radiusLink,
77 pOff=[0,0,-0.5*(innerWidthLink+thicknessLink)], color=colorLink,
78 nTiles=nTiles, addEdges=addEdges, addFaces=addFaces)]
79 gcl += [graphics.Cylinder([0,0,-0.5*(innerWidthLink)],[0,0,innerWidthLink], radius=radiusRoller,
80 color=colorRoller, nTiles=nTiles, addEdges=addEdges, addFaces=addFaces)]
81 return gcl
82
83#unique drawing for roller chain link:
84
85
86useGraphics = True #without test
87SC = exu.SystemContainer()
88mbs = SC.AddSystem()
89
90pi = np.pi
91mass = 1
92nTeeth = 16
93nLinks3 = 7 #free links between sprockets
94phiTooth = 2*pi/nTeeth
95rRoller = 0.005
96lLink = 0.03 #Length of link; l = 2*r*sin(phi/2)
97wLink = 0.012 #inner with of link
98tpLink = 0.002 #thickness of plate
99hpLink = 0.012 #height of plate
100rpLink = 1.5*rRoller #radius of plate
101rSprocketPitch = lLink/(2*sin(phiTooth/2)) #radius to center of holes in gear, given by chain link length and nTeeth!
102rSprocketInner = rSprocketPitch - rRoller
103wSprocket = wLink*0.95
104addSecondWheel = True
105contactStiffness = 1e4*5
106contactDamping = 1e1
107
108exu.Print('Chain geometry:')
109exu.Print(' lLink =',lLink)
110exu.Print(' rRoller=',rRoller)
111exu.Print(' rPitch =',rSprocketPitch)
112exu.Print(' rInner =',rSprocketPitch)
113
114tEnd = 10 #end time of simulation
115stepSize = 1e-4
116g = 9.81
117
118gGround = graphics.CheckerBoard(point = [0,-rSprocketInner*1.5,-wSprocket*2], size = 0.5)
119nGround = mbs.AddNode(NodePointGround())
120mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
121
122oGround = mbs.CreateGround(referencePosition=[0,0,0],graphicsDataList=[gGround])
123mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
124
125#ground node
126gListSprocket = [graphics.Cylinder([0,0,-0.5*wSprocket],[0,0,wSprocket],
127 radius=rSprocketPitch, color=graphics.color.grey, nTiles=32)]
128
129pList = []
130nSprocket = nTeeth*4 #8
131for i in range(nSprocket):
132 angle = (i-0.5)/(nSprocket)*2*pi
133 r = rSprocketInner + 2*rRoller*(i%4 > 1)
134 x =-r*cos(angle)
135 y =-r*sin(angle)
136 pList.append([x,y])
137
138segments = np.vstack((np.arange(nSprocket),
139 np.arange(nSprocket)+1)).T
140segments[-1,-1] = 0
141segmentsData = np.zeros((nSprocket,4))
142segmentsData[:,0:2] = pList
143segmentsData[:,2:4] = np.roll(pList,2) #roll is element wise on rows and columns, therefore 2>shift one row
144
145print(segmentsData[:5,:])
146
147gListSprocket = [
148 graphics.SolidExtrusion(pList, segments.tolist(), wSprocket, pOff=[0,0,-0.5*wSprocket],
149 color=graphics.color.lightgrey,addEdges=True,addFaces=True),
150 graphics.Cylinder(pAxis=[0,0,-wSprocket], vAxis=[0,0,2*wSprocket], radius = rRoller, color=graphics.color.grey),
151 graphics.Brick(size=[0.5*wSprocket,wSprocket*2,wSprocket*1.05], color=graphics.color.red)
152]
153
154
155
156#sprocket wheel:
157sprocket = mbs.CreateRigidBody(referencePosition=[0,0,0],
158 inertia=InertiaCylinder(density=100,length=wSprocket,outerRadius=rSprocketPitch,axis=2),
159 gravity = [0,-g,0],
160 graphicsDataList=gListSprocket,
161 create2D=True
162 )
163mSprocket = mbs.AddMarker(MarkerBodyRigid(bodyNumber=sprocket))
164mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGround, mSprocket],
165 visualization=VRevoluteJoint2D(drawSize=rRoller)) )
166
167#add prescribed velocity:
168if True:
169 def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
170 vStart = 0
171 vEnd = pi
172 return SmoothStep(t, 0, 0.5, vStart, vEnd)
173
174 nodeSprocket = mbs.GetObject(sprocket)['nodeNumber']
175 mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodeSprocket, coordinate=2))
176 velControl = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
177 velocityLevel=True, offsetUserFunction_t= UFvelocityDrive,
178 visualization=VCoordinateConstraint(show = False)))#UFvelocityDrive
179else:
180 mbs.AddLoad(Torque(markerNumber=mSprocket, loadVector=[0,0,-0.1]))
181
182sVel1 = mbs.AddSensor(SensorBody(bodyNumber=sprocket, storeInternal=True, outputVariableType=exu.OutputVariableType.AngularVelocity))
183
184if addSecondWheel:
185 oGround2 = mbs.CreateGround(referencePosition=[0,-nLinks3*lLink,0],graphicsDataList=[])
186 mGround2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround2))
187 sprocket2 = mbs.CreateRigidBody(referencePosition=[0,-nLinks3*lLink,0],
188 inertia=InertiaCylinder(density=100,length=wSprocket,outerRadius=rSprocketPitch,axis=2),
189 gravity = [0,-g,0],
190 graphicsDataList=gListSprocket,
191 create2D=True
192 )
193 # nn = mbs.GetObject(sprocket)['nodeNumber']
194 mSprocket2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=sprocket2))
195 mbs.AddObject(SphericalJoint(markerNumbers=[mGround2, mSprocket2],
196 constrainedAxes=[1,0,0],
197 visualization=VSphericalJoint(jointRadius=0.25*rRoller)) )
198
199 mbs.AddObject(LinearSpringDamper(markerNumbers=[mGround2, mSprocket2],
200 stiffness = 1000, damping=100,
201 axisMarker0=[0,1,0], force = 10,
202 visualization=VLinearSpringDamper(drawSize=0.01),
203 ))
204
205 dTSD = 0.01
206 mbs.AddObject(TorsionalSpringDamper(markerNumbers = [mGround, mSprocket2], damping = dTSD))
207 #mbs.AddLoad(Torque(markerNumber=mSprocket2, loadVector=[0,0,0.2]))
208 sVel2 = mbs.AddSensor(SensorBody(bodyNumber=sprocket2, storeInternal=True, outputVariableType=exu.OutputVariableType.AngularVelocity))
209
210
211#graphics for chain link:
212gChainLink1 = GraphicsChainLink(lLink, hpLink, tpLink, rpLink, wLink,
213 rRoller, 0.3*rRoller, wLink+4.5*tpLink, pOff=[0,0,0],
214 nTiles=32, addEdges = True)
215gChainLink2 = GraphicsChainLink(lLink, hpLink, tpLink, rpLink, wLink+tpLink*2,
216 rRoller, 0.3*rRoller, wLink+4.5*tpLink, pOff=[0,0,0],
217 nTiles=32, addEdges = True)
218
219#create geometry of hanging chain:
220posList = []
221rotList = []
222for i in range(nLinks3):
223 x=-rSprocketPitch
224 y=-(nLinks3-i)*lLink
225 posList.append([x,y])
226 rotList.append(pi/2)
227
228for i in range(int(nTeeth/2)):
229 angle = i/(nTeeth/2)*pi
230 x=-rSprocketPitch*cos(angle)
231 y= rSprocketPitch*sin(angle)
232 posList.append([x,y])
233 rotList.append(pi/2-angle-0.5/(nTeeth/2)*pi)
234
235for i in range(nLinks3):
236 x= rSprocketPitch
237 y=-(i+0)*lLink
238 posList.append([x,y])
239 rotList.append(-pi/2)
240
241for i in range(int(nTeeth/2)):
242 angle = -(i+0)/(nTeeth/2)*pi
243 x= rSprocketPitch*cos(angle)
244 y= rSprocketPitch*sin(angle)
245 posList.append([x,y-nLinks3*lLink])
246 rotList.append(-pi/2+angle-0.5/(nTeeth/2)*pi)
247
248oLinksList = []
249mLinksList = []
250
251prevMarker = None
252firstMarker = None
253for i in range(len(posList)):
254 pos = posList[i]
255 #note: in 2D case, the reference point = COM, which is not correct for chain; for correct representation, move ref. point to COM of link!
256 link = mbs.CreateRigidBody(referencePosition=pos+[0],
257 inertia=InertiaCuboid(1000,sideLengths=[lLink,wLink,wLink]),
258 referenceRotationMatrix=RotationMatrixZ(rotList[i]),
259 gravity = [0,-g,0],
260 graphicsDataList=gChainLink1 if i%2==0 else gChainLink2,
261 create2D=True
262 )
263 # nn = mbs.GetObject(sprocket)['nodeNumber']
264 oLinksList += [link]
265 mLinksList += [mbs.AddMarker(MarkerBodyRigid(bodyNumber=link))]
266
267 if prevMarker is not None:
268 mbs.AddObject(RevoluteJoint2D(markerNumbers=[prevMarker, mLinksList[-1]],
269 visualization=VRevoluteJoint2D(drawSize=rRoller)) )
270 else:
271 firstMarker = mLinksList[-1]
272
273 prevMarker = mbs.AddMarker(MarkerBodyRigid(bodyNumber=link, localPosition=[lLink,0,0]))
274
275#close chain:
276mbs.AddObject(RevoluteJoint2D(markerNumbers=[prevMarker, firstMarker],
277 visualization=VRevoluteJoint2D(drawSize=rRoller)) )
278
279
280if False: #put one link to ground:
281 posLast = mbs.GetMarkerOutput(prevMarker, variableType=exu.OutputVariableType.Position,
282 configuration=exu.ConfigurationType.Reference)
283
284 oGround2 = mbs.CreateGround(referencePosition=posLast)
285 mGround2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround2))
286 mbs.AddObject(RevoluteJoint2D(markerNumbers=[prevMarker, mGround2],
287 visualization=VRevoluteJoint2D(drawSize=rRoller)) )
288
289nGenericData = mbs.AddNode(NodeGenericData(initialCoordinates=[-1,0,0]*nSprocket,
290 numberOfDataCoordinates=3*nSprocket))
291mbs.AddObject(ObjectContactCurveCircles(markerNumbers=[mSprocket]+mLinksList, nodeNumber=nGenericData,
292 circlesRadii=[rRoller]*len(mLinksList), segmentsData=exu.MatrixContainer(segmentsData),
293 contactStiffness=contactStiffness, contactDamping=contactDamping))
294
295if addSecondWheel:
296 nGenericData2 = mbs.AddNode(NodeGenericData(initialCoordinates=[-1,0,0]*nSprocket,
297 numberOfDataCoordinates=3*nSprocket))
298 mbs.AddObject(ObjectContactCurveCircles(markerNumbers=[mSprocket2]+mLinksList, nodeNumber=nGenericData2,
299 circlesRadii=[rRoller]*len(mLinksList), segmentsData=exu.MatrixContainer(segmentsData),
300 contactStiffness=contactStiffness, contactDamping=contactDamping))
301
302
303mbs.Assemble()
304
305simulationSettings = exu.SimulationSettings()
306simulationSettings.solution.file.write = True
307simulationSettings.solution.file.writePeriod = 0.01
308simulationSettings.solution.sensors.writePeriod = stepSize #output interval
309simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
310simulationSettings.timeIntegration.endTime = tEnd
311#simulationSettings.timeIntegration.realtime.active = True
312# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
313# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
314
315simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
316simulationSettings.timeIntegration.newton.useModifiedNewton = True
317#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
318
319simulationSettings.show.statistics = True
320simulationSettings.timeIntegration.verboseMode = 1
321
322SC.visualizationSettings.view0.window.renderWindowSize=[1600,2000]
323SC.visualizationSettings.openGL.multiSampling=4
324#SC.visualizationSettings.view0.scene.facesTransparent=True
325SC.visualizationSettings.openGL.light0.shadow=0.3
326SC.visualizationSettings.loads.show = False
327SC.visualizationSettings.connectors.contactPointsDefaultSize = 0.001
328SC.visualizationSettings.connectors.showContact = True
329
330SC.renderer.Start() #start graphics visualization
331SC.renderer.DoIdleTasks() #wait for pressing SPACE bar to continue
332
333mbs.SolveDynamic(simulationSettings)
334
335SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
336SC.renderer.Stop() #safely close rendering window!
337
338#plot results:
339if False:
340 mbs.PlotSensor([sVel1,sVel2], components=[2,2], closeAll=True)
341 import matplotlib.pyplot as plt
342 plt.show(block=True) #for figures to stay open at end of plot routines
343
344#+++++++++++++++++++++++++++++++++++++++++++++++++++++
345
346mbs.SolutionViewer()