ANCFcontactFrictionTest.py
You can view and download this file on Github: ANCFcontactFrictionTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test model for cable with contact and friction; test model for ObjectContactFrictionCircleCable2D,
5# which models frictional contact between 2D ANCF element and circular object, using contact and stick-slip friction;
6#
7# Author: Johannes Gerstmayr
8# Date: 2019-08-15 (created)
9# Date: 2022-07-20 (last modified)
10#
11# 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.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14
15import exudyn as exu
16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
17import exudyn.graphics as graphics #only import if it does not conflict
18
19useGraphics = True #without test
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
22try: #only if called from test suite
23 from modelUnitTests import exudynTestGlobals #for globally storing test results
24 useGraphics = exudynTestGlobals.useGraphics
25except:
26 class ExudynTestGlobals:
27 pass
28 exudynTestGlobals = ExudynTestGlobals()
29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
30
31SC = exu.SystemContainer()
32mbs = SC.AddSystem()
33
34
35#background
36rect = [-0.5,-1,2.5,1] #xmin,ymin,xmax,ymax
37background0 = {'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
38background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
39oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0, background1])))
40
41#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
42#cable:
43mypi = 3.141592653589793
44
45L=2 # length of ANCF element in m
46#L=mypi # length of ANCF element in m
47E=2.07e11 # Young's modulus of ANCF element in N/m^2
48rho=7800 # density of ANCF element in kg/m^3
49b=0.001*10 # width of rectangular ANCF element in m
50h=0.001*10 # height of rectangular ANCF element in m
51A=b*h # cross sectional area of ANCF element in m^2
52I=b*h**3/12 # second moment of area of ANCF element in m^4
53f=3*E*I/L**2 # tip load applied to ANCF element in N
54
55exu.Print("load f="+str(f))
56exu.Print("EI="+str(E*I))
57
58nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
59mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
60
61cableList=[] #for cable elements
62nodeList=[] #for nodes of cable
63markerList=[] #for nodes
64
65#%%+++++++++++++++++++++
66#create nodes and cable elements;
67#alternatively, GenerateStraightLineANCFCable2D from exudyn.beams could be used
68nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
69nodeList+=[nc0]
70nElements = 8 #32*4
71lElem = L / nElements
72for i in range(nElements):
73 nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
74 nodeList+=[nLast]
75 elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
76 physicsBendingStiffness=E*I, physicsAxialStiffness=E*A*0.1,
77 physicsBendingDamping=E*I*0.025*0, physicsAxialDamping=E*A*0.1,
78 nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
79 cableList+=[elem]
80
81#%%+++++++++++++++++++++
82mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
83mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
84mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
85
86mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
87mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
88mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
89
90#add gravity:
91markerList=[]
92for i in range(len(nodeList)):
93 m = mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeList[i]))
94 markerList+=[m]
95 fact = 1 #add (half) weight of two elements to node
96 if (i==0) | (i==len(nodeList)-1): fact = 0.5 # first and last node only weighted half
97 mbs.AddLoad(Force(markerNumber = m, loadVector = [0, -0.1*400*rho*A*fact*lElem, 0])) #will be changed in load steps
98
99
100cStiffness = 1e3
101cDamping = 0.02*cStiffness*2
102r1 = 0.1
103r2 = 0.3
104
105posRoll1 = [0.25*L,-0.15,0]
106posRoll2 = [0.75*L,-0.5,0]
107
108useFriction = 1
109useCircleContact = True
110if useCircleContact:
111 nSegments = 2*4 #4; number of contact segments; must be consistent between nodedata and contact element
112 initialGapList = [0.1]*nSegments #initial gap of 0.1
113 if (useFriction):
114 initialGapList += [0]*(2*nSegments)
115
116 mGroundCircle = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=posRoll1))
117 mGroundCircle2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=posRoll2))
118
119
120 rGraphics = GraphicsDataRectangle(0.,0.,0.1*r2,r2)
121 vRigidBody = VObjectRigidBody2D(graphicsData = [rGraphics])
122 nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=posRoll2))
123 oRigid = mbs.AddObject(RigidBody2D(nodeNumber = nRigid, physicsMass = 1, physicsInertia=0.001, visualization=vRigidBody))
124 mRigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRigid, localPosition=[0,0,0]))
125 # mRigid = mbs.AddMarker(MarkerNodeRigid(nodeNumber = nRigid)) #gives identical result
126
127 mbs.AddLoad(Torque(markerNumber=mRigid, loadVector=[0,0,-0.1]))
128
129 #fix position of roll:
130 for i in range(0,2):
131 mRigidC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=i))
132 #mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mRigidC]))
133 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRigidC], stiffness = 100, damping=5, visualization= VObjectConnectorCoordinateSpringDamper(show = False)))
134
135
136 for i in range(len(cableList)):
137 #exu.Print("cable="+str(cableList[i]))
138 mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[i], numberOfSegments = nSegments))
139
140 nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments*(1+2*useFriction)))
141 mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mRigid, mCable], nodeNumber = nodeDataContactCable,
142 numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
143 frictionVelocityPenalty = 10, frictionCoefficient=2,
144 useSegmentNormals=False, #for this test
145 circleRadius = r2))
146
147#%%++++++++++++++++++++++++++++++++++
148#finally assemble and start computation
149mbs.Assemble()
150#exu.Print(mbs)
151
152simulationSettings = exu.SimulationSettings() #takes currently set values or default values
153
154fact = 300
155simulationSettings.timeIntegration.numberOfSteps = fact
156simulationSettings.timeIntegration.endTime = 0.0005*fact
157simulationSettings.solutionSettings.writeSolutionToFile = True
158simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
159#simulationSettings.solutionSettings.outputPrecision = 4
160#simulationSettings.displayComputationTime = True
161simulationSettings.timeIntegration.verboseMode = 1
162
163simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8 #10000
164simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
165
166#simulationSettings.timeIntegration.discontinuous.maxIterations = 5
167# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 0.001 #with
168
169simulationSettings.timeIntegration.newton.useModifiedNewton = False
170simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
171#simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
172#simulationSettings.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 = True
176simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
177simulationSettings.displayStatistics = True #just in this example ...
178
179SC.visualizationSettings.bodies.showNumbers = False
180SC.visualizationSettings.nodes.defaultSize = 0.01
181SC.visualizationSettings.markers.defaultSize = 0.01
182SC.visualizationSettings.connectors.defaultSize = 0.01
183SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
184SC.visualizationSettings.connectors.showContact = 1
185
186simulationSettings.solutionSettings.solutionInformation = "ANCF cable with rigid contact"
187
188# useGraphics=False
189if useGraphics:
190 exu.StartRenderer()
191 mbs.WaitForUserToContinue()
192
193solveDynamic = True
194if solveDynamic:
195
196 mbs.SolveDynamic(simulationSettings)
197
198 sol = mbs.systemData.GetODE2Coordinates()
199 u = sol[len(sol)-3]
200 exu.Print('tip displacement: y='+str(u))
201else:
202 simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-10*100 #can be quite small; WHY?
203 simulationSettings.staticSolver.verboseMode = 1
204 simulationSettings.staticSolver.numberOfLoadSteps = 20*2
205 simulationSettings.staticSolver.loadStepGeometric = False;
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.staticSolver.pauseAfterEachStep = False
215 simulationSettings.staticSolver.stabilizerODE2term = 50
216
217 mbs.SolveStatic(simulationSettings)
218
219 sol = mbs.systemData.GetODE2Coordinates()
220 n = len(sol)
221 u=sol[n-4]
222 exu.Print('static tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
223
224#put outside if
225exudynTestGlobals.testError = u - (-0.014187561328096003) #until 2022-03-09 (old ObjectContactFrictionCircleCable2D): -0.014188649931870346 #2019-12-26: -0.014188649931870346; 2019-12-16: (-0.01418281035370442);
226exudynTestGlobals.testResult = u
227exu.Print("test result=",exudynTestGlobals.testResult)
228
229if useGraphics:
230 SC.WaitForRenderEngineStopFlag()
231 exu.StopRenderer() #safely close rendering window!