ANCFswitchingSlidingJoint2D.py
You can view and download this file on Github: ANCFswitchingSlidingJoint2D.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: ANCF Cable2D element with SlidingJoint2D; after the object reaches a certain position, it is reset to the origin
5#
6# Author: Johannes Gerstmayr
7# Date: 2019-10-01
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.itemInterface import *
15from exudyn.utilities import *
16
17SC = exu.SystemContainer()
18mbs = SC.AddSystem()
19
20print('EXUDYN version='+exu.GetVersionString())
21
22#testInterface = TestInterface(exudyn = exu, systemContainer = SC, useGraphics=False)
23#RunAllModelUnitTests(mbs, testInterface)
24
25
26SC = exu.SystemContainer()
27mbs = SC.AddSystem()
28
29background = GraphicsDataRectangle(-0.1,-1.5,2.5,0.25, color=[0.9,0.9,0.9,1.])
30oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
31#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32#cable:
33mypi = 3.141592653589793
34
35L=2 # length of ANCF element in m
36#L=mypi # length of ANCF element in m
37E=2.07e11*0.2 # Young's modulus of ANCF element in N/m^2
38rho=7800 # density of ANCF element in kg/m^3
39b=0.001 # width of rectangular ANCF element in m
40h=0.001 # height of rectangular ANCF element in m
41A=b*h # cross sectional area of ANCF element in m^2
42I=b*h**3/12 # second moment of area of ANCF element in m^4
43f=3*E*I/L**2 # tip load applied to ANCF element in N
44g=9.81
45
46print("load f="+str(f))
47print("EI="+str(E*I))
48
49nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
50mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
51
52cableList=[] #for cable elements
53nodeList=[] #for nodes of cable
54nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
55nodeList+=[nc0]
56nElements = 8*32
57lElem = L / nElements
58for i in range(nElements):
59 nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
60 nodeList+=[nLast]
61 elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
62 physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
63 cableList+=[elem]
64 mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber = elem))
65 mbs.AddLoad(Gravity(markerNumber=mBody, loadVector=[0,-g,0]))
66
67
68mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=0))
69mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=1))
70mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=3))
71
72mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
73mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
74mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
75
76nGroundTip = mbs.AddNode(NodePointGround(referenceCoordinates=[L,0,0])) #ground node for coordinate constraint
77mGroundTip = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGroundTip, coordinate=0)) #Ground node ==> no action
78
79mANCF3 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=0))
80#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF3]))
81k=1e3
82mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGroundTip,mANCF3], stiffness = k, damping = k*0.02))
83
84mANCF4 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=1))
85mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF4]))
86mANCF5 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=2))
87mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF5]))
88
89a = 0.1 #y-dim/2 of gondula
90b = 0.001 #x-dim/2 of gondula
91massRigid = 12*0.01
92inertiaRigid = massRigid/12*(2*a)**2
93g = 9.81 # gravity
94
95slidingCoordinateInit = 0*0.25*lElem #0*lElem*1.5 #0.75*L
96initialLocalMarker = 0 #1 .. second element
97if nElements<2:
98 slidingCoordinateInit /= 3.
99 initialLocalMarker = 0
100
101addRigidBody = True
102if addRigidBody:
103 vSliding = 2
104 #rigid body which slides:
105 graphicsRigid = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-b,-a,0, b,-a,0, b,a,0, -b,a,0, -b,-a,0]} #drawing of rigid body
106 nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[slidingCoordinateInit,-a,0], initialVelocities=[vSliding,0,0]));
107 oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid])))
108
109 markerRigidTop = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,a,0.])) #support point
110 mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
111
112 #constant velocity driving:
113 mNCRigid = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=0)) #BaseException-coordinate
114 mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNCRigid], velocityLevel = True, offset = vSliding))
115
116
117 #mbs.AddLoad(Force(markerNumber = mR2, loadVector = [massRigid*g*0.1, -massRigid*g, 0]))
118
119
120#slidingJoint:
121addSlidingJoint = True
122if addSlidingJoint:
123 cableMarkerList = []#list of Cable2DCoordinates markers
124 offsetList = [] #list of offsets counted from first cable element; needed in sliding joint
125 offset = 0 #first cable element has offset 0
126 for item in cableList: #create markers for cable elements
127 m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = item))
128 cableMarkerList += [m]
129 offsetList += [offset]
130 offset += lElem
131
132 nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker,slidingCoordinateInit],numberOfDataCoordinates=2)) #initial index in cable list
133 slidingJoint = mbs.AddObject(ObjectJointSliding2D(name='slider', markerNumbers=[markerRigidTop,cableMarkerList[initialLocalMarker]],
134 slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList,
135 nodeNumber=nodeDataSJ))
136
137
138mbs.Assemble()
139print(mbs)
140
141simulationSettings = exu.SimulationSettings() #takes currently set values or default values
142#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
143
144
145
146fact = 2000
147deltaT = 0.0005*fact
148simulationSettings.timeIntegration.numberOfSteps = 1*fact
149simulationSettings.timeIntegration.endTime = deltaT
150simulationSettings.solutionSettings.writeSolutionToFile = True
151simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
152#simulationSettings.solutionSettings.outputPrecision = 4
153simulationSettings.displayComputationTime = False
154simulationSettings.timeIntegration.verboseMode = 1
155
156simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
157simulationSettings.timeIntegration.newton.useModifiedNewton = True
158simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
159simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-5
160simulationSettings.timeIntegration.discontinuous.maxIterations = 2 #only two for selection of correct sliding cable element
161
162useIndex2 = False
163simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = useIndex2
164simulationSettings.timeIntegration.generalizedAlpha.useNewmark = useIndex2
165simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
166simulationSettings.displayStatistics = False
167simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
168
169#SC.visualizationSettings.nodes.showNumbers = True
170SC.visualizationSettings.bodies.showNumbers = False
171SC.visualizationSettings.loads.show = False
172#SC.visualizationSettings.connectors.showNumbers = True
173SC.visualizationSettings.nodes.defaultSize = 0.002
174SC.visualizationSettings.markers.defaultSize = 0.002
175SC.visualizationSettings.connectors.defaultSize = 0.01
176SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
177SC.visualizationSettings.connectors.showContact = 1
178
179#SC.visualizationSettings.general.minSceneSize = 4
180SC.visualizationSettings.openGL.initialCenterPoint = [0.5*L,-0.25*L,0]
181#SC.visualizationSettings.openGL.lineWidth=2
182
183simulationSettings.solutionSettings.solutionInformation = "ANCF cable with sliding joint"
184
185#mbs.systemData.Info()
186
187
188def gondulaReset(oRigid, oSlidingJoint, maxL, vSliding):
189 u = mbs.GetObjectOutput(oSlidingJoint, exu.OutputVariableType.SlidingCoordinate)
190
191 if u > maxL: #reset rigid body to start of rope
192 print('active connector = ', mbs.GetObjectParameter(slidingJoint, 'activeConnector'))
193 coordsODE2 = mbs.systemData.GetODE2Coordinates()
194 coordsODE2_t = mbs.systemData.GetODE2Coordinates_t()
195 coordsData = mbs.systemData.GetDataCoordinates()
196
197 LTG = mbs.systemData.GetObjectLTGODE2(oRigid)
198 LTGdata = mbs.systemData.GetObjectLTGData(oSlidingJoint)
199
200 #set new data coordinates:
201 coordsODE2[LTG[0]] = 0
202 coordsODE2[LTG[1]] = 0
203 coordsODE2[LTG[2]] = 0
204 coordsODE2_t[LTG[0]] = vSliding
205 coordsODE2_t[LTG[1]] = 0
206 coordsODE2_t[LTG[2]] = 0
207 coordsData[LTGdata[0]] = 0 #initial sliding marker index
208 coordsData[LTGdata[1]] = 0 #initial (start of step) sliding coordinate
209
210 #fill into system coordinates:
211 mbs.systemData.SetODE2Coordinates(coordsODE2)
212 mbs.systemData.SetODE2Coordinates_t(coordsODE2_t)
213 mbs.systemData.SetDataCoordinates(coordsData)
214 mbs.systemData.SetDataCoordinates(coordsData, configuration=exu.ConfigurationType.StartOfStep)
215
216
217maxL = 0.9999*L
218
219#new user function executed at every beginning of time steps
220def UFgondulaReset(mbs, t):
221 gondulaReset(oRigid, slidingJoint, maxL, vSliding)
222 return True #True, means that everything is alright, False=stop simulation
223
224mbs.SetPreStepUserFunction(UFgondulaReset)
225
226
227exu.StartRenderer()
228mbs.SolveDynamic(simulationSettings)
229
230if False:
231 for i in range(5000): #2500
232 mbs.SolveDynamic(simulationSettings)
233
234 if mbs.GetRenderEngineStopFlag():
235 print('stopped by user')
236 break
237
238 u = mbs.GetObjectOutput(slidingJoint, exu.OutputVariableType.SlidingCoordinate)
239 #print('STEP ',i, ', t =', i*deltaT, ', sliding coordinate =',u)
240
241 coordsODE2 = mbs.systemData.GetODE2Coordinates()
242 coordsODE2_t = mbs.systemData.GetODE2Coordinates_t()
243 coordsAE = mbs.systemData.GetAECoordinates()
244 coordsData = mbs.systemData.GetDataCoordinates()
245 LTG = mbs.systemData.GetObjectLTGODE2(oRigid)
246 LTGAE = mbs.systemData.GetObjectLTGAE(slidingJoint)
247 LTGdata = mbs.systemData.GetObjectLTGData(slidingJoint)
248
249 if i*deltaT > 10:
250 print('coordsODE2 =', coordsODE2[LTG[0:3]])
251 print('coordsODE2_t=', coordsODE2_t[LTG[0:3]])
252 print('coordsAE =', coordsAE[LTGAE[0:3]])
253 print('coordsData =', coordsData[LTGdata[0:2]])
254
255 if u > 0.99*L: #reset rigid body to start of rope
256 print('active connector = ', mbs.GetObjectParameter(slidingJoint, 'activeConnector'))
257 #simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
258 #simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
259 #mbs.SetObjectParameter(slidingJoint, 'activeConnector', False)
260 #set parameters back to origin
261 coordsODE2[LTG[0]] = 0
262 coordsODE2[LTG[1]] = 0
263 coordsODE2[LTG[2]] = 0
264 coordsODE2_t[LTG[0]] = vSliding
265 coordsODE2_t[LTG[1]] = 0
266 coordsODE2_t[LTG[2]] = 0
267 coordsData[LTGdata[0]] = 0 #initial sliding marker index
268 coordsData[LTGdata[1]] = 0 #initial (start of step) sliding coordinate
269 mbs.systemData.SetDataCoordinates(coordsData,configuration = exu.ConfigurationType.Current) #is used as startOfStep for next step
270 #mbs.WaitForUserToContinue()
271
272
273
274 mbs.systemData.SetODE2Coordinates(coordsODE2,configuration = exu.ConfigurationType.Initial)
275 mbs.systemData.SetODE2Coordinates_t(coordsODE2_t,configuration = exu.ConfigurationType.Initial)
276 mbs.systemData.SetDataCoordinates(coordsData,configuration = exu.ConfigurationType.Initial)
277 mbs.systemData.SetAECoordinates(coordsAE,configuration = exu.ConfigurationType.Initial)
278
279
280SC.WaitForRenderEngineStopFlag()
281exu.StopRenderer() #safely close rendering window!