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