ANCFmovingRigidBodyTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for moving rigid body on two cables
  5#
  6# Author:   Andreas Zwölfer, Johannes Gerstmayr
  7# Date:     2019-12-16
  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 GenerateBeamElementsAlongLine
 16import exudyn.graphics as graphics
 17
 18testIsActive = exu.sys.get('testIsActive', False)
 19
 20SC = exu.SystemContainer()
 21mbs = SC.AddSystem()
 22
 23
 24#Import function  
 25#import GenerateStraightLineANCFCable2D as func
 26
 27#background
 28rect = [-2,-2,4,2] #xmin,ymin,xmax,ymax
 29background0 = {'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
 30background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 31oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 32
 33L=10
 34gravityFieldConstant=9.81 #9000
 35complianceFactBend = 0.1
 36complianceFactAxial = 0.1
 37nEl=10 #must be even number
 38vALE=2*5
 39offset=-1
 40 
 41
 42nGlobalGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) 
 43mGlobalGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGlobalGround, coordinate=0))
 44
 45fixANCFRotation = 1
 46
 47#######################SUSPENSION ROPE##################################################################################################################################################################
 48suspensionCableTemplate=Cable2D(massPerLength=20.87, 
 49                                bendingStiffness=78878*complianceFactBend, 
 50                                axialStiffness=398240000*complianceFactAxial)
 51
 52beamInfo=GenerateBeamElementsAlongLine(mbs=mbs, positionStart=[0,0,0], positionEnd=[L,0,0], numberOfElements=nEl, beamTemplate=suspensionCableTemplate,
 53                                                                  gravity=[0,-gravityFieldConstant,0], groundConstraintsStart=[1,1,0,fixANCFRotation], groundConstraintsEnd=[1,1,0,fixANCFRotation])
 54suspensionCableNodeList, suspensionCableObjectList, suspensionLoadList, suspensionCableNodePositionList, dummy = beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']
 55##################################################################################################################################################################
 56
 57 
 58
 59######################Haulage ROPE##################################################################################################################################################################
 60nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], initialCoordinates=[0], initialCoordinates_t=[vALE]))
 61mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity  marker
 62
 63haulageCableTemplate=ALECable2D(massPerLength=6.96, 
 64                                bendingStiffness=5956*complianceFactBend, 
 65                                axialStiffness=96725000*complianceFactAxial,
 66                                addALEvariation=False) #for compatibility with test suite results
 67haulageCableTemplate.nodeNumbers[2]=nALE #this will not be overwritten!
 68
 69beamInfo=GenerateBeamElementsAlongLine(mbs=mbs, 
 70                     positionStart=[0,offset,0], positionEnd=[L,offset,0], numberOfElements=nEl, beamTemplate=haulageCableTemplate,
 71                     gravity=[0,-gravityFieldConstant,0], groundConstraintsStart=[1,1,0,fixANCFRotation], groundConstraintsEnd=[1,1,0,fixANCFRotation])
 72haulageCableNodeList, haulageCableObjectList, haulageLoadList, haulageCableNodePositionList, dummy = beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']
 73
 74cAleConstraint=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mALE]))
 75##################################################################################################################################################################################################################################################
 76
 77
 78#slack carrier test
 79slackCarrierWheelRadius=0.75
 80mSuspensionRopeAttachmentNodeX=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = suspensionCableNodeList[int(nEl/2)], coordinate=0))
 81mSuspensionRopeAttachmentNodeY=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = suspensionCableNodeList[int(nEl/2)], coordinate=1))
 82
 83graphicsSlackCarrier={'type':'Circle', 'color':[.1,0.1,0.8,1], 'position':[0,0,0], 'radius': slackCarrierWheelRadius}
 84nSlackCarrierRigidBody = mbs.AddNode(Rigid2D(referenceCoordinates=[5,offset-slackCarrierWheelRadius,0]))
 85oSlackCarrierRigidBody = mbs.AddObject(RigidBody2D(mass=1, inertia=1, nodeNumber=nSlackCarrierRigidBody,visualization=VObjectRigidBody2D(graphicsData= [graphicsSlackCarrier]))) #, visualization=VObjectRigidBody2D(graphicsData= [graphicsSupportWheels])
 86
 87mSlackCarrierX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlackCarrierRigidBody,coordinate=0))
 88mSlackCarrierY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlackCarrierRigidBody,coordinate=1))
 89mSlackCarrierRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlackCarrierRigidBody,coordinate=2))
 90
 91mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeX,mSlackCarrierX])) 
 92mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeY,mSlackCarrierY])) 
 93mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mSlackCarrierRot]))
 94               
 95nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element
 96useFriction = False
 97nFactFriction = 1
 98if useFriction: nFactFriction = 3
 99
100initialGapList = [0.1]*(nSegments*nFactFriction) #initial gap of 0.1
101cStiffness = 1e7
102mContactSlackCarrier=mbs.AddMarker(MarkerBodyRigid(bodyNumber = oSlackCarrierRigidBody))  
103
104
105for i in haulageCableObjectList:
106    mContactCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=i, numberOfSegments = nSegments))  
107    nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments*nFactFriction)) 
108    if useFriction: 
109        mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mContactSlackCarrier, mContactCable], nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, contactStiffness = cStiffness, circleRadius = slackCarrierWheelRadius, offset = 0))  
110    else:
111        mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mContactSlackCarrier, mContactCable], nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, contactStiffness = cStiffness, circleRadius = slackCarrierWheelRadius, offset = 0))  
112      
113
114
115
116
117##################################################################################
118    
119a = 0.8     #y-dim/2 of gondula
120b = 0.02   #x-dim/2 of gondula
121yCOM = a    #COM distance to attachment point; in vertical direction
122massRigid = 60 #12
123#vInit=40
124inertiaRigid = massRigid/12*(2*a)**2
125g = 9.81    # gravity
126    
127refPos = [0,offset,0]
128#    refPos = [fieldData['stationData'][0]['referencePointCoordinates'][1][0],fieldData['maxVerticalPositionSuspensionRopeShoes'][0],0]
129
130#rigid body which slides:
131graphicsRigid1 = graphics.Lines([[-b,0,0], [b,0,0], [b,a,0], [-b,a,0], [-b,0,0]]) #drawing of rigid body
132graphicsRigid2 = graphics.Lines([[-a,-a,0], [a,-a,0], [a,0,0], [-a,0,0], [-a,-a,0]]) #drawing of rigid body
133
134nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[refPos[0],refPos[1]-yCOM,0], initialVelocities=[vALE,0,0]));
135oRigid = mbs.AddObject(RigidBody2D(mass=massRigid, inertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid1,graphicsRigid2])))
136
137
138markerRigidTopAle = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,yCOM,0.])) #support point
139mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
140mbs.AddLoad(Force(markerNumber=mR2, loadVector=[0,-massRigid*g,0]))
141
142aleCableMarkerList = []#list of Cable2DCoordinates markers
143aleOffsetList = []     #list of offsets counted from first cable element; needed in sliding joint
144aleOffset = 0          #first cable element has offset 0
145aleSlidingCoordinateInit = 0 
146aleInitialLocalMarker = 0 
147
148
149for i in range(len(haulageCableObjectList)): #create markers for cable elements
150    m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = haulageCableObjectList[i]))
151    aleCableMarkerList += [m] #list containing 'MarkerBodyCable2DCoordinates' marker for sliding joint
152    aleOffsetList += [aleOffset] #list of relative (arclength) coordinates of the starting point of a cable
153    aleOffset += L/nEl
154 
155nodeDataAleSlidingJoint = mbs.AddNode(NodeGenericData(initialCoordinates=[aleInitialLocalMarker],numberOfDataCoordinates=1)) #initial index in cable list   
156aleSlidingJoint = mbs.AddObject(ObjectJointALEMoving2D(name='aleSlider', markerNumbers=[markerRigidTopAle,aleCableMarkerList[aleInitialLocalMarker]], 
157                                                  slidingMarkerNumbers=aleCableMarkerList, slidingMarkerOffsets=aleOffsetList,nodeNumbers=[nodeDataAleSlidingJoint, nALE], activeConnector = False))
158 
159mnRigid0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=0)) #add rigid body marker
160mnRigid1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=1)) #add rigid body marker
161mnRigid2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=2)) #add rigid body marker
162nCCRigid0 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mnRigid0])) 
163nCCRigid1 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mnRigid1])) 
164nCCRigid2 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mnRigid2])) 
165
166#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
167# Assemble multibody system
168#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
169mbs.Assemble()
170#exu.Print(mbs)
171
172#SC.renderer.DoIdleTasks()
173
174#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
175# Simualtion settings:
176#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
177
178simulationSettings = exu.SimulationSettings()
179simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
180simulationSettings.staticSolver.numberOfLoadSteps  = 2
181 
182SC.visualizationSettings.general.circleTiling = 64
183SC.visualizationSettings.nodes.defaultSize=0.125
184
185SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
186SC.visualizationSettings.contour.outputVariableComponent = 1 # plot y-component
187
188
189SC.visualizationSettings.contact.contactPointsDefaultSize = .005
190SC.visualizationSettings.connectors.showContact = True
191 
192
193
194if not testIsActive: 
195    SC.renderer.Start()
196
197#get initial velocities
198vInit = mbs.systemData.GetODE2Coordinates_t(configuration = exu.ConfigurationType.Initial)
199
200#SC.renderer.DoIdleTasks()
201
202mbs.SolveStatic(simulationSettings) 
203
204#prolong solution for next computation
205u = mbs.systemData.GetODE2Coordinates()
206#v = mbs.systemData.GetODE2Coordinates_t()
207data = mbs.systemData.GetDataCoordinates()
208mbs.systemData.SetODE2Coordinates(u,configuration = exu.ConfigurationType.Initial)
209mbs.systemData.SetODE2Coordinates_t(vInit,configuration = exu.ConfigurationType.Initial)
210mbs.systemData.SetDataCoordinates(data,configuration = exu.ConfigurationType.Initial)
211
212#store some reference value:
213ncables = len(suspensionCableNodeList)
214sol = mbs.systemData.GetODE2Coordinates(); 
215u = sol[int(ncables/4)*4+1]; #y-displacement of node at midpoint of rope
216     
217
218#SC.renderer.DoIdleTasks()
219
220mbs.SetObjectParameter(aleSlidingJoint, 'activeConnector', True)
221mbs.SetObjectParameter(nCCRigid0, 'activeConnector', False)
222mbs.SetObjectParameter(nCCRigid1, 'activeConnector', False)
223mbs.SetObjectParameter(nCCRigid2, 'activeConnector', False)
224mbs.SetObjectParameter(cAleConstraint, 'activeConnector', False)
225
226
227solveDynamic = True
228if solveDynamic:
229    # time related settings:
230    steps=200
231    tend=0.1
232    h=tend/steps
233    
234    #fact = 15000
235    simulationSettings.timeIntegration.numberOfSteps = steps #1*fact
236    simulationSettings.timeIntegration.endTime = tend #0.002*fact
237    # Integrator related settings:
238    # simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
239    # simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
240    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.3
241    simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
242
243    simulationSettings.timeIntegration.verboseMode = 1
244    simulationSettings.timeIntegration.verboseModeFile = 0
245
246    
247    mbs.SolveDynamic(simulationSettings)
248    
249
250if not testIsActive: 
251    #SC.renderer.DoIdleTasks()
252    SC.renderer.Stop() #safely close rendering window!
253
254#compute error for test suite:
255ncables = len(suspensionCableNodeList)
256sol2 = mbs.systemData.GetODE2Coordinates(); 
257u2 = sol2[int(ncables/4)*4+1]; #y-displacement of node in first quater of rope
258exu.Print('static deflection  =',u)      #2020-03-05(corrected Cable2DshapeMarker): -0.06446474690480661    2019-12-17(new static solver): -0.06446474690512931;  2019-12-16: -0.06446474679809994
259exu.Print('dynamic deflection =',u2)       #2020-03-05(corrected Cable2DshapeMarker):0.06446627698400298; 2020-01-09: -0.06446627698121662(computeInitialAccelerations = False) 2020-01-09: -0.06446627843202835; 2019-12-26: -0.06446627698104967; 2019-12-17(update residual): -0.06446627698121662;  2019-12-16 (late): -0.06446627699890756; 2019-12-16: -0.06446610364603222
260exu.Print('ANCFmovingRigidBodyTest=',u+u2)       
261exu.sys['testResult'] = u + u2