rigidRotor3Dnutation.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example with 3D rotor, test nutation with point force
  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
 27m = 2                   #mass in kg
 28r = 0.5                 #radius for disc mass distribution
 29lRotor = 0.2            #length of rotor disk
 30k = 8000                 #stiffness of (all/both) springs in rotor in N/m
 31Jxx = 0.5*m*r**2        #polar moment of inertia
 32Jyyzz = 0.25*m*r**2 + 1/12.*m*lRotor**2      #moment of inertia for y and z axes
 33
 34omega0=np.sqrt(2*k/m) #linear system
 35
 36D0 = 0.002              #dimensionless damping
 37d = 2*omega0*D0*m       #damping constant in N/(m/s)
 38
 39omegaInitial = 0.1*omega0 #initial rotation speed in rad/s
 40
 41print('resonance frequency (rad/s)= '+str(omega0))
 42
 43tEnd = 50               #end time of simulation
 44steps = 20000         #number of steps
 45
 46
 47#user function for load
 48def userLoad(mbs, t, load):
 49    #time.sleep(0.005) #make simulation slower
 50    if t<0.01: print(load)
 51    if t>10 and t<10.05:
 52        return load
 53    else:
 54        return [0,0,0]
 55
 56
 57#draw RGB-frame at origin
 58p=[0,0,0]
 59lFrame = 0.8
 60tFrame = 0.01
 61backgroundX = graphics.Cylinder(p,[lFrame,0,0],tFrame,[0.9,0.3,0.3,1],12)
 62backgroundY = graphics.Cylinder(p,[0,lFrame,0],tFrame,[0.3,0.9,0.3,1],12)
 63backgroundZ = graphics.Cylinder(p,[0,0,lFrame],tFrame,[0.3,0.3,0.9,1],12)
 64#mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [backgroundX, backgroundY, backgroundZ])))
 65
 66#rotor is rotating around x-axis
 67ep0 = eulerParameters0 #no rotation
 68ep_t0 = AngularVelocity2EulerParameters_t([omegaInitial,0,0], ep0)
 69print(ep_t0)
 70
 71p0 = [0,0,0] #reference position
 72v0 = [0.,0.,0.] #initial translational velocity
 73
 74#node for Rigid2D body: px, py, phi:
 75n1=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates = p0+ep0, initialVelocities=v0+list(ep_t0)))
 76
 77#ground nodes
 78nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 79
 80#add mass point (this is a 3D object with 3 coordinates):
 81gRotor = graphics.Cylinder([-lRotor*0.5,0,0],[lRotor*0.5,0,0],r,[0.3,0.3,0.9,1],32)
 82gRotor3 = [backgroundX, backgroundY, backgroundZ]
 83rigid = mbs.AddObject(RigidBody(physicsMass=m, physicsInertia=[Jxx,Jyyzz,Jyyzz,0,0,0], nodeNumber = n1,
 84                                visualization=VObjectRigidBody2D(graphicsData=[gRotor]+gRotor3)))
 85
 86#marker for ground (=fixed):
 87groundMarker0=mbs.AddMarker(MarkerNodePosition(nodeNumber= nGround0))
 88
 89#marker for rotor axis and support:
 90rotorAxisMarker0 =mbs.AddMarker(MarkerBodyPosition(bodyNumber=rigid, localPosition=[0,0,0]))
 91
 92
 93#++++++++++++++++++++++++++++++++++++
 94mbs.AddObject(CartesianSpringDamper(markerNumbers=[groundMarker0, rotorAxisMarker0],
 95                                    stiffness=[k,k,k], damping=[d, d, d]))
 96
 97#add force/torque:
 98rotorRigidMarker =mbs.AddMarker(MarkerBodyRigid(bodyNumber=rigid, localPosition=[0,r,0]))
 99mbs.AddLoad(Force(markerNumber=rotorRigidMarker, loadVector=[0.3,0.2,0.1], loadVectorUserFunction = userLoad))
100#mbs.AddLoad(Torque(markerNumber=rotorRigidMarker, loadVector=[torque,0,0]))
101
102#print(mbs)
103mbs.Assemble()
104#mbs.systemData.Info()
105
106simulationSettings = exu.SimulationSettings()
107simulationSettings.solutionSettings.solutionWritePeriod = 1e-2  #output interval
108simulationSettings.timeIntegration.numberOfSteps = steps
109simulationSettings.timeIntegration.endTime = tEnd
110simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
111simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
112
113simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
114
115
116#start solver:
117mbs.SolveDynamic(simulationSettings)
118
119exu.StartRenderer()              #start graphics visualization
120mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
121
122exu.StopRenderer()               #safely close rendering window!
123
124mbs.SolutionViewer()
125
126###+++++++++++++++++++++++++++++++++++++++++++++++++++++
127#import matplotlib.pyplot as plt
128#import matplotlib.ticker as ticker
129#
130#if True:
131#    data = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
132#    n=steps
133#    #plt.plot(data[:,2], data[:,3], 'r-') #numerical solution
134#    #plt.plot(data[:,0], data[:,2], 'b-') #numerical solution
135#    plt.plot(data[:,0], data[:,3], 'g-') #numerical solution
136#    #plt.plot(data[n-500:n-1,1], data[n-500:n-1,2], 'r-') #numerical solution
137#
138#    ax=plt.gca() # get current axes
139#    ax.grid(True, 'major', 'both')
140#    ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
141#    ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
142#    plt.tight_layout()
143#    plt.show()