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!