ANCFALEtest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF ALE with under gravity
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-02-17
  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.utilities import * #includes itemInterface and rigidBodyUtilities
 15import exudyn.graphics as graphics #only import if it does not conflict
 16
 17import numpy as np
 18from math import sqrt, sin, cos
 19
 20import matplotlib.pyplot as plt
 21import matplotlib.ticker as ticker
 22
 23#plt.clear('all')
 24#plt.rcParams['text.usetex'] = True #slows down figures
 25
 26#%%++++++++++++++++++++++++++++++++++++++++
 27useGraphics = True
 28plotResults=False
 29
 30tEnd = 2
 31h= 1e-3
 32
 33SC = exu.SystemContainer()
 34mbs = SC.AddSystem()
 35
 36#++++++++++++++++++++++++++++++++++
 37#initialize variables
 38
 39useGraphics = True
 40
 41nElements = 8
 42vALE0=1 #initial velocity
 43h= 2e-3
 44tEnd = 2 #fails at higher times ... check if this is just unstable due to very flexible beam
 45
 46damper = 0.01 #0.1: standard for parameter variation; 0.001: almost no damping, but solution is still oscillating at evaluation period
 47
 48
 49L=1.        #length of ANCF element in m
 50rhoA=10     #beam + discrete masses
 51
 52EA=1e5
 53EI=10
 54
 55movingMassFactor = 1 #factor for beam;1=axially moving beam, <1: pipe
 56
 57useCoordinateSpringDamper=True #use damping for every node use this for Yang Example
 58
 59# #additional bending and axial damping
 60bendingDamping=0 # for ALE Element
 61axialDamping=0 # for ALE Element
 62
 63#generate coordinate marker
 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
 67#++++++++++++++++++++++++++++++++++++++++
 68#create ALE node
 69#start rope moving upwards
 70nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0],
 71                                   initialCoordinates=[0], initialCoordinates_t=[vALE0]))
 72mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity
 73mbs.variables['nALE'] = nALE
 74
 75if useGraphics:
 76    mbs.variables['sALEpos'] = mbs.AddSensor(SensorNode(nodeNumber=nALE, fileName='solution/nodeALEpos.txt',
 77                             outputVariableType=exu.OutputVariableType.Coordinates))
 78    mbs.variables['sALEvel'] = mbs.AddSensor(SensorNode(nodeNumber=nALE, fileName='solution/nodeALEvel.txt',
 79                             outputVariableType=exu.OutputVariableType.Coordinates_t))
 80
 81oCCvALE=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mALE], offset=vALE0*0,  #for static computation
 82                                           velocityLevel = False,
 83                                           activeConnector = True,
 84                                           visualization=VCoordinateConstraint(show=False))) # False for static computation
 85
 86#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 87#create one beam template
 88cable = ALECable2D(#physicsLength=L,
 89                    physicsMassPerLength=rhoA,
 90                    physicsBendingStiffness=EI,
 91                    physicsAxialStiffness=EA,
 92                    physicsBendingDamping=bendingDamping,
 93                    physicsAxialDamping=axialDamping,
 94                    physicsMovingMassFactor=movingMassFactor,
 95                    nodeNumbers=[0,0,nALE],
 96                    # physicsUseCouplingTerms = True,
 97                    # useReducedOrderIntegration = True, #faster
 98                    )
 99
100phi = 0.25*pi/2
101#alternative to mbs.AddObject(ALECable2D(...)) with nodes:
102ancf=GenerateStraightLineANCFCable2D(mbs=mbs,
103                positionOfNode0=[0,0,0], positionOfNode1=[L*cos(phi),L*sin(phi),0],
104                numberOfElements=nElements,
105                cableTemplate=cable, #this defines the beam element properties
106                massProportionalLoad = [0,-9.81,0], #add larger gravity for larger deformation
107                # fixedConstraintsNode0 = [1,1,1,1], #fixed
108                fixedConstraintsNode0 = [1,1,1*0,1*0], #fixed
109                fixedConstraintsNode1 = [1,1,1*0,1*0]) #fixed
110
111ancfNodes = ancf[0]
112ancfObjects = ancf[1]
113for oCC in ancf[4]:
114    mbs.SetObjectParameter(oCC,'VdrawSize',0.005)
115
116
117if useCoordinateSpringDamper:
118    for node in ancfNodes:
119        mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = node, coordinate=0))
120        mbs.AddObject(CoordinateSpringDamper(markerNumbers = [mGround , mANCF0],
121                                             stiffness = 0, damping = 1*damper,
122                                             visualization=VCoordinateSpringDamper(show=False)))
123
124        mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = node, coordinate=1))
125        mbs.AddObject(CoordinateSpringDamper(markerNumbers = [mGround, mANCF0],
126                                             stiffness = 0, damping = damper,
127                                             visualization=VCoordinateSpringDamper(show=False)))
128
129#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
130midNode = ancfNodes[int(len(ancfNodes)/4)] #gives correct result for odd node numbers / even nElements
131sensorFileName = 'solution/beamALEmidPoint.txt'
132sMid = mbs.AddSensor(SensorNode(nodeNumber=midNode, fileName=sensorFileName,
133                            outputVariableType=exu.OutputVariableType.Displacement))
134
135
136mbs.Assemble()
137# print(mbs)
138#mbs.systemData.Info()
139
140simulationSettings = exu.SimulationSettings() #takes currently set values or default values
141if useGraphics:
142    verboseMode = 1
143else:
144    verboseMode = 0
145
146
147simulationSettings.solutionSettings.writeSolutionToFile = False
148simulationSettings.solutionSettings.sensorsWritePeriod = h
149#simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6 #10000
150simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10 #default:1e-10
151simulationSettings.timeIntegration.verboseMode = verboseMode
152simulationSettings.staticSolver.verboseMode = verboseMode
153
154simulationSettings.timeIntegration.newton.useModifiedNewton = True
155# simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
156simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
157simulationSettings.timeIntegration.adaptiveStep = True #disable adaptive step reduction
158
159simulationSettings.displayStatistics = True
160SC.visualizationSettings.loads.show = False
161
162if useGraphics:
163    exu.StartRenderer()
164    mbs.WaitForUserToContinue()
165
166#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
167#static step
168simulationSettings.staticSolver.numberOfLoadSteps=10
169
170success = mbs.SolveStatic(simulationSettings, updateInitialValues=True)
171
172
173#turn on moving beam:
174mbs.SetObjectParameter(oCCvALE, 'activeConnector', True)
175mbs.SetObjectParameter(oCCvALE, 'velocityLevel', True)
176mbs.SetObjectParameter(oCCvALE, 'offset', vALE0)
177
178#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
179#turn on vALE velocity (could also be done in modifying coordinates):
180#rope decelerates due to gravity and then runs backwards
181simulationSettings.timeIntegration.numberOfSteps = int(1/h)
182simulationSettings.timeIntegration.endTime = 1
183success = mbs.SolveDynamic(simulationSettings,
184                            exudyn.DynamicSolverType.TrapezoidalIndex2,
185                            updateInitialValues=True)
186mbs.systemData.SetODE2Coordinates_tt(coordinates = mbs.systemData.GetODE2Coordinates_tt(),
187                                     configuration = exudyn.ConfigurationType.Initial)
188
189if useGraphics:
190    mbs.WaitForUserToContinue()
191
192#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
193#finally: solve dynamic problem under self weight
194mbs.SetObjectParameter(oCCvALE, 'activeConnector', False) #rope under self-weight
195mbs.SetObjectParameter(oCCvALE, 'velocityLevel', False)
196mbs.SetObjectParameter(oCCvALE, 'offset', 0)
197
198simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
199simulationSettings.timeIntegration.startTime = 1
200simulationSettings.solutionSettings.appendToFile = True #continue solution
201simulationSettings.timeIntegration.endTime = tEnd
202
203success = mbs.SolveDynamic(simulationSettings,
204                           exudyn.DynamicSolverType.TrapezoidalIndex2
205                           )
206
207if useGraphics:
208    SC.WaitForRenderEngineStopFlag()
209    exu.StopRenderer() #safely close rendering window!
210
211    plt.close('all')
212    if True:
213
214        plt.figure("ALE pos/vel")
215        mbs.PlotSensor(sensorNumbers=[mbs.variables['sALEpos'],mbs.variables['sALEvel']], components=[0,0])
216
217    plt.figure("midpoint")
218    data0 = np.loadtxt('solution/beamALEmidPoint.txt', comments='#', delimiter=',')
219    y0 = data0[0,2]
220    plt.plot(data0[:,0],data0[:,2]-y0,'b-',label='midPointDeflection')
221    ax=plt.gca()
222    ax.grid(True,'major','both')
223    plt.tight_layout()
224    plt.legend()
225    plt.show()