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