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