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