rigidRotor3Drunup.py

You can view and download this file on Github: rigidRotor3Drunup.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example with 3D rotor, showing runup
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-12-05
  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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12import sys
 13sys.path.append('../TestModels')            #for modelUnitTest as this example may be used also as a unit test
 14
 15import exudyn as exu
 16from exudyn.itemInterface import *
 17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 18import exudyn.graphics as graphics #only import if it does not conflict
 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.2            #driving torque; Nm ; 0.1Nm does not surpass critical speed; 0.2Nm works
 50eps = 2e-3*0.74         # 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
 54omegaInitial = 0*4*omega0 #initial rotation speed in rad/s
 55
 56tEnd = 200              #end time of simulation
 57steps = 40000           #number of steps
 58
 59fRes = omega0/(2*np.pi)
 60print('symmetric rotor resonance frequency (Hz)= '+str(fRes))
 61#print('runup over '+str(tEnd)+' seconds, fStart='+str(f0)+'Hz, fEnd='+str(f1)+'Hz')
 62
 63#linear frequency sweep evaluated at time t, for time interval [0, t1] and frequency interval [f0,f1];
 64def Sweep(t, t1, f0, f1):
 65    k = (f1-f0)/t1
 66    return np.sin(2*np.pi*(f0+k*0.5*t)*t) #take care of factor 0.5 in k*0.5*t, in order to obtain correct frequencies!!!
 67def SweepCos(t, t1, f0, f1):
 68    k = (f1-f0)/t1
 69    return np.cos(2*np.pi*(f0+k*0.5*t)*t) #take care of factor 0.5 in k*0.5*t, in order to obtain correct frequencies!!!
 70
 71# #user function for load
 72# def userLoad(t, load):
 73#     #return load*np.sin(0.5*omega0*t) #gives resonance
 74#     if t>40: time.sleep(0.02) #make simulation slower
 75#     return load*Sweep(t, tEnd, f0, f1)
 76#     #return load*Sweep(t, tEnd, f1, f0) #backward sweep
 77
 78# #backward whirl excitation:
 79# amp = 0.10  #in resonance: *0.01
 80# def userLoadBWy(t, load):
 81#     return load*SweepCos(t, tEnd, f0, f1) #negative sign: BW, positive sign: FW
 82# def userLoadBWz(t, load):
 83#     return load*Sweep(t, tEnd, f0, f1)
 84#def userLoadBWx(t, load):
 85#    return load*np.sin(omegaInitial*t)
 86#def userLoadBWy(t, load):
 87#    return -load*np.cos(omegaInitial*t) #negative sign: FW, positive sign: BW
 88
 89#background1 = graphics.BrickXYZ(0,0,0,.5,0.5,0.5,[0.3,0.3,0.9,1])
 90
 91#draw RGB-frame at origin
 92p=[0,0,0]
 93lFrame = 0.8
 94tFrame = 0.01
 95backgroundX = graphics.Cylinder(p,[lFrame,0,0],tFrame,[0.9,0.3,0.3,1],12)
 96backgroundY = graphics.Cylinder(p,[0,lFrame,0],tFrame,[0.3,0.9,0.3,1],12)
 97backgroundZ = graphics.Cylinder(p,[0,0,lFrame],tFrame,[0.3,0.3,0.9,1],12)
 98#mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [backgroundX, backgroundY, backgroundZ])))
 99
100#rotor is rotating around x-axis
101ep0 = eulerParameters0 #no rotation
102ep_t0 = AngularVelocity2EulerParameters_t([omegaInitial,0,0], ep0)
103print(ep_t0)
104
105p0 = [L0-0.5*L,eps,0] #reference position
106v0 = [0.,0.,0.] #initial translational velocity
107
108#node for Rigid2D body: px, py, phi:
109n1=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates = p0+ep0, initialVelocities=v0+list(ep_t0)))
110
111#ground nodes
112nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [-L/2,0,0]))
113nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [ L/2,0,0]))
114
115#add mass point (this is a 3D object with 3 coordinates):
116gRotor = graphics.Cylinder([-lRotor*0.5,0,0],[lRotor,0,0],r,[0.3,0.3,0.9,1],128)
117gRotor2 = graphics.Cylinder([-L0,0,0],[L,0,0],r*0.05,[0.3,0.3,0.9,1],16)
118gRotor3 = [backgroundX, backgroundY, backgroundZ]
119rigid = mbs.AddObject(RigidBody(physicsMass=m, physicsInertia=[Jxx,Jyyzz,Jyyzz,0,0,0], nodeNumber = n1, visualization=VObjectRigidBody2D(graphicsData=[gRotor, gRotor2]+gRotor3)))
120
121mbs.AddSensor(SensorBody(bodyNumber=rigid,
122                         fileName='solution/runupDisplacement.txt',
123                         outputVariableType=exu.OutputVariableType.Displacement))
124mbs.AddSensor(SensorBody(bodyNumber=rigid,
125                         fileName='solution/runupAngularVelocity.txt',
126                         outputVariableType=exu.OutputVariableType.AngularVelocity))
127
128#marker for ground (=fixed):
129groundMarker0=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround0))
130groundMarker1=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround1))
131
132#marker for rotor axis and support:
133rotorAxisMarker0 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[-L0,-eps,0]))
134rotorAxisMarker1 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[ L1,-eps,0]))
135
136
137#++++++++++++++++++++++++++++++++++++
138mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker0, rotorAxisMarker0],
139                                    stiffness=[k,k,k], damping=[d, d, d]))
140mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker1, rotorAxisMarker1],
141                                   stiffness=[0,k,k], damping=[0, d, d])) #do not constrain x-axis twice
142
143#coordinate markers for loads:
144rotorMarkerUy=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate=1))
145rotorMarkerUz=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate=2))
146
147#add torque:
148rotorRigidMarker =mbs.AddMarker(MarkerBodyRigid(bodyNumber=rigid, localPosition=[0,0,0]))
149mbs.AddLoad(Torque(markerNumber=rotorRigidMarker, loadVector=[torque,0,0]))
150
151#print(mbs)
152mbs.Assemble()
153#mbs.systemData.Info()
154
155simulationSettings = exu.SimulationSettings()
156simulationSettings.solutionSettings.solutionWritePeriod = 1e-5  #output interval
157simulationSettings.solutionSettings.sensorsWritePeriod = 1e-5  #output interval
158
159if isSymmetric:
160    simulationSettings.solutionSettings.solutionInformation = "Runup of Laval rotor, resonance="+str(round(fRes,3))+"Hz at 80-90 seconds"
161else:
162    simulationSettings.solutionSettings.solutionInformation = "Runup of unsymmetric rotor"
163
164simulationSettings.timeIntegration.numberOfSteps = steps
165simulationSettings.timeIntegration.endTime = tEnd
166simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
167simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
168
169simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
170
171#create animations (causes slow simulation):
172createAnimation=True
173if createAnimation:
174    simulationSettings.solutionSettings.recordImagesInterval = 0.2
175    SC.visualizationSettings.exportImages.saveImageFileName = "images/frame"
176    SC.visualizationSettings.window.renderWindowSize = [1600,1080]
177
178
179exu.StartRenderer()              #start graphics visualization
180mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
181
182#start solver:
183mbs.SolveDynamic(simulationSettings)
184
185#SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
186exu.StopRenderer()               #safely close rendering window!
187
188#evaluate final (=current) output values
189u = mbs.GetNodeOutput(n1, exu.OutputVariableType.AngularVelocity)
190print('omega=',u)
191#print('displacement=',u[0])
192
193#+++++++++++++++++++++++++++++++++++++++++++++++++++++
194import matplotlib.pyplot as plt
195import matplotlib.ticker as ticker
196
197if False:
198    plt.close('all') #close all plots
199
200    dataDisp = np.loadtxt('solution/runupDisplacement.txt', comments='#', delimiter=',')
201    dataOmega = np.loadtxt('solution/runupAngularVelocity.txt', comments='#', delimiter=',')
202
203    plt.plot(dataDisp[:,0], dataDisp[:,3], 'b-') #numerical solution
204    plt.xlabel("time (s)")
205    plt.ylabel("z-displacement (m)")
206
207    plt.figure()
208    plt.plot((1/(2*np.pi))*dataOmega[:,1], dataDisp[:,3], 'b-') #numerical solution
209    plt.xlabel("angular velocity (1/s)")
210    plt.ylabel("z-displacement (m)")
211
212    plt.figure()
213    plt.plot(dataOmega[:,0], (1/(2*np.pi))*dataOmega[:,1], 'b-') #numerical solution
214    plt.xlabel("time (s)")
215    plt.ylabel("angular velocity (1/s)")
216
217    plt.figure()
218    plt.plot(dataDisp[:,2], dataDisp[:,3], 'r-') #numerical solution
219    plt.xlabel("y-displacement (m)")
220    plt.ylabel("z-displacement (m)")
221
222    #plt.plot(data[n-500:n-1,1], data[n-500:n-1,2], 'r-') #numerical solution
223
224    ax=plt.gca() # get current axes
225    ax.grid(True, 'major', 'both')
226    ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
227    ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
228    plt.tight_layout()
229    plt.show()