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