ANCFcontactCircleTest.py

You can view and download this file on Github: ANCFcontactCircleTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF cable element contact with circle; test model fo ObjectContactCircleCable2D
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-09-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
 16useGraphics = True #without test
 17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 18#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 19try: #only if called from test suite
 20    from modelUnitTests import exudynTestGlobals #for globally storing test results
 21    useGraphics = exudynTestGlobals.useGraphics
 22except:
 23    class ExudynTestGlobals:
 24        pass
 25    exudynTestGlobals = ExudynTestGlobals()
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27
 28SC = exu.SystemContainer()
 29mbs = SC.AddSystem()
 30
 31#background
 32rect = [-2,-2,4,2] #xmin,ymin,xmax,ymax
 33background0 = {'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
 34background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 35oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0, background1])))
 36
 37#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 38#cable:
 39mypi = 3.141592653589793
 40
 41L=2                     # length of ANCF element in m
 42#L=mypi                 # length of ANCF element in m
 43E=2.07e11               # Young's modulus of ANCF element in N/m^2
 44rho=7800                # density of ANCF element in kg/m^3
 45b=0.001                 # width of rectangular ANCF element in m
 46h=0.001                 # height of rectangular ANCF element in m
 47A=b*h                   # cross sectional area of ANCF element in m^2
 48I=b*h**3/12             # second moment of area of ANCF element in m^4
 49f=3*E*I/L**2            # tip load applied to ANCF element in N
 50
 51exu.Print("load f="+str(f))
 52exu.Print("EI="+str(E*I))
 53
 54nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 55mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 56
 57cableList=[]        #for cable elements
 58nodeList=[]  #for nodes of cable
 59markerList=[]       #for nodes
 60nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 61nodeList+=[nc0]
 62nElements = 8 #8 original in test
 63lElem = L / nElements
 64for i in range(nElements):
 65    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 66    nodeList+=[nLast]
 67    elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
 68                               physicsBendingStiffness=E*I, physicsAxialStiffness=E*A,
 69                               nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 70    cableList+=[elem]
 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, -40*2*rho*A*fact*lElem, 0])) #will be changed in load steps
 88
 89#mANCFend = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nodeList[-1], coordinate=1)) #last marker
 90#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCFend]))
 91
 92#mGroundTip = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[L,0,0]))
 93#mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundTip,markerList[-1]], stiffness=[10,10,10], damping=[0.1,0.1,0.1]))
 94
 95#mGroundTip2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[L,0.2,0]))
 96#mbs.AddObject(SpringDamper(markerNumbers=[mGroundTip2,markerList[-1]], stiffness=0.1, referenceLength=0.2))
 97
 98#mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
 99#mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -1e8, 0])) #will be changed in load steps
100
101#mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
102#mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*1*mypi]))
103
104#mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
105#mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 3*E*I*1*mypi]))
106
107cStiffness = 1e3
108cDamping = 0.02*cStiffness
109useContact = False
110if useContact:
111    tipContact = False
112    if tipContact:
113        nodeData = mbs.AddNode(NodeGenericData(initialCoordinates=[0],numberOfDataCoordinates=1))
114        mbs.AddObject(ObjectContactCoordinate(markerNumbers=[mGround, mANCFend],nodeNumber = nodeData, contactStiffness = cStiffness, contactDamping=0*cDamping, offset = -0.8))
115    else:
116        for i in range(len(nodeList)):
117            mNC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nodeList[i], coordinate=1))
118            nodeData = mbs.AddNode(NodeGenericData(initialCoordinates=[1],numberOfDataCoordinates=1)) #start with gap!
119            mbs.AddObject(ObjectContactCoordinate(markerNumbers=[mGround, mNC], nodeNumber = nodeData, contactStiffness = cStiffness, contactDamping=0*cDamping, offset = -1))
120
121nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element
122initialGapList = [0.1]*nSegments #initial gap of 0.1
123
124mGroundCircle = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.75*L,-0.5,0]))
125mGroundCircle2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.25*L,-0.15,0]))
126
127#mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=elem, numberOfSegments = nSegments))
128#nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
129#mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle, mCable], nodeNumber = nodeDataContactCable,
130#                                         numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=cDamping,
131#                                         circleRadius = 0.4, offset = 0))
132for i in range(len(cableList)):
133    mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[i], numberOfSegments = nSegments))
134    nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
135    mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle, mCable], nodeNumber = nodeDataContactCable,
136                                             numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=0*cDamping,
137                                             circleRadius = 0.2, offset = 0))
138    nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments))
139    mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mGroundCircle2, mCable], nodeNumber = nodeDataContactCable,
140                                             numberOfContactSegments=nSegments, contactStiffness = cStiffness, contactDamping=0*cDamping,
141                                             circleRadius = 0.1, offset = 0))
142
143
144#mbs.systemData.Info()
145
146mbs.Assemble()
147#exu.Print(mbs)
148
149simulationSettings = exu.SimulationSettings() #takes currently set values or default values
150
151simulationSettings.solutionSettings.writeSolutionToFile = True
152#simulationSettings.solutionSettings.outputPrecision = 4
153simulationSettings.displayComputationTime = False
154
155simulationSettings.displayStatistics = False
156
157#SC.visualizationSettings.nodes.showNumbers = True
158SC.visualizationSettings.bodies.showNumbers = False
159#SC.visualizationSettings.connectors.showNumbers = True
160SC.visualizationSettings.nodes.defaultSize = 0.01
161SC.visualizationSettings.markers.defaultSize = 0.01
162SC.visualizationSettings.connectors.defaultSize = 0.01
163SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
164SC.visualizationSettings.connectors.showContact = 1
165
166simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
167
168simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-10 #can be quite small; WHY?
169simulationSettings.staticSolver.verboseMode = 0 #otherwise, load steps are shown ...
170simulationSettings.staticSolver.numberOfLoadSteps  = 40
171simulationSettings.staticSolver.loadStepGeometric = True;
172simulationSettings.staticSolver.loadStepGeometricRange = 1e4;
173simulationSettings.staticSolver.adaptiveStep = False
174
175simulationSettings.staticSolver.newton.relativeTolerance = 1e-7 #10000
176simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
177simulationSettings.staticSolver.newton.maxIterations = 30 #50 for bending into circle
178
179simulationSettings.staticSolver.discontinuous.iterationTolerance = 1
180simulationSettings.staticSolver.stabilizerODE2term = 2 #may only act on position degrees of freedom
181
182if useGraphics:
183    simulationSettings.staticSolver.verboseMode = 1 #otherwise, load steps are shown ...
184    simulationSettings.staticSolver.verboseModeFile = 0 #otherwise, load steps are shown ...
185    simulationSettings.displayStatistics = True
186
187    exu.StartRenderer()
188
189#mbs.WaitForUserToContinue()
190mbs.SolveStatic(simulationSettings) #183 Newton iterations, 0.114 seconds
191
192sol = mbs.systemData.GetODE2Coordinates()
193n = len(sol)
194exu.Print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
195
196if useGraphics:
197    SC.WaitForRenderEngineStopFlag()
198    exu.StopRenderer() #safely close rendering window!
199
200exudynTestGlobals.testError = sol[n-3] - (-0.4842698420787613) #-0.4842698420787613 ; 2021-05-07 (deactivated StaticSolveOldSolver):-0.4842656133238705  #2019-12-17(relTol=1e-7 / up to 7 digits accurate): -0.4842656547442095;  2019-11-22: (-0.4844812763485709) (with relTol=1e-5);  y-displacement
201exudynTestGlobals.testResult = sol[n-3]