ANCFslidingJoint2D.py
You can view and download this file on Github: ANCFslidingJoint2D.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 = 3
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
97#mANCFend = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nodeList[-1], coordinate=1)) #last marker
98#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCFend]))
99
100#mGroundTip = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[L,0,0]))
101#mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundTip,markerList[-1]], stiffness=[10,10,10], damping=[0.1,0.1,0.1]))
102
103#mGroundTip2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[L,0.2,0]))
104#mbs.AddObject(SpringDamper(markerNumbers=[mGroundTip2,markerList[-1]], stiffness=0.1, referenceLength=0.2))
105
106#mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
107#mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -1e8, 0])) #will be changed in load steps
108
109#mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
110#mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*1*mypi]))
111
112#mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
113#mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 3*E*I*1*mypi]))
114
115a = 0.1 #y-dim/2 of gondula
116b = 0.001 #x-dim/2 of gondula
117massRigid = 12*0.01
118inertiaRigid = massRigid/12*(2*a)**2
119g = 9.81 # gravity
120
121slidingCoordinateInit = lElem*1.5 #0.75*L
122initialLocalMarker = 1 #second element
123if nElements<2:
124 slidingCoordinateInit /= 3.
125 initialLocalMarker = 0
126
127addRigidBody = True
128if addRigidBody:
129 #rigid body which slides:
130 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
131 nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[slidingCoordinateInit,-a,0], initialVelocities=[0,0,0]));
132 oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid])))
133
134 markerRigidTop = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,a,0.])) #support point
135 mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
136
137 mbs.AddLoad(Force(markerNumber = mR2, loadVector = [massRigid*g*0.1, -massRigid*g, 0]))
138
139
140
141#slidingJoint:
142addSlidingJoint = True
143if addSlidingJoint:
144 cableMarkerList = []#list of Cable2DCoordinates markers
145 offsetList = [] #list of offsets counted from first cable element; needed in sliding joint
146 offset = 0 #first cable element has offset 0
147 for item in cableList: #create markers for cable elements
148 m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = item))
149 cableMarkerList += [m]
150 offsetList += [offset]
151 offset += lElem
152
153 #mGroundSJ = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.*lElem+0.75*L,0.,0.]))
154 nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker,slidingCoordinateInit],numberOfDataCoordinates=2)) #initial index in cable list
155 slidingJoint = mbs.AddObject(ObjectJointSliding2D(name='slider', markerNumbers=[markerRigidTop,cableMarkerList[initialLocalMarker]],
156 slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList,
157 nodeNumber=nodeDataSJ))
158
159 #print(offsetList)
160
161
162
163#cStiffness = 1e3
164#cDamping = 0.02*cStiffness
165#useCircleContact = True
166#if useCircleContact:
167# nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element
168# initialGapList = [0.1]*nSegments #initial gap of 0.1
169
170# mGroundCircle = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.65*L,-0.5,0]))
171# mGroundCircle2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.25*L,-0.15,0]))
172
173# for i in range(len(cableList)):
174# #print("cable="+str(cableList[i]))
175# mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[i], numberOfSegments = nSegments))
176# #print("mCable="+str(mCable))
177# nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
178# mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle, mCable], nodeNumber = nodeDataContactCable,
179# numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
180# circleRadius = 0.3, offset = 0))
181# nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
182# mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle2, mCable], nodeNumber = nodeDataContactCable,
183# numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
184# circleRadius = 0.1, offset = 0))
185
186
187#mbs.systemData.Info()
188
189mbs.Assemble()
190print(mbs)
191
192simulationSettings = exu.SimulationSettings() #takes currently set values or default values
193#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
194
195fact = 1000
196simulationSettings.timeIntegration.numberOfSteps = 1*fact
197simulationSettings.timeIntegration.endTime = 0.001*fact*0.5
198simulationSettings.solutionSettings.writeSolutionToFile = True
199simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
200#simulationSettings.solutionSettings.outputPrecision = 4
201simulationSettings.displayComputationTime = True
202simulationSettings.timeIntegration.verboseMode = 1
203
204simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
205simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
206
207simulationSettings.timeIntegration.newton.useModifiedNewton = False
208simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
209simulationSettings.timeIntegration.newton.numericalDifferentiation.addReferenceCoordinatesToEpsilon = False
210simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1.e-3
211simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-8 #6.055454452393343e-06*0.0001 #eps^(1/3)
212simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1e8
213simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
214simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
215simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
216simulationSettings.pauseAfterEachStep = False
217simulationSettings.displayStatistics = True
218
219#SC.visualizationSettings.nodes.showNumbers = True
220SC.visualizationSettings.bodies.showNumbers = False
221#SC.visualizationSettings.connectors.showNumbers = True
222SC.visualizationSettings.nodes.defaultSize = 0.01
223SC.visualizationSettings.markers.defaultSize = 0.01
224SC.visualizationSettings.connectors.defaultSize = 0.01
225SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
226SC.visualizationSettings.connectors.showContact = 1
227
228simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
229
230solveDynamic = True
231if solveDynamic:
232 exu.StartRenderer()
233
234 mbs.SolveDynamic(simulationSettings)
235
236 SC.WaitForRenderEngineStopFlag()
237 exu.StopRenderer() #safely close rendering window!
238
239else:
240 simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-10*100 #can be quite small; WHY?
241 simulationSettings.staticSolver.verboseMode = 3
242 simulationSettings.staticSolver.numberOfLoadSteps = 20*2
243 simulationSettings.staticSolver.loadStepGeometric = False;
244 simulationSettings.staticSolver.loadStepGeometricRange = 5e3;
245
246 simulationSettings.staticSolver.newton.relativeTolerance = 1e-5*100 #10000
247 simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
248 simulationSettings.staticSolver.newton.maxIterations = 30 #50 for bending into circle
249
250 simulationSettings.staticSolver.discontinuous.iterationTolerance = 0.1
251 #simulationSettings.staticSolver.discontinuous.maxIterations = 5
252 simulationSettings.staticSolver.pauseAfterEachStep = False
253 simulationSettings.staticSolver.stabilizerODE2term = 100
254
255 exu.StartRenderer()
256
257 mbs.SolveStatic(simulationSettings)
258
259 #sol = mbs.systemData.GetODE2Coordinates()
260 #n = len(sol)
261 #print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
262
263 SC.WaitForRenderEngineStopFlag()
264 exu.StopRenderer() #safely close rendering window!
265
266# exu.InfoStat();
267
268
269
270
271#class MyDialog:
272# def __init__(self, parent):
273# top = self.top = Toplevel(parent)
274# Label(top, text="Value").pack()
275# self.e = Entry(top)
276# self.e.pack(padx=5)
277# b = Button(top, text="OK", command=self.ok)
278# b.pack(pady=5)
279# def ok(self):
280# #print("value is " + self.e.get())
281# exec(self.e.get())
282# self.top.destroy()
283
284#root = Tk()
285#Button(root, text="Exudyn").pack()
286#root.update()
287#d = MyDialog(root)
288#root.wait_window(d.top)