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