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