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