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