ANCFcontactCircle2.py
You can view and download this file on Github: ANCFcontactCircle2.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: ANCF Cable2D contact test
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.utilities import *
15
16
17SC = exu.SystemContainer()
18mbs = SC.AddSystem()
19
20L=2 # length of ANCF element in m
21pCircle = [0.65*L,-0.5,0]
22pCircle2 = [0.25*L,-0.15,0]
23circleRadius=0.3
24circleRadius2=0.1
25
26#background
27rect = [-0.5,-1,2.5,1] #xmin,ymin,xmax,ymax
28background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
29
30
31background = [GraphicsDataRectangle(-0.5,-1,2.5,1, color=color4blue)]
32background += [GraphicsDataLine([[0,-1,0], [2,-1,0]], color=color4green)]
33background += [GraphicsDataCircle(point=pCircle, radius = circleRadius-0.002, color=color4blue)] #not necessary, as it is drawn by connector
34background += [GraphicsDataCircle(point=pCircle2, radius = circleRadius2-0.002, color=color4blue)] #not necessary, as it is drawn by connector
35background += [GraphicsDataText(point=[0.,0.2,0], text = 'ANCF contact with circle', color=color4black)]
36
37oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= background)))
38
39
40#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41#cable:
42mypi = 3.141592653589793
43
44#L=mypi # length of ANCF element in m
45E=2.07e11 # Young's modulus of ANCF element in N/m^2
46rho=7800 # density of ANCF element in kg/m^3
47b=0.001 # width of rectangular ANCF element in m
48h=0.001 # height of rectangular ANCF element in m
49A=b*h # cross sectional area of ANCF element in m^2
50I=b*h**3/12 # second moment of area of ANCF element in m^4
51f=3*E*I/L**2 # tip load applied to ANCF element in N
52
53print("load f="+str(f))
54print("EI="+str(E*I))
55
56nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
57mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
58
59cableList=[] #for cable elements
60nodeList=[] #for nodes of cable
61markerList=[] #for nodes
62nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
63nodeList+=[nc0]
64nElements = 8*4 #8,16,32,64
65lElem = L / nElements
66for i in range(nElements):
67 nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
68 nodeList+=[nLast]
69 elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
70 physicsBendingStiffness=E*I, physicsAxialStiffness=E*A*0.1,
71 nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
72 cableList+=[elem]
73
74mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
75mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
76mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
77
78mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
79mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
80mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
81
82#add gravity:
83markerList=[]
84for i in range(len(nodeList)):
85 m = mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeList[i]))
86 markerList+=[m]
87 fact = 1 #add (half) weight of two elements to node
88 if (i==0) | (i==len(nodeList)-1): fact = 0.5 # first and last node only weighted half
89 mbs.AddLoad(Force(markerNumber = m, loadVector = [0, -400*rho*A*fact*lElem, 0])) #will be changed in load steps
90
91#mANCFend = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nodeList[-1], coordinate=1)) #last marker
92#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCFend]))
93
94#mGroundTip = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[L,0,0]))
95#mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundTip,markerList[-1]], stiffness=[10,10,10], damping=[0.1,0.1,0.1]))
96
97#mGroundTip2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[L,0.2,0]))
98#mbs.AddObject(SpringDamper(markerNumbers=[mGroundTip2,markerList[-1]], stiffness=0.1, referenceLength=0.2))
99
100#mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
101#mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -1e8, 0])) #will be changed in load steps
102
103#mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
104#mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*1*mypi]))
105
106#mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
107#mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 3*E*I*1*mypi]))
108
109cStiffness = 1e3
110cDamping = 0.02*cStiffness*0
111useContact = False
112if useContact:
113 tipContact = False
114 if tipContact:
115 nodeData = mbs.AddNode(NodeGenericData(initialCoordinates=[0],numberOfDataCoordinates=1))
116 mbs.AddObject(ObjectContactCoordinate(markerNumbers=[mGround, mANCFend],nodeNumber = nodeData, contactStiffness = cStiffness, contactDamping=cDamping, offset = -0.8))
117 else:
118 for i in range(len(nodeList)):
119 mNC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nodeList[i], coordinate=1))
120 nodeData = mbs.AddNode(NodeGenericData(initialCoordinates=[1],numberOfDataCoordinates=1)) #start with gap!
121 mbs.AddObject(ObjectContactCoordinate(markerNumbers=[mGround, mNC], nodeNumber = nodeData, contactStiffness = cStiffness, contactDamping=cDamping, offset = -1))
122
123useCircleContact = True
124if useCircleContact:
125 nSegments = 8 #4; number of contact segments; must be consistent between nodedata and contact element
126 initialGapList = [0.1]*nSegments #initial gap of 0.1
127 mGroundCircle = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=pCircle))
128 mGroundCircle2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=pCircle2))
129
130 #mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=elem, numberOfSegments = nSegments))
131 #nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
132 #mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle, mCable], nodeNumber = nodeDataContactCable,
133 # numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
134 # circleRadius = 0.4, offset = 0))
135 for i in range(len(cableList)):
136 #print("cable="+str(cableList[i]))
137 mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[i], numberOfSegments = nSegments))
138 #print("mCable="+str(mCable))
139 nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
140 mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle, mCable], nodeNumber = nodeDataContactCable,
141 numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
142 circleRadius = circleRadius, offset = 0))
143 nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
144 mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle2, mCable], nodeNumber = nodeDataContactCable,
145 numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
146 circleRadius = circleRadius2, offset = 0))
147
148
149#mbs.systemData.Info()
150
151mbs.Assemble()
152print(mbs)
153
154simulationSettings = exu.SimulationSettings() #takes currently set values or default values
155#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
156
157fact = 10000
158simulationSettings.timeIntegration.numberOfSteps = 1*fact
159simulationSettings.timeIntegration.endTime = 0.001*fact
160simulationSettings.solutionSettings.writeSolutionToFile = True
161simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
162#simulationSettings.solutionSettings.outputPrecision = 4
163simulationSettings.displayComputationTime = True
164simulationSettings.timeIntegration.verboseMode = 1
165
166simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*10 #10000
167simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
168
169simulationSettings.timeIntegration.newton.useModifiedNewton = False
170simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
171simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
172simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 6.055454452393343e-06*0.1 #eps^(1/3)
173simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1e8
174simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
175simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
176simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
177simulationSettings.displayStatistics = True
178
179#SC.visualizationSettings.nodes.showNumbers = True
180SC.visualizationSettings.bodies.showNumbers = False
181#SC.visualizationSettings.connectors.showNumbers = True
182SC.visualizationSettings.nodes.defaultSize = 0.01
183SC.visualizationSettings.markers.defaultSize = 0.01
184SC.visualizationSettings.connectors.defaultSize = 0.01
185SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
186SC.visualizationSettings.connectors.showContact = 1
187SC.visualizationSettings.general.circleTiling = 64
188
189simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
190
191solveDynamic = False
192if solveDynamic:
193 exu.StartRenderer()
194
195 mbs.SolveDynamic(simulationSettings)
196
197 SC.WaitForRenderEngineStopFlag()
198 exu.StopRenderer() #safely close rendering window!
199
200else:
201 simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-10*100 #can be quite small; WHY?
202 simulationSettings.staticSolver.newton.numericalDifferentiation.doSystemWideDifferentiation = False
203 simulationSettings.staticSolver.verboseMode = 1
204 simulationSettings.staticSolver.numberOfLoadSteps = 20*2
205 simulationSettings.staticSolver.loadStepGeometric = True;
206 simulationSettings.staticSolver.loadStepGeometricRange = 5e3;
207
208 simulationSettings.staticSolver.newton.relativeTolerance = 1e-5*100 #10000
209 simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
210 simulationSettings.staticSolver.newton.maxIterations = 30 #50 for bending into circle
211
212 simulationSettings.staticSolver.discontinuous.iterationTolerance = 0.1
213 #simulationSettings.staticSolver.discontinuous.maxIterations = 5
214 simulationSettings.pauseAfterEachStep = False
215 simulationSettings.staticSolver.stabilizerODE2term = 50
216
217 exu.StartRenderer()
218
219 mbs.WaitForUserToContinue()
220 mbs.SolveStatic(simulationSettings)
221
222 #sol = mbs.systemData.GetODE2Coordinates()
223 #n = len(sol)
224 #print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
225
226 SC.WaitForRenderEngineStopFlag()
227 exu.StopRenderer() #safely close rendering window!
228
229# exu.InfoStat();