reevingSystem.py
You can view and download this file on Github: reevingSystem.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example to calculate rope line of reeving system
5#
6# Author: Johannes Gerstmayr
7# Date: 2022-02-02
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
14from exudyn.itemInterface import *
15from exudyn.utilities import *
16from exudyn.beams import *
17
18import numpy as np
19from math import sin, cos, pi, sqrt , asin, acos, atan2
20import copy
21
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24
25
26
27#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
28#settings:
29useGraphics= True
30useContact = True
31tEnd = 2 #end time of dynamic simulation
32h = 1e-3 #step size
33useFriction = True
34dryFriction = 0.5
35contactStiffness = 2e5
36contactDamping = 1e-3*contactStiffness
37
38wheelMass = 1
39wheelInertia = 0.01
40
41rotationDampingWheels = 0.00 #proportional to rotation speed
42torque = 1
43
44#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
45#create circles
46#complicated shape:
47nANCFnodes = 200
48h = 0.25e-3
49preStretch=-0.001
50circleList = [[[0,0],0.3,'L'],
51 [[1,0],0.3,'R'],
52 [[1,1],0.3,'R'],
53 [[2,1],0.4,'R'],
54 [[3,1],0.4,'R'],
55 [[4,1],0.4,'L'],
56 [[5,1],0.4,'R'],
57 [[5,-1],0.3,'L'],
58 [[0,-1],0.4,'L'],
59 [[0,0],0.3,'L'],
60 [[1,0],0.3,'R'],
61 ]
62
63#simple shape:
64# nANCFnodes = 50
65# preStretch=-0.005
66# circleList = [[[0,0],0.2,'L'],
67# [[0,1],0.2,'L'],
68# [[0.8,0.8],0.4,'L'],
69# [[1,0],0.2,'L'],
70# [[0,0],0.2,'L'],
71# [[0,1],0.2,'L'],
72# ]
73
74#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
75#create geometry:
76reevingDict = CreateReevingCurve(circleList, drawingLinesPerCircle = 64,
77 removeLastLine=True, #allows closed curve
78 numberOfANCFnodes=nANCFnodes, graphicsNodeSize= 0.01)
79
80del circleList[-1] #remove circles not needed for contact/visualization
81del circleList[-1] #remove circles not needed for contact/visualization
82
83gList=[]
84if False: #visualize reeving curve, without simulation
85 gList = reevingDict['graphicsDataLines'] + reevingDict['graphicsDataCircles']
86
87oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
88 visualization=VObjectGround(graphicsData= gList)))
89
90#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
91#create ANCF elements:
92
93gVec = [0,-9.81,0] # gravity
94E=1e9 # Young's modulus of ANCF element in N/m^2
95rhoBeam=1000 # density of ANCF element in kg/m^3
96b=0.002 # width of rectangular ANCF element in m
97h=0.002 # height of rectangular ANCF element in m
98A=b*h # cross sectional area of ANCF element in m^2
99I=b*h**3/12 # second moment of area of ANCF element in m^4
100dEI = 1e-3*E*I #bending proportional damping
101dEA = 1e-2*E*A #axial strain proportional damping
102
103dimZ = b #z.dimension
104
105cableTemplate = Cable2D(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
106 physicsMassPerLength = rhoBeam*A,
107 physicsBendingStiffness = E*I,
108 physicsAxialStiffness = E*A,
109 physicsBendingDamping = dEI,
110 physicsAxialDamping = dEA,
111 physicsReferenceAxialStrain = preStretch, #prestretch
112 visualization=VCable2D(drawHeight=h),
113 )
114
115ancf = PointsAndSlopes2ANCFCable2D(mbs, reevingDict['ancfPointsSlopes'], reevingDict['elementLengths'],
116 cableTemplate, massProportionalLoad=gVec,
117 #fixedConstraintsNode0=[1,1,1,1], fixedConstraintsNode1=[1,1,1,1],
118 firstNodeIsLastNode=True, graphicsSizeConstraints=0.01)
119
120
121#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122#add contact:
123if useContact:
124
125 gContact = mbs.AddGeneralContact()
126 gContact.verboseMode = 1
127 gContact.frictionProportionalZone = 1
128 gContact.ancfCableUseExactMethod = False
129 gContact.ancfCableNumberOfContactSegments = 8
130 ssx = 16#32 #search tree size
131 ssy = 16#32 #search tree size
132 ssz = 1 #search tree size
133 gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
134 #gContact.SetSearchTreeBox(pMin=np.array([-1,-1,-1]), pMax=np.array([4,1,1])) #automatically computed!
135
136 halfHeight = 0.5*h*0
137 dimZ= 0.01 #for drawing
138 # wheels = [{'center':wheelCenter0, 'radius':rWheel0-halfHeight, 'mass':mWheel},
139 # {'center':wheelCenter1, 'radius':rWheel1-halfHeight, 'mass':mWheel},
140 # {'center':rollCenter1, 'radius':rRoll-halfHeight, 'mass':mRoll}, #small mass for roll, not to influence beam
141 # ]
142 sWheelRot = [] #sensors for angular velocity
143
144 nGround = mbs.AddNode(NodePointGround())
145 mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
146
147 for i, wheel in enumerate(circleList):
148 p = [wheel[0][0], wheel[0][1], 0] #position of wheel center
149 r = wheel[1]
150
151 rot0 = 0 #initial rotation
152 pRef = [p[0], p[1], rot0]
153 gList = [GraphicsDataCylinder(pAxis=[0,0,-dimZ],vAxis=[0,0,-dimZ], radius=r,
154 color= color4dodgerblue, nTiles=64),
155 GraphicsDataArrow(pAxis=[0,0,0], vAxis=[0.9*r,0,0], radius=0.01*r, color=color4orange)]
156
157 omega0 = 0 #initial angular velocity
158 v0 = np.array([0,0,omega0])
159
160 nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
161 visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
162 oMass = mbs.AddObject(ObjectRigidBody2D(physicsMass=wheelMass, physicsInertia=wheelInertia,
163 nodeNumber=nMass, visualization=
164 VObjectRigidBody2D(graphicsData=gList)))
165 mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
166 mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p))
167
168 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode]))
169 sWheelRot += [mbs.AddSensor(SensorNode(nodeNumber=nMass,
170 fileName='solution/wheel'+str(i)+'angVel.txt',
171 outputVariableType=exu.OutputVariableType.AngularVelocity))]
172
173 def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
174 v = 10*t
175 vMax = 5
176 return max(v, vMax)
177
178 velocityControl = True
179 if i == 0:
180 if velocityControl:
181 mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
182 mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
183 velocityLevel=True, offsetUserFunction_t=UFvelocityDrive))
184
185 # else:
186 # mbs.AddLoad(LoadTorqueVector(markerNumber=mNode, loadVector=[0,0,torque]))
187 if i > 0: #friction on rolls:
188 mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
189 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel],
190 damping=rotationDampingWheels,
191 visualization=VCoordinateSpringDamper(show=False)))
192
193 frictionMaterialIndex=0
194 gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=contactStiffness,
195 contactDamping=contactDamping, frictionMaterialIndex=frictionMaterialIndex)
196
197
198
199 for oIndex in ancf[1]:
200 gContact.AddANCFCable(objectIndex=oIndex, halfHeight=halfHeight, #halfHeight should be h/2, but then cylinders should be smaller
201 contactStiffness=contactStiffness, contactDamping=contactDamping,
202 frictionMaterialIndex=0)
203
204 frictionMatrix = np.zeros((2,2))
205 frictionMatrix[0,0]=int(useFriction)*dryFriction
206 frictionMatrix[0,1]=0 #no friction between some rolls and cable
207 frictionMatrix[1,0]=0 #no friction between some rolls and cable
208 gContact.SetFrictionPairings(frictionMatrix)
209
210
211mbs.Assemble()
212
213simulationSettings = exu.SimulationSettings() #takes currently set values or default values
214
215simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
216simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
217simulationSettings.solutionSettings.writeSolutionToFile = True
218simulationSettings.solutionSettings.solutionWritePeriod = 0.005
219simulationSettings.solutionSettings.sensorsWritePeriod = 0.001
220#simulationSettings.displayComputationTime = True
221simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
222simulationSettings.displayStatistics = True
223
224simulationSettings.timeIntegration.endTime = tEnd
225simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
226simulationSettings.timeIntegration.stepInformation= 3+128+256
227
228simulationSettings.timeIntegration.verboseMode = 1
229
230simulationSettings.timeIntegration.newton.useModifiedNewton = True
231simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
232simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
233simulationSettings.displayStatistics = True
234
235SC.visualizationSettings.general.circleTiling = 24
236SC.visualizationSettings.loads.show=False
237SC.visualizationSettings.nodes.defaultSize = 0.01
238SC.visualizationSettings.openGL.multiSampling = 4
239
240exu.SetWriteToConsole(False)
241
242if False:
243 SC.visualizationSettings.contour.outputVariableComponent=0
244 SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
245
246#visualize contact:
247if False:
248 SC.visualizationSettings.contact.showSearchTree =True
249 SC.visualizationSettings.contact.showSearchTreeCells =True
250 SC.visualizationSettings.contact.showBoundingBoxes = True
251
252if useGraphics:
253 exu.StartRenderer()
254 mbs.WaitForUserToContinue()
255
256doDynamic = True
257if doDynamic :
258 mbs.SolveDynamic(simulationSettings) #183 Newton iterations, 0.114 seconds
259else:
260 mbs.SolveStatic(simulationSettings) #183 Newton iterations, 0.114 seconds
261
262
263if useGraphics and True:
264 SC.visualizationSettings.general.autoFitScene = False
265 SC.visualizationSettings.general.graphicsUpdateInterval=0.02
266
267 sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
268 mbs.SolutionViewer(sol)
269
270
271if useGraphics:
272 SC.WaitForRenderEngineStopFlag()
273 exu.StopRenderer() #safely close rendering window!
274
275 # if True:
276 #
277 # mbs.PlotSensor(sensorNumbers=[sAngVel[0],sAngVel[1]], components=2, closeAll=True)
278 # mbs.PlotSensor(sensorNumbers=sMeasureRoll, components=1)