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)