ANCFslidingJoint2Drigid.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D element with sliding joint test
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-09-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#background
 20rect = [-2.5,-2,2.5,1] #xmin,ymin,xmax,ymax
 21background0 = {'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
 22background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 23oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0])))
 24
 25
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27#cable:
 28mypi = 3.141592653589793
 29
 30L=2                     # length of ANCF element in m
 31#L=mypi                 # length of ANCF element in m
 32E=2.07e11               # Young's modulus of ANCF element in N/m^2
 33rho=7800                # density of ANCF element in kg/m^3
 34b=0.001                 # width of rectangular ANCF element in m
 35h=0.001                 # height of rectangular ANCF element in m
 36A=b*h                   # cross sectional area of ANCF element in m^2
 37I=b*h**3/12             # second moment of area of ANCF element in m^4
 38f=3*E*I/L**2            # tip load applied to ANCF element in N
 39g=9.81
 40
 41print("load f="+str(f))
 42print("EI="+str(E*I))
 43
 44nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 45mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 46
 47cableList=[]        #for cable elements
 48nodeList=[]  #for nodes of cable
 49markerList=[]       #for nodes
 50nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 51nodeList+=[nc0]
 52nElements = 32
 53lElem = L / nElements
 54for i in range(nElements):
 55    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 56    nodeList+=[nLast]
 57    elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
 58                               physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 59    cableList+=[elem]
 60    mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber = elem))
 61    mbs.AddLoad(Gravity(markerNumber=mBody, loadVector=[0,-g,0]))
 62
 63addPointMass = False
 64if addPointMass:
 65    massTip = 0.01 #tip mass
 66    nMass = mbs.AddNode(Point2D(referenceCoordinates=[L,0],visualization=VNodePoint2D(drawSize=0.3)))
 67    mTip0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 68    mTip1 = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast))
 69    mbs.AddObject(MassPoint2D(physicsMass = massTip, nodeNumber=nMass))
 70    mbs.AddLoad(Force(markerNumber=mTip0, loadVector=[0,-massTip*g,0]))
 71    mbs.AddObject(RevoluteJoint2D(markerNumbers=[mTip0,mTip1]))
 72
 73
 74mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=0))
 75mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=1))
 76mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=3))
 77
 78mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 79mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 80mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 81
 82#mANCF3 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=1))
 83#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF3]))
 84#mANCF4 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=0))
 85#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF4]))
 86
 87#add gravity:
 88markerList=[]
 89for i in range(len(nodeList)):
 90    m = mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeList[i]))
 91    markerList+=[m]
 92    #fact = 1 #add (half) weight of two elements to node
 93    #if (i==0) | (i==len(nodeList)-1):
 94    #    fact = 0.5 # first and last node only weighted half
 95    #mbs.AddLoad(Force(markerNumber = m, loadVector = [0., -rho*A*fact*lElem*g, 0])) #will be changed in load steps
 96
 97a = 0.1     #y-dim/2 of gondula
 98b = 0.001    #x-dim/2 of gondula
 99massRigid = 12*0.01
100inertiaRigid = massRigid/12*(2*a)**2
101g = 9.81    # gravity
102
103slidingCoordinateInit = lElem*1.5 #0.75*L
104initialLocalMarker = 1 #second element
105if nElements<2:
106    slidingCoordinateInit /= 3.
107    initialLocalMarker = 0
108
109addRigidBody = True
110if addRigidBody:
111    #rigid body which slides:
112    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
113    nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[slidingCoordinateInit,-a,0], initialVelocities=[0,0,0]));
114    oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid])))
115
116    markerRigidTop = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,a,0.])) #support point
117    mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
118
119    mbs.AddLoad(Force(markerNumber = mR2, loadVector = [massRigid*g*0.1, -massRigid*g, 0]))
120
121
122
123#slidingJoint:
124addSlidingJoint = True
125if addSlidingJoint:
126    cableMarkerList = []#list of Cable2DCoordinates markers
127    offsetList = []     #list of offsets counted from first cable element; needed in sliding joint
128    offset = 0          #first cable element has offset 0
129    for item in cableList: #create markers for cable elements
130        m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = item))
131        cableMarkerList += [m]
132        offsetList += [offset]
133        offset += lElem
134
135    #mGroundSJ = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.*lElem+0.75*L,0.,0.]))
136    nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker,slidingCoordinateInit],numberOfDataCoordinates=2)) #initial index in cable list
137    slidingJoint = mbs.AddObject(ObjectJointSliding2D(name='slider', markerNumbers=[markerRigidTop,cableMarkerList[initialLocalMarker]],
138                                                      slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList,
139                                                      nodeNumber=nodeDataSJ))
140
141
142mbs.Assemble()
143print(mbs)
144
145simulationSettings = exu.SimulationSettings() #takes currently set values or default values
146#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
147
148h=5e-4
149tEnd = 0.6
150simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
151simulationSettings.timeIntegration.endTime = tEnd
152simulationSettings.solutionSettings.writeSolutionToFile = True
153simulationSettings.solutionSettings.solutionWritePeriod = h
154#simulationSettings.solutionSettings.outputPrecision = 4
155simulationSettings.displayComputationTime = True
156simulationSettings.timeIntegration.verboseMode = 1
157
158# simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
159# simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
160
161simulationSettings.timeIntegration.newton.useModifiedNewton = True
162simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
163simulationSettings.pauseAfterEachStep = False
164simulationSettings.displayStatistics = True
165
166#SC.visualizationSettings.nodes.showNumbers = True
167SC.visualizationSettings.bodies.showNumbers = False
168SC.visualizationSettings.loads.show = False
169#SC.visualizationSettings.connectors.showNumbers = True
170SC.visualizationSettings.nodes.defaultSize = 0.01
171SC.visualizationSettings.markers.defaultSize = 0.01
172SC.visualizationSettings.connectors.defaultSize = 0.01
173SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
174SC.visualizationSettings.connectors.showContact = 1
175
176simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
177
178solveDynamic = True
179if solveDynamic:
180    exu.StartRenderer()
181    mbs.WaitForUserToContinue()
182
183    mbs.SolveDynamic(simulationSettings)
184
185    SC.WaitForRenderEngineStopFlag()
186    exu.StopRenderer() #safely close rendering window!
187
188else:
189    simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-10*100 #can be quite small; WHY?
190    simulationSettings.staticSolver.newton.numericalDifferentiation.doSystemWideDifferentiation = False
191    simulationSettings.staticSolver.newton.useNumericalDifferentiation = False
192    simulationSettings.staticSolver.verboseMode = 3
193    simulationSettings.staticSolver.numberOfLoadSteps  = 20*2
194    simulationSettings.staticSolver.loadStepGeometric = False;
195    simulationSettings.staticSolver.loadStepGeometricRange = 5e3;
196
197    simulationSettings.staticSolver.newton.relativeTolerance = 1e-5*100 #10000
198    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
199    simulationSettings.staticSolver.newton.maxIterations = 30 #50 for bending into circle
200
201    simulationSettings.staticSolver.discontinuous.iterationTolerance = 0.1
202    #simulationSettings.staticSolver.discontinuous.maxIterations = 5
203    simulationSettings.staticSolver.pauseAfterEachStep = False
204    simulationSettings.staticSolver.stabilizerODE2term = 100
205
206    exu.StartRenderer()
207
208    mbs.SolveStatic(simulationSettings)
209
210    #sol = mbs.systemData.GetODE2Coordinates()
211    #n = len(sol)
212    #print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
213
214    SC.WaitForRenderEngineStopFlag()
215    exu.StopRenderer() #safely close rendering window!
216
217# exu.InfoStat();
218
219
220
221
222#class MyDialog:
223#    def __init__(self, parent):
224#        top = self.top = Toplevel(parent)
225#        Label(top, text="Value").pack()
226#        self.e = Entry(top)
227#        self.e.pack(padx=5)
228#        b = Button(top, text="OK", command=self.ok)
229#        b.pack(pady=5)
230#    def ok(self):
231#        #print("value is " + self.e.get())
232#        exec(self.e.get())
233#        self.top.destroy()
234
235#root = Tk()
236#Button(root, text="Exudyn").pack()
237#root.update()
238#d = MyDialog(root)
239#root.wait_window(d.top)