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