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