Torsional dynamics

[LeifParr]

I am currently struggling because our commercial software which we use to model torsional natural frequencies and forced response does not support time transient analysis.

It seems that there is a possibility in ROSS for a 6dof shaft element, but I do not immediately see how I could run a torsional analysis using the ROSS package. Is it possible today? Could it be possible in the future?

I'm trying to set up a generator short-circuit simulation of a directly coupled generator-flywheel-turbine-unit. So it could be a very simple model, maybe a 5-inertia lumped model would do. Any thoughts here?

[raphaeltimbo]

Hi @LeifParr !
We currently have the shaft and a general bearing element implemented with 6 dof.
The process to instantiate these elements is very similar to the process for 4 dof.
To instantiate the rotor you just pass the lists with the elements to the Rotor class (just like 4 dof), making sure that all elements that you are passing are 6dof (ROSS will raise an error if you mix them).
A colleague at work caried out a torsional analysis using ROSS, but I don’t know all details and how complex the analysis was and in which parts of ROSS he might have found problems.
I know that our Campbell Diagram plot for example does not work now in the 6 dof rotor.
It is in our plan to improve the support for 6dof, but I don’t know when we will have the time to do it.

[atbrandao]

Hi, @LeifParr !

In the torsional analysis I carried out using Ross I was only interested in the system's harmonic response, so I used only the Transfer Matrix to calculate the system output.
But you can use the .time_transient() method to perform this time transient analysis. The only limitation, I think, is that there are no native method for plotting the results. You will need to get the output data and make your own plots.

I made a simple code (attached) with a model similar to what you proposed, with 4 shaft elements connecting 5 lumped inertia elements. The code inputs sine excitation at node 0 and plots the torsional response at node 4, in radians. With that result you could then do any analysis you wish.

import ross as rs
import matplotlib.pyplot as plt
import numpy as np

Nel = 4 # Number of shaft elements
shaft = [
    rs.ShaftElement6DoF(
        n=i,
        L=0.05,
        odl=0.05,
        idl=0,
        odr=0.05,
        idr=0,
        material=rs.materials.steel,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(Nel)
]

disks = [rs.DiskElement6DoF(Id=1.7809,Ip=0.32956,m=32.59,n=0),
         rs.DiskElement6DoF(Id=1.7809,Ip=0.32956,m=32.59,n=1),
         rs.DiskElement6DoF(Id=1.7809,Ip=0.32956,m=32.59,n=2),
         rs.DiskElement6DoF(Id=1.7809,Ip=0.32956,m=32.59,n=3),
         rs.DiskElement6DoF(Id=1.7809,Ip=0.32956,m=32.59,n=4),]

# Bearings to be created only to fix X Y coordinates. Arbitrary coefficients
brg1 = rs.BearingElement6DoF(n=0,
                        kxx=[1e9], kxy=[0],
                        kyx=[0], kyy=[1e9],
                        cxx=[0], cxy=[0],
                        cyx=[0], cyy=[0],
                        frequency=None,)
brg2 = rs.BearingElement6DoF(n=4,
                        kxx=[1e9], kxy=[0],
                        kyx=[0], kyy=[1e9],
                        cxx=[0], cxy=[0],
                        cyx=[0], cyy=[0],
                        frequency=None,)

rotor = rs.Rotor(shaft,disks,[brg1,brg2])

f = 600.0 
time_samples = 1001
node_in = 0 # node where torsional force will be applied
t = np.linspace(0, 16, time_samples)

F = np.zeros((time_samples, rotor.ndof))

# Input force time array
F[:, 6 * node_in + 5] = 10 * np.cos(f * t) # The "5" indicates the torsional DoF

response = rotor.time_response(0, F, t) # "0" stands for the speed, which is irrelevant for a purely torsional analysis

theta = response[2][:,5::6] # Gets the response only of the torsional DoF

node_out = 4 # Node to plot the response 
plt.plot(t,theta[:,node_out])
# Use 'node_out + 5' to plot angular velocity vs time

[LeifParr]

Thank you for feedback here, and the example, it got me going pretty quick.
I ran through my 3-node example, and some more questions came up.

Firstly, I wanted to verify the torsional natural frequencies versus other software. So I ran rotor.run_modal(speed=0.0)
I had results at 1e-5 Hz, 27 Hz and 90 Hz which were pretty close to output of other software. However, there were also other natural frequencies of non-torsional modes (35 Hz in this case).

How would I go about chasing down the purely torsional modes?

The second question would be to produce similar output as other software, both angle and moment. Is the moment a nodal property of the output, or do we have to calculate moment for the shafts based on angle of neighboring nodes and stiffness of shaft? I tried to chase down the calculation of the rotor.time_response, but it seems to run scipy.lti almost directly, and it kind of got overwhelming for me.

Then, there is the notion of damping. Using RAPPID for harmonic response, damping per mode is set in %. Other software uses damping per node. I do not see how damping could be included here?

The ultimate goal for me here would be to come up with simple routines to get rid of my existing software for torsional analysis and replace with Python + ROSS:

• torsional natural frequencies (I can make plotting routines)
• harmonic response
• transient response

[atbrandao]

Hi, @LeifParr
Thank you for your input! Good to see that you were able to move forward with your simulation.

The points that you raised are really interesting and certainly will drive improvements on the software. I think that there are no built-in ways that solve your problems (@raphaeltimbo can correct me if I'm wrong), specially because most of the software was developed for 4DoF elements, and the torsional analysis capability was a later development.

However, we can work a few ways to solve this issues right before further improvements are implemented.

Separating modes in a multi-DoF is fundamentaly complicated, because mathematically they are all simply eigenvectors of the system, with no distiction between torsional or lateral modes.
The only way to separate the modes is to take a look at the mode shapes and choose those with greater displacement of the torsional DoF.
However, to make the separation of the lateral modes more clear, you can add bearing with very high stiffness at every node.
I would do something like the code below to create the rotor and plot the desired mode:

brgs = [rs.BearingElement6DoF(n=i,
                        kxx=[1e20], kxy=[0],
                        kyx=[0], kyy=[1e20],
                        cxx=[0], cxy=[0],
                        cyx=[0], cyy=[0],
                        frequency=None,) for i in range(Nel+1)]

rotor2 = rs.Rotor(shaft,disks,brgs)

modal = rotor2.run_modal(0,num_modes=rotor2.ndof)

tors_modes_id = []

for i in range(len(modal.modes[0,:])):
    aux = []
    for j in range(6):
        aux.append(np.sum(np.abs(modal.modes[j::6,i])))
      
    if  np.argmax(aux) == 5:   
        tors_modes_id.append(i)
    
i_sort = np.argsort(np.abs(np.imag(modal.evalues[tors_modes_id])))

wn = np.abs(np.imag(modal.evalues[tors_modes_id]))[i_sort]
tors_modes = modal.modes[5::6,tors_modes_id]
tors_modes = tors_modes[:,i_sort]

n_mode = 5

re = np.real(tors_modes[:,n_mode])
imag = np.imag(tors_modes[:,n_mode])

plt.plot(np.sign(re+imag)*np.abs(tors_modes[:,n_mode]))

For the second question, about plotting the moment, I don't know of a built-in way of calculation that for the torsional DoF.
However, it is pretty easy to do using the .K() method for either the ShaftElement object or the Rotor object.
In the code below you get the Torque at the chosen element or node for a specific theta vector, as defined in the transient analysis in my previous code suggestion.

node = 2
T_node = rotor2.K(0)[node*6+5,5::6] @ theta[:,:len(rotor2.nodes)].transpose()

el = 3
T_element = (rotor.shaft_elements[el].K()[5::6,5::6] @ theta[:,el:el+2].transpose())[0,:]

The third question, about damping, is a bit trickier. In the torsional analysis I made, I artificially manipulates the M, K and C matrices in the rotor_assembly.py file to build my model, and then I added damping in this way.

The best way to do it for a typical rotating machinery would be to add a concentrated damping on the coupling element until it gets your first vibration mode to a damping factor around 1-2%, a value suggested by the API STD 684 depending on the coupling type.
To do that I would alter the .C() method in the rotor_assembly.py file in the following way, which would add a proportinal damping to the element corresponding to the coupling.

C0 = np.zeros((self.ndof, self.ndof))

coup_el = 3 # Element number corresponding to the coupling
coef = 1e-4 # Adjust this value until you get the desired damping factor

for i, elm in enumerate(self.elements):
    dofs = elm.dof_global_index
    if i==coup_el:
        try:
            C0[np.ix_(dofs, dofs)] += coef*elm.K(frequency)
        except TypeError:
            C0[np.ix_(dofs, dofs)] += coef*elm.K()
    else:
        try:
            C0[np.ix_(dofs, dofs)] += elm.C(frequency)
        except TypeError:
            C0[np.ix_(dofs, dofs)] += elm.C()

return C0

Thank you again for your input, and I hope these tips can help you.

[LeifParr]

I am moving forward with the transient torsional, however at a very slow pace. I am starting to get a fairly good match between ROSS and XLrotor for my 3-inertia simple example. Natural frequencies are equal, I have added damping to get about 1% damping ratio at the first natural frequency as per your suggestion above.
XLrotor transient:

Ross forcing function -

Ross response in turbine node.

However, there is one thing I am left wondering about - In this line where we do time_response
response = rotor.time_response(0, F, t) # "0" stands for the speed, which is irrelevant for a purely torsional analysis
theta = response[2][:,5::6] # Gets the response only of the torsional DoF

Here it seems the index response[2] gives the xout from time_response, while we use the yout for other plotting routines. Is this correct, or should we use the yout? e.g.
t_, yout, xout = rotor.time_response(0, F, t)
theta = yout[:,5::6]

I seem to be getting the same results anyway, so maybe it isn't important.

So far, running transient analysis for short-circuit gives a significantly greater moment in the coupling than running a harmonic response at 50 Hz and adding nominal torque.

[atbrandao]

Good to see you are getting good results with the ROSS model!

As for the question of yout vs xout, the difference between them is that xout shows the evolution of the system's whole state vector, i.e. displacements and velocities at each point in time. The yout variable only gives you the displacements.

So, if you have a 18 DoF system, and a simulation of 100 points in time, yout will be a 18x100 matrix and xout will be a 36x100 matrix, and yout = xout[ :18 , : ]. The second half of xout rows represents the velocities of each node.

If you are interested on the displacements only, and not in nodal angular velocities, you can use only yout, no problem.

The result of much greater moments in the transient analysis, if compared to the harmonic analysis, makes sense for me. The transient excitation ends up excitation the system's free vibration response, which in weakly damped systems such as this one can cause high displacement amplitudes.