ANCFcontactCircle.py

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