rigidRotor3DFWBW.py
You can view and download this file on Github: rigidRotor3DFWBW.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example with 3D rotor; rotates around X-axis
5# showing forward (FW) and backward (BW) whirls and FW/BW whirl resonances
6#
7# Author: Johannes Gerstmayr
8# Date: 2019-12-05
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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13import sys
14sys.path.append('../TestModels') #for modelUnitTest as this example may be used also as a unit test
15
16import exudyn as exu
17from exudyn.itemInterface import *
18from exudyn.utilities import *
19from math import cos, sin
20
21import time
22import numpy as np
23
24SC = exu.SystemContainer()
25mbs = SC.AddSystem()
26print('EXUDYN version='+exu.GetVersionString())
27
28L=1 #rotor axis length
29isSymmetric = True
30symStr = "General"
31if isSymmetric:
32 L0 = 0.5 #0.5 (symmetric rotor); position of rotor on x-axis
33 symStr = "Laval"
34 s = "symmetric"
35else :
36 L0 = 0.9 #default: 0.9m; position of rotor on x-axis
37 s = "unsymmetric"
38
39L1 = L-L0 #
40m = 2 #mass in kg
41r = 0.5*1.5 #radius for disc mass distribution
42lRotor = 0.2 #length of rotor disk
43k = 800 #stiffness of (all/both) springs in rotor in N/m
44Jxx = 0.5*m*r**2 #polar moment of inertia
45Jyyzz = 0.25*m*r**2 + 1/12.*m*lRotor**2 #moment of inertia for y and z axes
46
47omega0=np.sqrt(2*k/m) #linear system
48
49D0 = 0.002 #dimensionless damping
50d = 2*omega0*D0*m #damping constant in N/(m/s)
51
52f0 = 0*omega0/(2*np.pi) #frequency start (Hz)
53f1 = 2.*omega0/(2*np.pi) #frequency end (Hz)
54
55torque = 0.2*2 #driving torque; Nm ; 0.1Nm does not surpass critical speed; 0.2Nm works
56eps = 2e-3*0. #0.74; excentricity of mass in y-direction
57 #symmetric rotor: 2e-3 gives large oscillations;
58 #symmetric rotor: 0.74*2e-3 shows kink in runup curve
59
60omegaInitial = 0*4*omega0 #initial rotation speed in rad/s
61
62tEnd = 20 #end time of simulation
63steps = int(tEnd*1000) #number of steps
64
65modeStr = "BW" #"FW" or "BW"
66sign = 1
67if modeStr == "BW": sign = -1
68
69fRes = omega0/(2*np.pi)
70print('symmetric rotor resonance frequency (Hz)= '+str(fRes))
71#print('runup over '+str(tEnd)+' seconds, fStart='+str(f0)+'Hz, fEnd='+str(f1)+'Hz')
72
73#3D load vector
74loadFact = 1
75n1 = 0 #rotor node number, will be defined later
76def UFload(mbs, t, load):
77 #compute angle of rotation (only true if rotation about 1 axis)
78 phi = mbs.GetNodeOutput(nodeNumber=n1, variableType=exu.OutputVariableType.Rotation)[0]
79 return [0, sign*loadFact*cos(phi), loadFact*sin(phi)]
80
81#background1 = GraphicsDataOrthoCube(0,0,0,.5,0.5,0.5,[0.3,0.3,0.9,1])
82
83#draw RGB-frame at origin
84p=[0,0,0]
85lFrame = 0.9
86tFrame = 0.01
87backgroundX = GraphicsDataCylinder(p,[lFrame,0,0],tFrame,[0.9,0.3,0.3,1],12)
88backgroundY = GraphicsDataCylinder(p,[0,lFrame,0],tFrame,[0.3,0.9,0.3,1],12)
89backgroundZ = GraphicsDataCylinder(p,[0,0,lFrame],tFrame,[0.3,0.3,0.9,1],12)
90
91#visualize forces:
92lVec = r*0.8
93f1=0.08 #arrow size
94f2=0.04 #arrow size 2
95rigid = 0 #will be set later!
96#ground user function, but draws load on body
97def UFgraphics(mainSystem, bodyNr):
98 t=mainSystem.systemData.GetTime()
99 #get rotation for rigid body:
100 rot = mbs.GetObjectOutputBody(rigid, exu.OutputVariableType.Rotation,
101 [0,0,0], exu.ConfigurationType.Visualization)
102 omega = mbs.GetObjectOutputBody(rigid, exu.OutputVariableType.AngularVelocity,
103 [0,0,0], exu.ConfigurationType.Visualization)
104 phi = rot[0]
105 phi += 0.9*0.05*omega[0] #graphics is always delayed ...
106 #print("rot=", phi)
107
108 #g0=GraphicsDataSphere(point=p0,radius=0.01, color=color4green)
109 #g1=GraphicsDataSphere(point=p1,radius=0.01, color=color4red)
110 xVec = 0.5+lRotor
111 p0 = [xVec,0,0]
112 p1 = [xVec,sign*lVec*cos(phi), lVec*sin(phi)]
113 p2 = [xVec,sign*(1-f1)*lVec*cos(phi)-f2*lVec*sin(phi), (1-f1)*lVec*sin(phi)+f2*sign*lVec*cos(phi)]
114 p3 = [xVec,sign*(1-f1)*lVec*cos(phi)+f2*lVec*sin(phi), (1-f1)*lVec*sin(phi)-f2*sign*lVec*cos(phi)]
115 gLine = {'type':'Line', 'data': p0+p1+p2+p1+p3, 'color':color4red}
116 return [gLine]#,g0,g1] #return list of graphics data
117
118mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
119 visualization=VObjectGround(graphicsDataUserFunction=UFgraphics)))
120
121#rotor is rotating around x-axis
122ep0 = eulerParameters0 #no rotation
123ep_t0 = AngularVelocity2EulerParameters_t([omegaInitial,0,0], ep0)
124print(ep_t0)
125
126p0 = [L0-0.5*L,eps,0] #reference position
127v0 = [0.,0.,0.] #initial translational velocity
128
129#node for Rigid2D body: px, py, phi:
130n1=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates = p0+ep0, initialVelocities=v0+list(ep_t0)))
131
132#ground nodes
133nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [-L/2,0,0]))
134nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [ L/2,0,0]))
135
136#add mass point (this is a 3D object with 3 coordinates):
137gRotor = GraphicsDataCylinder([-lRotor*0.5,0,0],[lRotor,0,0],r,[0.3,0.3,0.9,1],128)
138gRotor2 = GraphicsDataCylinder([-L0,0,0],[L,0,0],r*0.05,[0.3,0.3,0.9,1],16)
139gRotor3 = [backgroundX, backgroundY, backgroundZ]
140rigid = mbs.AddObject(RigidBody(physicsMass=m, physicsInertia=[Jxx,Jyyzz,Jyyzz,0,0,0], nodeNumber = n1, visualization=VObjectRigidBody2D(graphicsData=[gRotor, gRotor2]+gRotor3)))
141
142mbs.AddSensor(SensorBody(bodyNumber=rigid,
143 fileName='solution/runupDisplacement'+modeStr+'.txt',
144 outputVariableType=exu.OutputVariableType.Displacement))
145mbs.AddSensor(SensorBody(bodyNumber=rigid,
146 fileName='solution/runupAngularVelocity'+modeStr+'.txt',
147 outputVariableType=exu.OutputVariableType.AngularVelocity))
148
149#marker for ground (=fixed):
150groundMarker0=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround0))
151groundMarker1=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround1))
152
153#marker for rotor axis and support:
154rotorAxisMarker0 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[-L0,-eps,0]))
155rotorAxisMarker1 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[ L1,-eps,0]))
156
157
158#++++++++++++++++++++++++++++++++++++
159mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker0, rotorAxisMarker0],
160 stiffness=[k,k,k], damping=[d, d, d]))
161mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker1, rotorAxisMarker1],
162 stiffness=[0,k,k], damping=[0, d, d])) #do not constrain x-axis twice
163
164#coordinate markers for loads:
165rotorMarkerUy=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate=1))
166rotorMarkerUz=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate=2))
167
168#add torque:
169rotorRigidMarker =mbs.AddMarker(MarkerBodyRigid(bodyNumber=rigid, localPosition=[0,0,0]))
170mbs.AddLoad(Torque(markerNumber=rotorRigidMarker,
171 loadVector=[torque,0,0],
172 visualization={'show':False}))
173
174#BW/FW load vector:
175mbs.AddLoad(Force(markerNumber=rotorRigidMarker,
176 loadVector=[0,0,0],
177 loadVectorUserFunction=UFload))
178
179#print(mbs)
180mbs.Assemble()
181#mbs.systemData.Info()
182
183simulationSettings = exu.SimulationSettings()
184simulationSettings.solutionSettings.solutionWritePeriod = 1e-5 #output interval
185simulationSettings.solutionSettings.sensorsWritePeriod = 1e-5 #output interval
186
187simulationSettings.solutionSettings.solutionInformation = "Runup of "+s+" rotor: "+modeStr + ", " + symStr
188
189simulationSettings.timeIntegration.numberOfSteps = steps
190simulationSettings.timeIntegration.endTime = tEnd
191simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
192simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
193
194simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
195
196#create animations (causes slow simulation):
197createAnimation=True
198if createAnimation:
199 simulationSettings.solutionSettings.recordImagesInterval = 0.05
200 SC.visualizationSettings.exportImages.saveImageFileName = "images/frame"
201 SC.visualizationSettings.window.renderWindowSize = [1600,1080]
202 SC.visualizationSettings.openGL.multiSampling = 4
203
204#visualize loads:
205# SC.visualizationSettings.loads.fixedLoadSize=False
206# SC.visualizationSettings.loads.loadSizeFactor=100
207# SC.visualizationSettings.loads.drawSimplified = False
208
209if True:
210 exu.StartRenderer() #start graphics visualization
211 mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
212
213 #start solver:
214 mbs.SolveDynamic(simulationSettings)
215
216 #SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
217 exu.StopRenderer() #safely close rendering window!
218
219 #evaluate final (=current) output values
220 u = mbs.GetNodeOutput(n1, exu.OutputVariableType.AngularVelocity)
221 print('omega final (Hz)=',(1./(2.*np.pi))*u)
222 #print('displacement=',u[0])
223
224
225##+++++++++++++++++++++++++++++++++++++++++++++++++++++
226import matplotlib.pyplot as plt
227import matplotlib.ticker as ticker
228
229if False:
230 plt.close('all') #close all plots
231
232
233 dataDispFW = np.loadtxt('solution/runupDisplacementFW.txt', comments='#', delimiter=',')
234 dataOmegaFW = np.loadtxt('solution/runupAngularVelocityFW.txt', comments='#', delimiter=',')
235 dataDispBW = np.loadtxt('solution/runupDisplacementBW.txt', comments='#', delimiter=',')
236 dataOmegaBW = np.loadtxt('solution/runupAngularVelocityBW.txt', comments='#', delimiter=',')
237 plt.rcParams.update({'font.size': 12})
238
239 if False:
240 plt.plot(dataDispFW[:,0], dataDispFW[:,3], 'r-', label='FW') #numerical solution
241 plt.plot(dataDispBW[:,0], dataDispBW[:,3], 'b-', label='BW') #numerical solution
242 plt.xlabel("time (s)")
243 plt.ylabel("z-displacement (m)")
244 ax=plt.gca() # get current axes
245 ax.grid(True, 'major', 'both')
246 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
247 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
248 plt.legend()
249 plt.tight_layout()
250
251 plt.figure()
252 plt.plot((1/(2*np.pi))*dataOmegaFW[:,1], dataDispFW[:,3], 'r-', label='FW', linewidth=0.7) #numerical solution
253 plt.plot((1/(2*np.pi))*dataOmegaBW[:,1], dataDispBW[:,3], 'b-', label='BW', linewidth=0.7) #numerical solution
254 plt.xlabel("angular velocity (1/s)")
255 plt.ylabel("z-displacement (m)")
256 ax=plt.gca() # get current axes
257 ax.grid(True, 'major', 'both')
258 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
259 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
260 plt.legend()
261 plt.tight_layout()
262
263 #save figure:
264 #plt.savefig("rotorBWFWrunup"+symStr+".pdf")
265
266 # plt.figure()
267 # plt.plot(dataOmegaFW[:,0], (1/(2*np.pi))*dataOmegaFW[:,1], 'r-') #numerical solution
268 # plt.xlabel("time (s)")
269 # plt.ylabel("angular velocity (1/s)")
270 # ax=plt.gca() # get current axes
271 # ax.grid(True, 'major', 'both')
272 # ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
273 # ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
274 # plt.tight_layout()
275
276 # plt.figure()
277 # plt.plot(dataDispFW[:,2], dataDispFW[:,3], 'r-') #numerical solution
278 # plt.xlabel("y-displacement (m)")
279 # plt.ylabel("z-displacement (m)")
280 # ax=plt.gca() # get current axes
281 # ax.grid(True, 'major', 'both')
282 # ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
283 # ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
284 # plt.tight_layout()
285
286 #plt.plot(data[n-500:n-1,1], data[n-500:n-1,2], 'r-') #numerical solution
287
288 plt.show()