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()