ANCFslidingAndALEjointTest.py

You can view and download this file on Github: ANCFslidingAndALEjointTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for rigid bodies sliding on cables
  5#
  6# Author:   Andreas Zwölfer, Johannes Gerstmayr
  7# Date:     2019-12-18
  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.utilities import * #includes itemInterface and rigidBodyUtilities
 15from exudyn.beams import GenerateAleSlidingJoint, GenerateSlidingJoint, GenerateBeamElementsAlongLine
 16import exudyn.graphics as graphics
 17
 18testIsActive = exu.sys.get('testIsActive', False)
 19
 20SC = exu.SystemContainer()
 21mbs = SC.AddSystem()
 22
 23
 24##################################################################################################################################################################
 25def AddBodyWithSlidingJoints(mbs,xPositionOfFirstNodes=0,referencePositionOfBodyAlongCable=0,gravityFieldConstant=0):
 26
 27    a = 0.8*2     #y-dim/2 of Body
 28    b = 0.02    #x-dim/2 of Body
 29    massRigid = 5000 #moving mass
 30    inertiaRigid = massRigid/12*(2*a)**2
 31    #rigid body which slides:
 32    graphicsRigid1 = graphics.Lines([[-b,-a,0], [b,-a,0], [b,a,0], [-b,a,0], [-b,-a,0]]) #drawing of rigid body
 33    yCOM = a    #COM distance to attachment point on suspension rope; in vertical direction      
 34    
 35    nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[xPositionOfFirstNodes+referencePositionOfBodyAlongCable,-yCOM,0]));
 36    oRigid = mbs.AddObject(RigidBody2D(mass=massRigid, inertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid1])))  
 37    markerRigidTop=mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,yCOM,0.])) #support point
 38    markerRigidTopAle = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,yCOM-offset,0.])) #support point
 39    
 40    if gravityFieldConstant:
 41        mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
 42        mbs.AddLoad(Force(markerNumber=mR2, loadVector=[0,-massRigid*gravityFieldConstant,0]))
 43    
 44    
 45    #find index of cable element associated with "referencePositionOfBodyAlongCable"
 46    cummulativeLength=0
 47    count=0
 48    while referencePositionOfBodyAlongCable >= cummulativeLength:
 49        cummulativeLength+=mbs.GetObjectParameter(cable2ObjectList[count],'length')
 50        count+=1
 51    count+=-1
 52    
 53    #AleSlidingJoint:
 54    oAleSlidingJoint=GenerateAleSlidingJoint(mbs,cable1ObjectList,markerRigidTopAle,AleNode=nALE,localMarkerIndexOfStartCable=count,AleSlidingOffset=referencePositionOfBodyAlongCable)[0]
 55      
 56    
 57    #slidingJoint:
 58    oSlidingJoint=GenerateSlidingJoint(mbs,cable2ObjectList,markerRigidTop,localMarkerIndexOfStartCable=count,slidingCoordinateStartPosition=referencePositionOfBodyAlongCable)[0] 
 59        
 60    return [oRigid,nRigid, oAleSlidingJoint, oSlidingJoint] 
 61
 62
 63#Background:
 64rect = [-2,-2,4,2] #xmin,ymin,xmax,ymax
 65background0 = {'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
 66background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 67oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 68
 69L=100  #length of ropes
 70gravityFieldConstant=9.81   
 71complianceFactBend = 1
 72complianceFactAxial = 1
 73nEl=6 #must be even number
 74vALE=0 
 75offset=0.35 #distance between the two ropes
 76offsetCarrier=-0.515 #distance between suspension rope and slack carrier wheel
 77
 78nGlobalGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) 
 79mGlobalGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGlobalGround, coordinate=0))
 80
 81fixANCFRotation = 0
 82
 83#######################ROPE2 (Carrier rope)########################################################################################################################
 84cable2Template=Cable2D(massPerLength=10, bendingStiffness=50000*complianceFactBend, axialStiffness=2e8*complianceFactAxial)
 85
 86beamInfo=GenerateBeamElementsAlongLine(mbs=mbs, 
 87                positionStart=[0,0,0], positionEnd=[L,0,0], numberOfElements=nEl, 
 88                beamTemplate=cable2Template, gravity=[0,-gravityFieldConstant,0], 
 89                groundConstraintsStart=[1,1,0,fixANCFRotation], groundConstraintsEnd=[1,1,0,fixANCFRotation])
 90cable2NodeList, cable2ObjectList, suspensionLoadList, cable2NodePositionList, dummy = beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']
 91##################################################################################################################################################################
 92
 93 
 94
 95######################ROPE1#######################################################################################################################################
 96nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], initialCoordinates=[0], initialCoordinates_t=[vALE]))
 97mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity  marker
 98
 99cable1Template=ALECable2D(massPerLength=3, bendingStiffness=4000*complianceFactBend, 
100                          axialStiffness=5e7*complianceFactAxial,useCouplingTerms=False,
101                          addALEvariation=False) #for compatibility with test suite results
102cable1Template.nodeNumbers[2]=nALE
103
104beamInfo=GenerateBeamElementsAlongLine(mbs=mbs, 
105            positionStart=[0,-offset,0], positionEnd=[L,-offset,0], numberOfElements=nEl, 
106            beamTemplate=cable1Template, gravity=[0,-gravityFieldConstant,0], 
107            groundConstraintsStart=[1,1,0,fixANCFRotation], groundConstraintsEnd=[1,1,0,fixANCFRotation])
108cable1NodeList, cable1ObjectList, haulageLoadList, cable1NodePositionList, dummy = beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']
109
110cAleConstraint=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mALE]))
111###################################################################################################################################################################
112
113
114#Slack carrier:
115carrierWheelRadius=0.42/2  
116mSuspensionRopeAttachmentNodeX=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = cable2NodeList[int(nEl/2)], coordinate=0))
117mSuspensionRopeAttachmentNodeY=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = cable2NodeList[int(nEl/2)], coordinate=1))
118
119graphicsCarrier={'type':'Circle', 'color':[.1,0.1,0.8,1], 'position':[0,0,0], 'radius': carrierWheelRadius}
120nCarrierRigidBody = mbs.AddNode(Rigid2D(referenceCoordinates=[L/2,offsetCarrier-carrierWheelRadius,0]))
121oCarrierRigidBody = mbs.AddObject(RigidBody2D(mass=200, inertia=1, 
122                                nodeNumber=nCarrierRigidBody,visualization=VObjectRigidBody2D(graphicsData= [graphicsCarrier])))  
123
124mCarrierX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nCarrierRigidBody,coordinate=0))
125mCarrierY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nCarrierRigidBody,coordinate=1))
126mCarrierRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nCarrierRigidBody,coordinate=2))
127
128mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeX,mCarrierX])) 
129mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeY,mCarrierY])) 
130mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mCarrierRot]))
131
132#mCarrierWheelLoad = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oCarrierRigidBody, localPosition=[ 0.,0.,0.])) #center of mass (for load)
133#mbs.AddLoad(Force(markerNumber=mConnectionWheelLoad, loadVector=[0,-CarrierMass*gravityFieldConstant,0]))
134               
135nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element
136useFriction = True
137nFactFriction = 1
138if useFriction: nFactFriction = 3
139
140initialGapList = [0.1]*(nSegments*nFactFriction) #initial gap of 0.1
141cStiffness = 1e7
142mContactCarrier=mbs.AddMarker(MarkerBodyRigid(bodyNumber = oCarrierRigidBody))  
143
144
145for i in cable1ObjectList:
146    mContactCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=i, numberOfSegments = nSegments))  
147    nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments*nFactFriction)) 
148    if useFriction: 
149        mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mContactCarrier, mContactCable], 
150                        nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, 
151                        contactStiffness = cStiffness, circleRadius = carrierWheelRadius))  
152    else:
153        mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mContactCarrier, mContactCable], 
154                    nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, 
155                    contactStiffness = cStiffness, circleRadius = carrierWheelRadius, offset = 0))  
156      
157
158
159#Add Bodys
160cRigidBodyRotList=[]
161cRigidBodyTranslationXlist=[]
162cRigidBodyTranslationYlist=[]
163distanceBetweenBodys=32
164numberOfBodys=2
165positionOfBody=0
166for i in range(numberOfBodys):
167
168    #Add Body
169    [oRigid,nRigid,oAleSlidingJoint,oSlidingJoint]=AddBodyWithSlidingJoints(mbs,xPositionOfFirstNodes=0,referencePositionOfBodyAlongCable=positionOfBody,gravityFieldConstant=gravityFieldConstant)
170    mbs.SetObjectParameter(oSlidingJoint, 'useClassicalFormulation', False) #test model computed with new sliding joint formulation
171
172    #fix rotation of rigid body
173    mRigidBodyRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=2)) #add rigid body marker
174    cRigidBodyRot = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mRigidBodyRot])) 
175    cRigidBodyRotList+=[cRigidBodyRot]
176    
177    #add damper for rotation of rigid body
178    mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGlobalGround,mRigidBodyRot],damping = 1e6, visualization=VObjectConnectorCoordinateSpringDamper(show=False)))  
179
180    positionOfBody+=distanceBetweenBodys
181    
182
183#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
184# Assemble multibody system
185#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
186mbs.Assemble()
187#exu.Print(mbs)
188
189#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
190# Simualtion settings:
191#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
192simulationSettings = exu.SimulationSettings()
193simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
194#simulationSettings.staticSolver.loadStepGeometric = True
195#simulationSettings.staticSolver.adaptiveStep = False
196simulationSettings.staticSolver.numberOfLoadSteps=10
197#simulationSettings.staticSolver.loadStepGeometricRange = 100
198simulationSettings.staticSolver.newton.relativeTolerance = 1e-7 #with this error tolerance, the adaptive step selection needs 4 steps
199
200#simulationSettings.staticSolver.verboseMode=1
201#simulationSettings.staticSolver.verboseModeFile=2
202simulationSettings.staticSolver.stabilizerODE2term = 2
203
204
205SC.visualizationSettings.general.circleTiling = 64
206SC.visualizationSettings.nodes.defaultSize=0.125
207SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
208SC.visualizationSettings.contour.outputVariableComponent = 1 # plot y-component
209SC.visualizationSettings.contact.contactPointsDefaultSize = .005
210SC.visualizationSettings.connectors.showContact = True
211
212
213if not testIsActive: 
214    SC.renderer.Start()
215
216#get initial velocities
217vInit = mbs.systemData.GetODE2Coordinates_t(configuration = exu.ConfigurationType.Initial)
218
219#start static calculation
220mbs.SolveStatic(simulationSettings) 
221
222#++++++++++++++++++++++++++++++++++++++++
223#compute error for test suite:
224ltgCable = mbs.systemData.GetObjectLTGODE2(cable1ObjectList[int(len(cable1ObjectList)/2)])
225nc = ltgCable[1] #vertical displacement
226exu.Print("select cable coordinate", nc)
227sol = mbs.systemData.GetODE2Coordinates(); 
228uStatic = sol[nc]; #y-displacement of first node of four bar mechanism
229exu.Print('static solution of cable1 =',uStatic)
230exu.sys['testResult'] = uStatic 
231
232#++++++++++++++++++++++++++++++++++++++++
233#store solution for next computation
234u = mbs.systemData.GetODE2Coordinates()
235data = mbs.systemData.GetDataCoordinates()
236
237#Add drive via ALE:
238def userLoadDriveAle(mbs, t, load):
239    #if t < 1:
240    return t*50000*5
241    #else: return 0
242
243mbs.AddLoad(LoadCoordinate(markerNumber = mALE, loadUserFunction=userLoadDriveAle))
244mbs.Assemble() #because of new load
245
246#++++++++++++++++++++++++++++++++++++++++
247#resuse old solution
248mbs.systemData.SetODE2Coordinates(u,configuration = exu.ConfigurationType.Initial)
249mbs.systemData.SetODE2Coordinates_t(vInit,configuration = exu.ConfigurationType.Initial)
250mbs.systemData.SetDataCoordinates(data,configuration = exu.ConfigurationType.Initial)
251
252
253    
254# remove rigid body motion constraints for dynamic analysis
255for item in cRigidBodyRotList:
256    mbs.SetObjectParameter(item, 'activeConnector', False)
257mbs.SetObjectParameter(cAleConstraint, 'activeConnector', False)
258
259
260solveDynamic = True
261if solveDynamic:
262    # time related settings:
263    steps=400
264    tend=0.8
265    h=tend/steps
266    SC.visualizationSettings.openGL.lineWidth = 2
267    simulationSettings.timeIntegration.numberOfSteps = steps 
268    simulationSettings.timeIntegration.endTime = tend 
269    simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
270    simulationSettings.timeIntegration.generalizedAlpha.useNewmark = simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints
271    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.3
272    simulationSettings.timeIntegration.verboseMode = 1
273    simulationSettings.show.statistics = True
274    
275    mbs.SolveDynamic(simulationSettings)
276    
277if not testIsActive: 
278    #SC.renderer.DoIdleTasks()
279    SC.renderer.Stop() #safely close rendering window!
280
281
282sol = mbs.systemData.GetODE2Coordinates(); 
283uDynamic = sol[nc]; #y-displacement of first node of four bar mechanism
284exu.Print('dynamic solution of cable1 =',uDynamic)
285
286exu.sys['testResult'] += uDynamic
287
288exu.Print('result of ANCFslidingAndALEjointTest=',exu.sys['testResult'])