ANCFmovingRigidbody.py

You can view and download this file on Github: ANCFmovingRigidbody.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF ALE Cable2D test including a moving rigid body
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-10-15
  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 *
 15
 16SC = exu.SystemContainer()
 17mbs = SC.AddSystem()
 18
 19#err = ANCFCable2D_bending_test(df, SC, mbs)
 20#print(err)
 21
 22
 23#background
 24rect = [-2.5,-2,2.5,1] #xmin,ymin,xmax,ymax
 25background0 = {'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
 26background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 27oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0])))
 28
 29
 30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 31#cable:
 32mypi = 3.141592653589793
 33
 34L=2                     # length of ANCF element in m
 35#L=mypi                 # length of ANCF element in m
 36Em=2.07e11               # Young's modulus of ANCF element in N/m^2
 37rho=7800                # density of ANCF element in kg/m^3
 38b=0.002                 # width of rectangular ANCF element in m
 39h=0.002                 # height of rectangular ANCF element in m
 40A=b*h                   # cross sectional area of ANCF element in m^2
 41I=b*h**3/12             # second moment of area of ANCF element in m^4
 42EI = Em*I
 43rhoA = rho*A
 44EA = Em*A
 45movingMassFactor = 1    #1 = whole cable is moving with vALE speed
 46vALE = 1.3
 47g=9.81
 48
 49
 50print("L="+str(L))
 51print("EI="+str(EI))
 52print("EA="+str(EA))
 53print("rhoA="+str(rhoA))
 54
 55nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 56mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 57
 58cableList=[]        #for cable elements
 59nodeList=[]  #for nodes of cable
 60markerList=[]       #for nodes
 61
 62useALE = True
 63if useALE:
 64    nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], initialCoordinates=[0], initialCoordinates_t=[vALE]))
 65    mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity ==> must implement JacobianAE_t in CoordinateConstraint or similar
 66    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mALE], offset=vALE, velocityLevel = True)) # for static computation
 67
 68nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 69nodeList+=[nc0]
 70nElements = 32
 71lElem = L / nElements
 72for i in range(nElements):
 73    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 74    nodeList+=[nLast]
 75    if useALE:
 76        elem=mbs.AddObject(ALECable2D(physicsLength=lElem, physicsMassPerLength=rhoA, physicsBendingStiffness=EI, physicsAxialStiffness=EA, physicsMovingMassFactor=movingMassFactor, nodeNumbers=[nodeList[i],nodeList[i+1],nALE]))
 77    else:
 78        elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rhoA, physicsBendingStiffness=EI, physicsAxialStiffness=EA, nodeNumbers=[nc0+i,nc0+i+1]))
 79
 80    cableList+=[elem]
 81    mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber = elem))
 82    #mbs.AddLoad(Gravity(markerNumber=mBody, loadVector=[0,-g,0]))
 83
 84
 85
 86
 87mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=0))
 88mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=1))
 89mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=3))
 90
 91mANCF3 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=1)) #tip constraint
 92mANCF4 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=2)) #tip constraint
 93
 94mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 95mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 96mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 97mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF3]))
 98mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF4]))
 99
100#mANCF3 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=1))
101#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF3]))
102#mANCF4 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=0))
103#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF4]))
104
105
106a = 0.1     #y-dim/2 of gondula
107b = 0.001    #x-dim/2 of gondula
108massRigid = 12*0.01
109inertiaRigid = massRigid/12*(2*a)**2
110g = 9.81    # gravity
111
112slidingCoordinateInit = lElem*1.5*0 #0.75*L
113initialLocalMarker = 1 #second element
114if nElements<2:
115    slidingCoordinateInit /= 3.
116    initialLocalMarker = 0
117
118addRigidBody = True
119if addRigidBody:
120    #rigid body which slides:
121    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
122    nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[slidingCoordinateInit,-a,0], initialVelocities=[vALE,0,0]));
123    oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid])))
124
125    markerRigidTop = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,a,0.])) #support point
126    mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
127    mbs.AddLoad(Force(markerNumber = mR2, loadVector = [0, -massRigid*g, 0]))
128
129
130
131#slidingJoint:
132addSlidingJoint = False
133if addSlidingJoint:
134    cableMarkerList = []#list of Cable2DCoordinates markers
135    offsetList = []     #list of offsets counted from first cable element; needed in sliding joint
136    offset = 0          #first cable element has offset 0
137    for item in cableList: #create markers for cable elements
138        m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = item))
139        cableMarkerList += [m]
140        offsetList += [offset]
141        offset += lElem
142
143    #mGroundSJ = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.*lElem+0.75*L,0.,0.]))
144    nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker,slidingCoordinateInit],numberOfDataCoordinates=2)) #initial index in cable list
145    slidingJoint = mbs.AddObject(ObjectJointSliding2D(name='slider', markerNumbers=[markerRigidTop,cableMarkerList[initialLocalMarker]],
146                                                      slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList,
147                                                      nodeNumber=nodeDataSJ))
148#ALEslidingJoint:
149addALESlidingJoint = True
150if addALESlidingJoint:
151    cableMarkerList = []#list of Cable2DCoordinates markers
152    offsetList = []     #list of offsets counted from first cable element; needed in sliding joint
153    offset = 0          #first cable element has offset 0
154    for item in cableList: #create markers for cable elements
155        m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = item))
156        cableMarkerList += [m]
157        offsetList += [offset]
158        offset += lElem
159
160    #mGroundSJ = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.*lElem+0.75*L,0.,0.]))
161    nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker],numberOfDataCoordinates=1)) #initial index in cable list
162    #nodeODE2ALE = mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0],initialCoordinates=[0],initialCoordinates_t=[0],numberOfODE2Coordinates=1)) #initial index in cable list
163    slidingJoint = mbs.AddObject(ObjectJointALEMoving2D(name='slider', markerNumbers=[markerRigidTop,cableMarkerList[initialLocalMarker]],
164                                                      slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList, slidingOffset= -0*0.5*lElem,
165                                                      nodeNumbers=[nodeDataSJ, nALE]))
166
167
168    #print(offsetList)
169
170
171
172#cStiffness = 1e3
173#cDamping = 0.02*cStiffness
174#useCircleContact = True
175#if useCircleContact:
176#    nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element
177#    initialGapList = [0.1]*nSegments #initial gap of 0.1
178
179#    mGroundCircle = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.65*L,-0.5,0]))
180#    mGroundCircle2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.25*L,-0.15,0]))
181
182#    for i in range(len(cableList)):
183#        #print("cable="+str(cableList[i]))
184#        mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[i], numberOfSegments = nSegments))
185#        #print("mCable="+str(mCable))
186#        nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
187#        mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle, mCable], nodeNumber = nodeDataContactCable,
188#                                                 numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
189#                                                 circleRadius = 0.3, offset = 0))
190#        nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
191#        mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle2, mCable], nodeNumber = nodeDataContactCable,
192#                                                 numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
193#                                                 circleRadius = 0.1, offset = 0))
194
195
196#mbs.systemData.Info()
197
198mbs.Assemble()
199print(mbs)
200
201simulationSettings = exu.SimulationSettings() #takes currently set values or default values
202#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
203#simulationSettings.outputPrecision = 16
204
205fact = 1500
206simulationSettings.timeIntegration.numberOfSteps = fact
207simulationSettings.timeIntegration.endTime = 0.001*fact
208simulationSettings.solutionSettings.writeSolutionToFile = True
209simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/2000
210#simulationSettings.solutionSettings.outputPrecision = 4
211simulationSettings.displayComputationTime = True
212simulationSettings.timeIntegration.verboseMode = 1
213
214simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*10 #10000
215simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
216
217simulationSettings.timeIntegration.newton.useModifiedNewton = False
218simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 8
219simulationSettings.timeIntegration.newton.numericalDifferentiation.forAE = True #True should not be used in general, slow&inaccurate!
220simulationSettings.timeIntegration.newton.numericalDifferentiation.addReferenceCoordinatesToEpsilon = False
221simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1.e-3
222simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-8*10 #6.055454452393343e-06*0.0001 #eps^(1/3)
223simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1e8
224simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
225simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
226simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.4 #0.6 works well
227simulationSettings.pauseAfterEachStep = False
228simulationSettings.displayStatistics = True
229
230#SC.visualizationSettings.nodes.showNumbers = True
231SC.visualizationSettings.bodies.showNumbers = False
232#SC.visualizationSettings.connectors.showNumbers = True
233SC.visualizationSettings.nodes.defaultSize = 0.01
234SC.visualizationSettings.markers.defaultSize = 0.01
235SC.visualizationSettings.connectors.defaultSize = 0.01
236SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
237SC.visualizationSettings.connectors.showContact = 1
238
239simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
240
241solveDynamic = True
242if solveDynamic:
243    exu.StartRenderer()
244
245    mbs.SolveDynamic(simulationSettings)
246
247    SC.WaitForRenderEngineStopFlag()
248    exu.StopRenderer() #safely close rendering window!
249
250else:
251    simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-8 #*100 #can be quite small; WHY?
252    simulationSettings.staticSolver.verboseMode = 2
253    simulationSettings.staticSolver.numberOfLoadSteps  = 20#20*2
254    simulationSettings.staticSolver.loadStepGeometric = True;
255    simulationSettings.staticSolver.loadStepGeometricRange = 1e3;
256
257    simulationSettings.staticSolver.newton.relativeTolerance = 1e-5 #1e-5*100
258    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
259    simulationSettings.staticSolver.newton.maxIterations = 20 #50 for bending into circle
260
261    simulationSettings.staticSolver.discontinuous.iterationTolerance = 0.1
262    #simulationSettings.staticSolver.discontinuous.maxIterations = 5
263    simulationSettings.staticSolver.pauseAfterEachStep = False
264    simulationSettings.staticSolver.stabilizerODE2term = 100*0.0
265
266    exu.StartRenderer()
267
268    mbs.SolveStatic(simulationSettings)
269
270    sol = mbs.systemData.GetODE2Coordinates()
271    n = len(sol)
272    print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
273    sol_t = mbs.systemData.GetODE2Coordinates_t()
274    print('vALE='+str(sol_t[0]))
275
276    #print('sol='+str(sol))
277    print('sol_t='+str(sol_t))
278
279
280    SC.WaitForRenderEngineStopFlag()
281    exu.StopRenderer() #safely close rendering window!
282
283# exu.InfoStat();
284
285
286
287
288#class MyDialog:
289#    def __init__(self, parent):
290#        top = self.top = Toplevel(parent)
291#        Label(top, text="Value").pack()
292#        self.e = Entry(top)
293#        self.e.pack(padx=5)
294#        b = Button(top, text="OK", command=self.ok)
295#        b.pack(pady=5)
296#    def ok(self):
297#        #print("value is " + self.e.get())
298#        exec(self.e.get())
299#        self.top.destroy()
300
301#root = Tk()
302#Button(root, text="Exudyn").pack()
303#root.update()
304#d = MyDialog(root)
305#root.wait_window(d.top)