flexibleRotor3Dtest.py
You can view and download this file on Github: flexibleRotor3Dtest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Flexible rotor test using two rigid bodies connected by 4 springs (corotating)
5# This test shows the unstable behavior if inner damping is larger than outer damping
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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14import exudyn as exu
15from exudyn.itemInterface import *
16from exudyn.utilities import *
17
18import time
19import numpy as np
20
21SC = exu.SystemContainer()
22mbs = SC.AddSystem()
23print('EXUDYN version='+exu.GetVersionString())
24
25L=1 #total rotor axis length
26m = 1 #mass of one disc in kg
27r = 0.5 #radius for disc mass distribution
28lRotor = 0.1 #length of one half rotor disk
29k = 800 #stiffness of (all/both) springs in rotor in N/m
30Jxx = 0.5*m*r**2 #polar moment of inertia
31Jyyzz = 0.25*m*r**2 + 1/12.*m*lRotor**2 #moment of inertia for y and z axes
32
33omega0=np.sqrt(2*k/(2*m)) #linear system; without flexibility of rotor
34
35#case 1: external damping: D0=0.002, D0int=0
36#case 2: external damping with small internal damping: D0=0.002, D0int=0.001
37#case 3: external damping with larger internal damping: D0=0.002, D0int=0.1
38#case 4: no external damping with small internal damping: D0=0, D0int=0.001
39attr = 'g-' #color in plot
40D0 = 0.002 #0.002 default; dimensionless damping
41D0int = 0.001*200 #*200 #default 0.001; dimensionless damping (not fully); value > 0.08 gives instability
42
43d = 2*omega0*D0*(2*m) #damping constant for external damping in N/(m/s)
44
45kInt = 4*800 #stiffness of (all/both) springs in rotor in N/m
46omega0int = np.sqrt(kInt/m)
47dInt = 2*omega0int*D0int*m #damping constant in N/(m/s)
48
49f0 = 0*omega0/(2*np.pi) #frequency start (Hz)
50f1 = 2.*omega0/(2*np.pi) #frequency end (Hz)
51
52torque = 0.5 #driving torque; Nm
53eps = 2e-3 #excentricity of mass in y-direction
54omegaInitial = 0*4*omega0 #initial rotation speed in rad/s
55
56print('resonance frequency (rad/s)= '+str(omega0))
57tEnd = 80 #end time of simulation
58steps = 10000 #number of steps
59
60
61#draw RGB-frame at origin
62p=[0,0,0]
63lFrame = 0.8
64tFrame = 0.01
65backgroundX = GraphicsDataCylinder(p,[lFrame,0,0],tFrame,[0.9,0.3,0.3,1],12)
66backgroundY = GraphicsDataCylinder(p,[0,lFrame,0],tFrame,[0.3,0.9,0.3,1],12)
67backgroundZ = GraphicsDataCylinder(p,[0,0,lFrame],tFrame,[0.3,0.3,0.9,1],12)
68#mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [backgroundX, backgroundY, backgroundZ])))
69
70#rotor is rotating around x-axis
71ep0 = eulerParameters0 #no rotation
72ep_t0 = AngularVelocity2EulerParameters_t([omegaInitial,0,0], ep0)
73print(ep_t0)
74
75p0 = [-lRotor*0.5,eps,0] #reference position
76p1 = [ lRotor*0.5,eps,0] #reference position
77v0 = [0.,0.,0.] #initial translational velocity
78
79#node for Rigid2D body: px, py, phi:
80n0=mbs.AddNode(RigidEP(referenceCoordinates = p0+ep0, initialVelocities=v0+list(ep_t0)))
81n1=mbs.AddNode(RigidEP(referenceCoordinates = p1+ep0, initialVelocities=v0+list(ep_t0)))
82
83#ground nodes
84nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [-L/2,0,0]))
85nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [ L/2,0,0]))
86
87#add mass point (this is a 3D object with 3 coordinates):
88transl = 0.9 #<1 gives transparent object
89gRotor0 = GraphicsDataCylinder([-lRotor*0.5,0,0],[lRotor,0,0],r,[0.3,0.3,0.9,transl],32)
90gRotor1 = GraphicsDataCylinder([-lRotor*0.5,0,0],[lRotor,0,0],r,[0.9,0.3,0.3,transl],32)
91gRotor0Axis = GraphicsDataCylinder([-L*0.5+0.5*lRotor,0,0],[L*0.5,0,0],r*0.05,[0.3,0.3,0.9,1],16)
92gRotor1Axis = GraphicsDataCylinder([-0.5*lRotor,0,0],[L*0.5,0,0],r*0.05,[0.3,0.3,0.9,1],16)
93gRotorCS = [backgroundX, backgroundY, backgroundZ]
94rigid0 = mbs.AddObject(RigidBody(physicsMass=m, physicsInertia=[Jxx,Jyyzz,Jyyzz,0,0,0], nodeNumber = n0, visualization=VObjectRigidBody2D(graphicsData=[gRotor0, gRotor0Axis]+gRotorCS)))
95rigid1 = mbs.AddObject(RigidBody(physicsMass=m, physicsInertia=[Jxx,Jyyzz,Jyyzz,0,0,0], nodeNumber = n1, visualization=VObjectRigidBody2D(graphicsData=[gRotor1, gRotor1Axis]+gRotorCS)))
96
97#marker for ground (=fixed):
98groundMarker0=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround0))
99groundMarker1=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround1))
100
101#marker for rotor axis and support:
102rotorAxisMarker0 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid0, localPosition=[-0.5*L+0.5*lRotor,-eps,0]))
103rotorAxisMarker1 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid1, localPosition=[ 0.5*L-0.5*lRotor,-eps,0]))
104
105
106#++++++++++++++++++++++++++++++++++++
107#supports:
108mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker0, rotorAxisMarker0],
109 stiffness=[k,k,k], damping=[d, d, d]))
110mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker1, rotorAxisMarker1],
111 stiffness=[0,k,k], damping=[0, d, d])) #do not constrain x-axis twice
112
113#++++++++++++++++++++++++++++++++++++
114#flexible rotor:
115nSprings = 4
116for i in range(nSprings):
117 #add corresponding markers
118 phi = 2*np.pi*i/nSprings
119 rSpring = 0.5
120 yPos = rSpring*np.sin(phi)
121 zPos = rSpring*np.cos(phi)
122 rotorM0 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid0, localPosition=[ 0.5*lRotor,yPos,zPos]))
123 rotorM1 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid1, localPosition=[-0.5*lRotor,yPos,zPos]))
124
125 mbs.AddObject(CartesianSpringDamper(markerNumbers=[rotorM0, rotorM1],
126 stiffness=[kInt,kInt,kInt], damping=[dInt, dInt, dInt]))
127
128#coordinate markers for loads:
129rotorMarkerUy=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate=1))
130rotorMarkerUz=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate=2))
131
132#add torque:
133rotorRigidMarker =mbs.AddMarker(MarkerBodyRigid(bodyNumber=rigid0, localPosition=[0,0,0]))
134mbs.AddLoad(Torque(markerNumber=rotorRigidMarker, loadVector=[torque,0,0]))
135
136#print(mbs)
137mbs.Assemble()
138#mbs.systemData.Info()
139
140simulationSettings = exu.SimulationSettings()
141simulationSettings.solutionSettings.solutionWritePeriod = 1e-5 #output interval
142simulationSettings.timeIntegration.numberOfSteps = steps
143simulationSettings.timeIntegration.endTime = 30#tEnd
144simulationSettings.timeIntegration.newton.useModifiedNewton=True
145simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
146simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
147simulationSettings.timeIntegration.verboseMode = 1
148simulationSettings.displayStatistics = True
149simulationSettings.displayComputationTime = True
150simulationSettings.linearSolverType = exu.LinearSolverType.EXUdense
151
152simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
153SC.visualizationSettings.general.useMultiThreadedRendering = False
154
155exu.StartRenderer() #start graphics visualization
156mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
157
158#start solver:
159mbs.SolveDynamic(simulationSettings)
160
161SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
162exu.StopRenderer() #safely close rendering window!
163
164#evaluate final (=current) output values
165u = mbs.GetNodeOutput(n1, exu.OutputVariableType.AngularVelocity)
166print('omega=',u)
167
168
169##+++++++++++++++++++++++++++++++++++++++++++++++++++++
170import matplotlib.pyplot as plt
171import matplotlib.ticker as ticker
172
173if True:
174 data = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
175 n=steps
176 plt.rcParams.update({'font.size': 24})
177
178 plt.plot(data[:,0], data[:,3], 'r-') #numerical solution
179
180 ax=plt.gca() # get current axes
181 ax.grid(True, 'major', 'both')
182 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
183 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
184 plt.tight_layout()
185 plt.show()