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