lavalRotor2Dtest.py
You can view and download this file on Github: lavalRotor2Dtest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example with 2D laval rotor; shows backward and forward whirl effects using
5# Includes user load and Sweep
6#
7# Author: Johannes Gerstmayr
8# Date: 2019-12-03
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 length
28mass = 1.6 #mass in kg
29r = 0.5 #radius for disc mass distribution
30spring = 4000 #stiffness of (all/both) springs in rotor in N/m
31omega0=np.sqrt(spring/mass) #linear system
32
33zeta = 0.001*5
34damper = 2*omega0*zeta*mass #damping constant in N/(m/s)
35
36f0 = 0.*omega0/(2*np.pi)
37f1 = 1.*omega0/(2*np.pi)
38
39torque = 0*1 #Nm
40eps = 1e-2 #excentricity in y-direction
41omegaInitial = 0.5*omega0
42
43print('resonance frequency = '+str(omega0/2/np.pi)+'Hz')
44tEnd = 50 #end time of simulation
45steps = 10000 #number of steps
46
47
48#linear frequency sweep in time interval [0, t1] and frequency interval [f0,f1];
49def Sweep(t, t1, f0, f1):
50 k = (f1-f0)/t1
51 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!!!
52
53#user function for load
54def userLoad(mbs, t, load):
55 #return load*np.sin(0.5*omega0*t) #gives resonance
56 if t>40: time.sleep(0.02) #make simulation slower
57 return load*Sweep(t, tEnd, f0, f1)
58 #return load*Sweep(t, tEnd, f1, f0) #backward sweep
59
60#backward whirl excitation:
61amp = 0*10 #in resonance: *0.01
62def userLoadBWx(mbs, t, load):
63 return load*np.sin(omegaInitial*t)
64def userLoadBWy(mbs, t, load):
65 return -load*np.cos(omegaInitial*t) #negative sign: FW, positive sign: BW
66
67#node for Rigid2D body: px, py, phi:
68n1=mbs.AddNode(Rigid2D(referenceCoordinates = [0,eps,0],
69 initialCoordinates=[0,0,0],
70 initialVelocities=[0,0,omegaInitial]))
71
72#ground node
73nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
74
75#add mass point (this is a 3D object with 3 coordinates):
76gRotor = GraphicsDataRectangle(-r*0.5,-r*0.5,r*0.5,r*0.5,[1,0,0,1])
77rigid2D = mbs.AddObject(RigidBody2D(physicsMass=mass, physicsInertia=mass*r**2, nodeNumber = n1, visualization=VObjectRigidBody2D(graphicsData=[gRotor])))
78
79#marker for ground (=fixed):
80groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
81
82#marker for rotor axis and support:
83rotorAxisMarker =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid2D, localPosition=[0,-eps,0]))
84groundBody= mbs.AddObject(ObjectGround(referencePosition=[0,0,0]))
85rotorSupportMarker=mbs.AddMarker(MarkerBodyPosition(bodyNumber=groundBody, localPosition=[0,0,0]))
86
87#marker for springDamper for first (x-)coordinate:
88coordXMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
89coordYMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 1))
90coordPhiMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 2))
91
92mbs.AddObject(CartesianSpringDamper(markerNumbers=[rotorSupportMarker, rotorAxisMarker], stiffness=[spring,spring,0], damping=[damper, damper,0]))
93
94#add load:
95mbs.AddLoad(LoadCoordinate(markerNumber = coordYMarker, load = 0, loadUserFunction=userLoad))
96mbs.AddLoad(LoadCoordinate(markerNumber = coordPhiMarker, load = torque))#, loadUserFunction=userLoad))
97
98mbs.AddLoad(LoadCoordinate(markerNumber = coordXMarker, load = amp, loadUserFunction=userLoadBWx))
99mbs.AddLoad(LoadCoordinate(markerNumber = coordYMarker, load = amp, loadUserFunction=userLoadBWy))
100
101print(mbs)
102mbs.Assemble()
103
104simulationSettings = exu.SimulationSettings()
105simulationSettings.solutionSettings.solutionWritePeriod = 1e-5 #output interval
106simulationSettings.timeIntegration.numberOfSteps = steps
107simulationSettings.timeIntegration.endTime = tEnd
108
109simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
110
111exu.StartRenderer() #start graphics visualization
112mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
113
114#start solver:
115mbs.SolveDynamic(simulationSettings)
116
117exu.StopRenderer() #safely close rendering window!
118
119#evaluate final (=current) output values
120u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
121print('displacement=',u[0])
122
123#+++++++++++++++++++++++++++++++++++++++++++++++++++++
124import matplotlib.pyplot as plt
125import matplotlib.ticker as ticker
126
127if True:
128 data = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
129 n=steps
130 #plt.plot(data[:,0], data[:,6], 'r-') #numerical solution
131 plt.plot(data[n-500:n-1,1], data[n-500:n-1,2], 'r-') #numerical solution
132
133 ax=plt.gca() # get current axes
134 ax.grid(True, 'major', 'both')
135 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
136 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
137 plt.tight_layout()
138 plt.show()