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!