ALEANCFpipe.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF ALE Cable2D 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#exu.SetOutputPrecision(16)
 19
 20#background
 21rect = [-2.5,-2,2.5,1] #xmin,ymin,xmax,ymax
 22background0 = {'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
 23background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 24oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0])))
 25
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27#cable:
 28mypi = 3.141592653589793
 29
 30L=2                     # length of ANCF element in m
 31#L=mypi                 # length of ANCF element in m
 32Em=2.07e11              # Young's modulus of ANCF element in N/m^2
 33rho=7800                # density of ANCF element in kg/m^3
 34b=0.1                   # width of rectangular ANCF element in m
 35h=0.1                   # height of rectangular ANCF element in m
 36A=b*h                   # cross sectional area of ANCF element in m^2
 37I=b*h**3/12             # second moment of area of ANCF element in m^4
 38EI = Em*I
 39rhoA = rho*A
 40EA = Em*A
 41movingMassFactor = 1
 42vALE = 2.3*1
 43
 44#f=3*E*I/L**2            # tip load applied to ANCF element in N
 45g=9.81
 46
 47#+++++++++++++++++++++++++++++++++++++++++++++++++
 48#paper pipe:
 49pipePaper=True
 50if pipePaper:
 51    L=1
 52    vALE = 10 #check sign (direction of fuild?)
 53    EI = 10 #*0.01
 54    rhoA=10                 #fluid+pipe
 55    EA = 100000*10 #*10        #not given in paper
 56    movingMassFactor = 1    #pipe has 8kg/m and fluid has 2kg/m
 57    g=0.1*9.81             #small perturbation
 58
 59print("L="+str(L))
 60print("EI="+str(EI))
 61print("EA="+str(EA))
 62print("rhoA="+str(rhoA))
 63
 64nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 65mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 66
 67cableList=[]        #for cable elements
 68nodeList=[]  #for nodes of cable
 69markerList=[]       #for nodes
 70
 71useALE = True
 72
 73
 74if useALE:
 75    nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], initialCoordinates=[0], initialCoordinates_t=[vALE]))
 76    mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity
 77    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mALE], offset=vALE, velocityLevel = True)) # for static computation
 78
 79nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 80nodeList+=[nc0]
 81nElements = 16
 82lElem = L / nElements
 83for i in range(nElements):
 84    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 85    #nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[L*2/3.1415926,L*2/3.1415926,0,1]))
 86    nodeList+=[nLast]
 87    if useALE:
 88        elem=mbs.AddObject(ALECable2D(physicsLength=lElem, physicsMassPerLength=rhoA,
 89                                      physicsBendingStiffness=EI, physicsAxialStiffness=EA, physicsMovingMassFactor=movingMassFactor,
 90                                      nodeNumbers=[nodeList[i],nodeList[i+1],nALE]))
 91    else:
 92        elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rhoA, physicsBendingStiffness=EI,
 93                                   physicsAxialStiffness=EA, nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 94
 95    cableList+=[elem]
 96    mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber = elem))
 97    mbs.AddLoad(Gravity(markerNumber=mBody, loadVector=[0,-g,0]))
 98
 99
100
101mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=0))
102mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=1))
103mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = int(nc0)+1*0, coordinate=3))
104
105mANCF3 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=0)) #tip constraint
106mANCF4 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLast, coordinate=1)) #tip constraint
107
108mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
109mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
110mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
111#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF3]))
112#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF4]))
113
114#add gravity:
115markerList=[]
116for i in range(len(nodeList)):
117    m = mbs.AddMarker(MarkerNodePosition(nodeNumber=nodeList[i]))
118    markerList+=[m]
119
120
121#a = 0.1     #y-dim/2 of gondula
122#b = 0.001    #x-dim/2 of gondula
123#massRigid = 12*0.01
124#inertiaRigid = massRigid/12*(2*a)**2
125#g = 9.81    # gravity
126#
127#slidingCoordinateInit = lElem*1.5 #0.75*L
128#initialLocalMarker = 1 #second element
129#if nElements<2:
130#    slidingCoordinateInit /= 3.
131#    initialLocalMarker = 0
132#
133
134mbs.Assemble()
135print(mbs)
136
137simulationSettings = exu.SimulationSettings() #takes currently set values or default values
138#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
139#simulationSettings.outputPrecision = 16
140
141fact = 20000
142simulationSettings.timeIntegration.numberOfSteps = 1*fact
143simulationSettings.timeIntegration.endTime = 0.001*fact
144simulationSettings.solutionSettings.writeSolutionToFile = True
145simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/2000
146#simulationSettings.solutionSettings.outputPrecision = 4
147simulationSettings.displayComputationTime = True
148simulationSettings.timeIntegration.verboseMode = 1
149
150simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8 #10000
151simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
152
153simulationSettings.timeIntegration.newton.useModifiedNewton = False
154simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
155simulationSettings.timeIntegration.newton.numericalDifferentiation.addReferenceCoordinatesToEpsilon = False
156simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1.e-3
157simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-8 #6.055454452393343e-06*0.0001 #eps^(1/3)
158simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1e8
159# simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
160# simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
161simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
162simulationSettings.pauseAfterEachStep = False
163simulationSettings.displayStatistics = True
164
165#SC.visualizationSettings.nodes.showNumbers = True
166SC.visualizationSettings.bodies.showNumbers = False
167#SC.visualizationSettings.connectors.showNumbers = True
168SC.visualizationSettings.nodes.defaultSize = 0.01
169SC.visualizationSettings.markers.defaultSize = 0.01
170SC.visualizationSettings.connectors.defaultSize = 0.01
171SC.visualizationSettings.contact.contactPointsDefaultSize = 0.005
172SC.visualizationSettings.connectors.showContact = 1
173
174simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
175
176solveDynamic = True
177if solveDynamic:
178    exu.StartRenderer()
179    #mbs.WaitForUserToContinue()
180
181    mbs.SolveDynamic(simulationSettings,
182                     solverType=exu.DynamicSolverType.TrapezoidalIndex2)
183
184    SC.WaitForRenderEngineStopFlag()
185    exu.StopRenderer() #safely close rendering window!
186
187else:
188    simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-8 #*100 #can be quite small; WHY?
189    simulationSettings.staticSolver.newton.numericalDifferentiation.doSystemWideDifferentiation = False
190    simulationSettings.staticSolver.verboseMode = 2
191    simulationSettings.staticSolver.numberOfLoadSteps  = 20#20*2
192    simulationSettings.staticSolver.loadStepGeometric = True;
193    simulationSettings.staticSolver.loadStepGeometricRange = 1e3;
194
195    simulationSettings.staticSolver.newton.relativeTolerance = 1e-5 #1e-5*100
196    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
197    simulationSettings.staticSolver.newton.maxIterations = 20 #50 for bending into circle
198
199    simulationSettings.staticSolver.discontinuous.iterationTolerance = 0.1
200    #simulationSettings.staticSolver.discontinuous.maxIterations = 5
201    simulationSettings.staticSolver.pauseAfterEachStep = False
202    simulationSettings.staticSolver.stabilizerODE2term = 100*0.0
203
204    exu.StartRenderer()
205
206    mbs.SolveStatic(simulationSettings)
207
208    sol = mbs.systemData.GetODE2Coordinates()
209    n = len(sol)
210    print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
211    sol_t = mbs.systemData.GetODE2Coordinates_t()
212    print('vALE='+str(sol_t[0]))
213
214    #print('sol='+str(sol))
215    print('sol_t='+str(sol_t))
216
217
218    SC.WaitForRenderEngineStopFlag()
219    exu.StopRenderer() #safely close rendering window!