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