rigidRotor3DbasicBehaviour.py
You can view and download this file on Github: rigidRotor3DbasicBehaviour.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example with 3D rotor, showing basic behaviour of rotor
5# show COM, unbalance for low, critical and high rotation speeds
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 *
19
20import time
21import numpy as np
22
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25print('EXUDYN version='+exu.GetVersionString())
26
27L=1 #rotor axis length
28isSymmetric = True
29if isSymmetric:
30 L0 = 0.5 #0.5 (symmetric rotor); position of rotor on x-axis
31else :
32 L0 = 0.9 #default: 0.9m; position of rotor on x-axis
33L1 = L-L0 #
34m = 2 #mass in kg
35r = 0.5*1.5 #radius for disc mass distribution
36lRotor = 0.2 #length of rotor disk
37k = 800 #stiffness of (all/both) springs in rotor in N/m
38Jxx = 0.5*m*r**2 #polar moment of inertia
39Jyyzz = 0.25*m*r**2 + 1/12.*m*lRotor**2 #moment of inertia for y and z axes
40
41omega0=np.sqrt(2*k/m) #linear system
42
43D0 = 0.002 #dimensionless damping
44d = 2*omega0*D0*m #damping constant in N/(m/s)
45
46f0 = 0*omega0/(2*np.pi) #frequency start (Hz)
47f1 = 2.*omega0/(2*np.pi) #frequency end (Hz)
48
49torque = 0*0.2 #driving torque; Nm ; 0.1Nm does not surpass critical speed; 0.2Nm works
50eps = 10e-3 # excentricity of mass in y-direction
51 #symmetric rotor: 2e-3 gives large oscillations;
52 #symmetric rotor: 0.74*2e-3 shows kink in runup curve
53#k*=1000
54
55modeStr=['slow (omega0/2)',
56 'critical (omega0)',
57 'fast (2*omega0)' ]
58mode = 2
59
60#add constraint on euler parameters or euler angles
61#add three cases
62
63if mode == 0:
64 omegaInitial = 0.5*omega0 #initial rotation speed in rad/s
65elif mode == 1:
66 omegaInitial = 1*omega0 #initial rotation speed in rad/s
67 eps *= 0.1
68 d *= 10
69elif mode == 2:
70 omegaInitial = 2*omega0 #initial rotation speed in rad/s
71
72tEnd = 50 #end time of simulation
73steps = 50000 #number of steps
74
75fRes = omega0/(2*np.pi)
76print('symmetric rotor resonance frequency (Hz)= '+str(fRes))
77print('omega intial (Hz)= '+str(omegaInitial/(2*np.pi)))
78#print('runup over '+str(tEnd)+' seconds, fStart='+str(f0)+'Hz, fEnd='+str(f1)+'Hz')
79
80
81# #user function for load
82# def userLoad(t, load):
83# #return load*np.sin(0.5*omega0*t) #gives resonance
84# if t>40: time.sleep(0.02) #make simulation slower
85# return load*Sweep(t, tEnd, f0, f1)
86# #return load*Sweep(t, tEnd, f1, f0) #backward sweep
87
88# #backward whirl excitation:
89# amp = 0.10 #in resonance: *0.01
90# def userLoadBWy(t, load):
91# return load*SweepCos(t, tEnd, f0, f1) #negative sign: BW, positive sign: FW
92# def userLoadBWz(t, load):
93# return load*Sweep(t, tEnd, f0, f1)
94#def userLoadBWx(t, load):
95# return load*np.sin(omegaInitial*t)
96#def userLoadBWy(t, load):
97# return -load*np.cos(omegaInitial*t) #negative sign: FW, positive sign: BW
98
99#background1 = GraphicsDataOrthoCube(0,0,0,.5,0.5,0.5,[0.3,0.3,0.9,1])
100
101#draw RGB-frame at origin
102p=[0,0,0]
103rDraw = 0.05*r
104lFrame = rDraw*1.2
105tFrame = 0.01*0.15
106backgroundX = GraphicsDataCylinder(p,[lFrame,0,0],tFrame,[0.9,0.3,0.3,1],12)
107backgroundY = GraphicsDataCylinder(p,[0,lFrame,0],tFrame*0.5,[0.3,0.9,0.3,1],12)
108backgroundZ = GraphicsDataCylinder(p,[0,0,lFrame],tFrame*0.5,[0.3,0.3,0.9,1],12)
109black=[0,0,0,1]
110textCOM = {'type':'Text', 'text': 'COM', 'color': black, 'position': [lFrame*1.1,0,0]}
111textSHAFT = {'type':'Text', 'text': 'SHAFT', 'color': black, 'position': [L-L0+0.1,-eps,0]}
112textY = {'type':'Text', 'text': 'Y', 'color': black, 'position': [0,lFrame*1.05,0]}
113textZ = {'type':'Text', 'text': 'Z', 'color': black, 'position': [0,0,lFrame*1.05]}
114
115#rotor is rotating around x-axis
116ep0 = eulerParameters0 #no rotation
117ep_t0 = AngularVelocity2EulerParameters_t([omegaInitial,0,0], ep0)
118print(ep_t0)
119
120p0 = [L0-0.5*L,eps,0] #reference position, displaced by eccentricity eps !
121v0 = [0.,0.,0.] #initial translational velocity
122
123#node for Rigid2D body: px, py, phi:
124n1=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates = p0+ep0,
125 initialVelocities=v0+list(ep_t0)))
126
127#ground nodes
128nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [-L/2,0,0]))
129nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [ L/2,0,0]))
130
131#add mass point (this is a 3D object with 3 coordinates):
132gRotor = GraphicsDataCylinder([-lRotor*0.2,0,0],[lRotor*0.4,0,0],rDraw,
133 [0.3,0.3,0.9,1],128)
134gRotor2 = GraphicsDataCylinder([-L0,-eps,0],[L,0,0],r*0.01*0.25,[0.6,0.6,0.6,1],16)
135gRotorCOM = GraphicsDataCylinder([-lRotor*0.1,0,0],[lRotor*0.6*0.1,0,0],r*0.01*0.5,
136 [0.3,0.9,0.3,1],16)
137gRotor3 = [backgroundX, backgroundY, backgroundZ, textCOM, textY, textZ, textSHAFT]
138rigid = mbs.AddObject(RigidBody(physicsMass=m,
139 physicsInertia=[Jxx,Jyyzz,Jyyzz,0,0,0],
140 nodeNumber = n1,
141 visualization=VObjectRigidBody2D(graphicsData=[gRotor, gRotor2, gRotorCOM]+gRotor3)))
142
143mbs.AddSensor(SensorBody(bodyNumber=rigid,
144 fileName='solution/rotorDisplacement.txt',
145 localPosition=[0,-eps,0],
146 outputVariableType=exu.OutputVariableType.Displacement))
147# mbs.AddSensor(SensorBody(bodyNumber=rigid,
148# fileName='solution/rotorAngularVelocity.txt',
149# outputVariableType=exu.OutputVariableType.AngularVelocity))
150
151#marker for ground (=fixed):
152groundMarker0=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround0))
153groundMarker1=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround1))
154
155#marker for rotor axis and support:
156rotorAxisMarker0 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[-L0,-eps,0]))
157rotorAxisMarker1 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[ L1,-eps,0]))
158
159
160#++++++++++++++++++++++++++++++++++++
161mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker0, rotorAxisMarker0],
162 stiffness=[k,k,k], damping=[d, d, d],
163 visualization=VCartesianSpringDamper(drawSize=0.002)))
164mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker1, rotorAxisMarker1],
165 stiffness=[0,k,k], damping=[0, d, d],
166 visualization=VCartesianSpringDamper(drawSize=0.002))) #do not constrain x-axis twice
167
168
169#add torque:
170# rotorRigidMarker =mbs.AddMarker(MarkerBodyRigid(bodyNumber=rigid, localPosition=[0,0,0]))
171# mbs.AddLoad(Torque(markerNumber=rotorRigidMarker, loadVector=[torque,0,0]))
172
173#constant velocity constraint:
174constantRotorVelocity = True
175if constantRotorVelocity :
176 mRotationAxis = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber = n1, rotationCoordinate=0))
177 mGroundCoordinate =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround0, coordinate=0))
178 mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundCoordinate, mRotationAxis],
179 offset=omegaInitial, velocityLevel=True,
180 visualization=VCoordinateConstraint(show=False))) #gives equation omegaMarker1 = offset
181
182
183#print(mbs)
184mbs.Assemble()
185#mbs.systemData.Info()
186
187simulationSettings = exu.SimulationSettings()
188simulationSettings.solutionSettings.solutionWritePeriod = 1e-5 #output interval
189simulationSettings.solutionSettings.sensorsWritePeriod = 1e-5 #output interval
190
191descrStr = "Laval rotor, resonance="+str(round(fRes,3))+", "+modeStr[mode]
192simulationSettings.solutionSettings.solutionInformation = descrStr
193
194simulationSettings.timeIntegration.numberOfSteps = steps
195simulationSettings.timeIntegration.endTime = tEnd
196simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
197simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
198
199simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
200SC.visualizationSettings.window.renderWindowSize = [1600,1080]
201SC.visualizationSettings.general.textSize = 22
202
203exu.StartRenderer() #start graphics visualization
204mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
205
206#simulate some time to get steady-state solution:
207mbs.SolveDynamic(simulationSettings)
208state = mbs.systemData.GetSystemState()
209
210#now simulate the steady state solution and record
211simulationSettings.timeIntegration.numberOfSteps = 10000
212simulationSettings.timeIntegration.endTime = 2.5
213
214#create animations (causes slow simulation):
215createAnimation=True
216if createAnimation:
217 mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
218 simulationSettings.solutionSettings.recordImagesInterval = 0.01
219 if mode == 1:
220 simulationSettings.timeIntegration.endTime = 1
221 simulationSettings.solutionSettings.recordImagesInterval = 0.0025
222 if mode == 2:
223 simulationSettings.timeIntegration.endTime = 0.5
224 simulationSettings.solutionSettings.recordImagesInterval = 0.001
225
226 SC.visualizationSettings.exportImages.saveImageFileName = "images/frame"
227
228 mbs.systemData.SetSystemState(state, configuration=exu.ConfigurationType.Initial)
229 mbs.SolveDynamic(simulationSettings)
230
231#SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
232exu.StopRenderer() #safely close rendering window!
233
234#evaluate final (=current) output values
235u = mbs.GetNodeOutput(n1, exu.OutputVariableType.AngularVelocity)
236print('omega final (Hz)=',u/(2*np.pi))
237#print('displacement=',u[0])
238c = mbs.GetNodeOutput(n1, exu.OutputVariableType.Coordinates)
239c_t = mbs.GetNodeOutput(n1, exu.OutputVariableType.Coordinates_t)
240print("nc=",c)
241print("nc_t=",c_t)
242
243##+++++++++++++++++++++++++++++++++++++++++++++++++++++
244import matplotlib.pyplot as plt
245import matplotlib.ticker as ticker
246
247if True:
248 plt.close('all') #close all plots
249
250 dataDisp = np.loadtxt('solution/rotorDisplacement.txt', comments='#', delimiter=',')
251
252 plt.plot(dataDisp[:,0], dataDisp[:,3], 'b-') #numerical solution
253 plt.xlabel("time (s)")
254 plt.ylabel("z-displacement (m)")
255
256 plt.figure()
257 plt.plot(dataDisp[:,2], dataDisp[:,3], 'r-') #numerical solution
258 plt.xlabel("y-displacement (m)")
259 plt.ylabel("z-displacement (m)")
260
261 #plt.plot(data[n-500:n-1,1], data[n-500:n-1,2], 'r-') #numerical solution
262
263 ax=plt.gca() # get current axes
264 ax.grid(True, 'major', 'both')
265 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
266 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
267 plt.tight_layout()
268 plt.show()