Tutorial - Running Analyzes#

This is the second part of a basic tutorial on how to use ROSS (rotordynamics open-source software), a Python library for rotordynamic analysis. In this tutorial, you will learn how to run several rotordynamic analyzes with your rotor model

To get results, we always have to use one of the .run_ methods available for a rotor object. These methods will return objects that store the analysis results and that also have plot methods available. These methods will use the plotly library to make graphs common to a rotordynamic analysis.

We can also use units when plotting results. For example, for a unbalance response plot we have the amplitude_units argument and we can choose between any length unit available in pint such as ‘meter’, ‘inch’, etc.

Rotor model#

First, let’s recover the rotor model built in the previous tutorial.

from pathlib import Path
import ross as rs
import numpy as np
import plotly.graph_objects as go
import plotly.io as pio

pio.renderers.default = "notebook"

shaft_file = Path("shaft_si.xls")
shaft = rs.ShaftElement.from_table(
    file=shaft_file, sheet_type="Model", sheet_name="Model"
)

file_path = Path("shaft_si.xls")
list_of_disks = rs.DiskElement.from_table(file=file_path, sheet_name="More")

bearing1 = rs.BearingElement.from_table(n=7, file="bearing_seal_si.xls")
bearing2 = rs.BearingElement.from_table(n=48, file="bearing_seal_si.xls")

bearings = [bearing1, bearing2]

rotor3 = rs.Rotor(shaft, list_of_disks, bearings)

node_increment = 5
rotor3.plot_rotor(nodes=node_increment)

Rotor Analyses#

In the last tutorial we have learnt how to create a rotor model with Rotor class. Now, we’ll use the same class to run the simulation. There’re some methods, most of them with the prefix run_ you can use to run the rotordynamics analyses.

For Most of the methods, ypu can use the command .plot() to display a graphical visualization of the results (e.g run_campbel().plot(), run_freq_response().plot()).

ROSS offers the following analyses:

  • Static analysis

  • Modal analysis

  • Campbell Diagram

  • Undamped Critical Speed Map

  • Frequency response

  • Unbalance response

  • Time response

Plotly library#

ROSS uses Plotly for plotting results. All the figures can be stored and manipulated following Plotly API.

The following sections presents the results and how to return the Plotly Figures.

1.1 Static Analysis#

This method runs the static analysis for the rotor. It calculate the static deformation due the gravity effects (shaft and disks weight). It also returns the bending moment and shearing force on each node, and you can return a free-body-diagram representation for the rotor, with the self weight, disks weight and reaction forces on bearings displayed.

1.1.1 Running static analysis#

To run the simulation, use the .run_static() method. You can define a variable to store the results.

Storing the results, it’s possible to return the following arrays:

  • disk_forces_nodal

  • disk_forces_tag

  • bearing_forces_nodal

  • bearing_forces_tag

  • disp_y

  • Vx

  • Bm

static = rotor3.run_static()

Returning forces#

Disk forces#

  • .disk_forces_nodal: Returns a dictionaty expliciting the node where the disk is located and the force value.

  • .disk_forces_tag: Returns a dictionaty expliciting the the disk tag is located and the force value.

Bearing forces#

  • .bearing_forces_nodal: Returns a dictionaty expliciting the node where the bearing is located and the force value.

  • .bearing_forces_tag: Returns a dictionaty expliciting the the bearing tag is located and the force value.

print("Disk forces - nodes")
print(rotor3.disk_forces_nodal)
print("")
print("Disk forces - tags")
print(rotor3.disk_forces_tag)

print("")
print("Bearing forces - nodes")
print(rotor3.bearing_forces_nodal)
print("")
print("Bearing forces - tags")
print(rotor3.bearing_forces_tag)

Disk forces - nodes
{'node_3': 148.2741036647235, 'node_20': 67.76283441291268, 'node_23': 67.95896417966495, 'node_26': 68.15509394641724, 'node_29': 68.44928859654566, 'node_32': 68.0570290630411, 'node_35': 68.25315882979338}

Disk forces - tags
{'Disk 0': 148.2741036647235, 'Disk 1': 67.76283441291268, 'Disk 2': 67.95896417966495, 'Disk 3': 68.15509394641724, 'Disk 4': 68.44928859654566, 'Disk 5': 68.0570290630411, 'Disk 6': 68.25315882979338}

Bearing forces - nodes
{'node_7': 1216.2827567762095, 'node_48': 1204.651471281892}

Bearing forces - tags
{'Bearing 0': 1216.2827567762095, 'Bearing 1': 1204.651471281892}

Other attributes from static analysis#

  • .Vx: Shearing force array

  • .Bm: Bending moment array

  • .deformation Displacement in Y direction

1.1.2 Plotting results#

With results stored, you can use some methods to plot the results. Currently, there’re four plots you can retrieve from static analysis:

  • .plot_free_body_diagram()

  • .plot_deformation()

  • .plot_shearing_force()

  • .plot_bending_moment()

Plotting free-body-diagram#

static.plot_free_body_diagram()

Plotting deformation#

static.plot_deformation()

Plotting shearing force diagram#

static.plot_shearing_force()

Plotting bending moment diagram#

static.plot_bending_moment()

1.3 Campbell Diagram#

Also called “whirl speed map” in rotordynamics, ROSS calculate and plots the modes’ damped eigenvalues and the logarithmic decrement as a function of rotor speed.

To run the Campbell Diagram, use the command .run_campbell(). The user must input an array of speeds, which will be iterated to calculate each point on the graph.

This method returns the damped natural frequencies, logarithmic decrement and the whirl values (values indicating the whirl direction: backward or forward).

1.3.1 Running campbell diagram#

In this example the whirl speed map is calculated for a speed range from 0 to 1000 rad/s (~9550 RPM).

Storing the results, it’s possible to return the following arrays:

  • wd

  • log_dec

  • whirl_values

Each value in these arrays is calculated for each speed value in speed_range

samples = 31
speed_range = np.linspace(315, 1150, samples)

campbell = rotor3.run_campbell(speed_range)

# results for each frequency
frequency_index = 0
print(campbell.wd[:, frequency_index])

[4.05696248e-03 1.13308778e-02 2.87443814e-02 6.81662768e-02
 1.53837840e-01 3.34102882e-01 7.04288340e-01 1.45788447e+00
 3.03663623e+00 5.25443046e+00 7.21720675e+00 1.06507330e+01
 1.85851449e+01 1.01602531e+02 4.67359814e+02 6.34570746e+02
 6.35176758e+02 6.35970906e+02 6.36891161e+02 6.37880817e+02
 6.38897245e+02 6.39918527e+02 6.40931104e+02 6.41925672e+02
 6.42889625e+02 6.43808455e+02 6.44672982e+02 6.45478063e+02
 6.46221340e+02 6.46902349e+02 6.47521893e+02]

# results for each frequency
frequency_index = 1
print(campbell.log_dec[:, frequency_index])

[1.29458651e+03 1.24034622e+03 1.18625615e+03 1.13071170e+03
 1.07242137e+03 1.01021334e+03 9.42865193e+02 8.68904405e+02
 7.86381133e+02 7.82375744e+02 3.44920711e+02 3.94207681e+02
 2.88161979e+02 7.77554137e+01 1.82932326e+01 4.83302025e-01
 4.33157440e-01 3.88197688e-01 3.48380244e-01 3.13495797e-01
 2.83129231e-01 2.56672664e-01 2.33605438e-01 2.13500543e-01
 1.95967979e-01 1.80635634e-01 1.67181122e-01 1.55328894e-01
 1.44844484e-01 1.35528792e-01 1.27212823e-01]

1.3.2 Plotting results#

Now that the results are stored, use .plot() method to display the Campbell Diagram plot.

For the Campbell Diagram, you can plot more than one harmonic. As default, the plot display only the 1x speed option. Input a list with more harmonics to display it at the graph.

campbell.plot(harmonics=[0.5, 1])

1.4 Frequency Response#

ROSS’ method to calculate the Frequency Response Function is run_freq_response(). This method returns the magnitude and phase in the frequency domain. The response is calculated for each node from the rotor model.

When plotting the results, you can choose to plot:

  • amplitude vs frequency: plot_magnitude()

  • phase vs frequency: plot_phase()

  • polar plot of amplitude vs phase: plot_polar_bode()

  • all: plot()

1.4.1 Clustering points#

The number of solution points is an important parameter to determine the computational cost of the simulation. Besides the classical method, using numpy.linspace, which creates an evenly spaced array over a specified interval, ROSS offers an automatic method to create an speed_range array.

The method clustering_points generates an automatic array to run frequency response analyses. The frequency points are calculated based on the damped natural frequencies and their respective damping ratios. The greater the damping ratio, the more spread the points are. If the damping ratio, for a given critical speed, is smaller than 0.005, it is redefined to be 0.005 (for this method only).

The main goal of this feature is getting a more accurate amplitude value for the respective critical frequencies and nearby frequency points.

1.4.2 Running frequency response#

To run the this analysis, use the command run_freq_response(). You can give a specific speed_range or let the program run with the default options. In this case, no arguments are needed to input.

First, let’s run an example with a “user-defined” speed_range. Setting an array to speed_range will disable all the frequency spacing parameters.

samples = 61
speed_range = np.linspace(315, 1150, samples) # rads/s
results1 = rotor3.run_freq_response(speed_range=speed_range)

results1.speed_range.size

61

Now, let’s run an example using clustering points array.

# results1_2 = rotor2.run_freq_response(cluster_points=True, num_points=5, num_modes=12)

# results1_2.speed_range.size

In the next section we’ll check the difference between both results.

1.4.3 Plotting results - Bode Plot#

We can plot the frequency response selecting the input and output degree of freedom.

  • Input is the degree of freedom to be excited;

  • Output is the degree of freedom to be observed.

Each shaft node has 4 local degrees of freedom (dof) \([x, y, \alpha, \beta]\), and each degree of freedom has it own index:

  • \(x\) -> index 0

  • \(y\) -> index 1

  • \(\alpha\) -> index 2

  • \(\beta\) -> index 3

To select a DoF to input and a DoF to the output, we have to use the following correlation:

\(global\_dof = dof\_per\_node \cdot node\_number + dof\_index\)

For example: node 26, local dof \(y\):

\(DoF = 26 \cdot 4 + 1 = 25\)

plot = results1.plot(inp=105, out=105)
plot

# converting the first plot yaxis to log scale
# plot = results1.plot(inp=105, out=105)
# plot.update_yaxes(type="log", row=1, col=1)
# plot

# plot1_2 = results1_2.plot(inp=105, out=105)
# plot1_2

# # converting the first plot yaxis to log scale
# plot1_2 = results1_2.plot(inp=105, out=105)
# plot1_2.update_yaxes(type="log", row=1, col=1)
# plot1_2

1.5 Unbalance Response#

ROSS’ method to simulate the reponse to an unbalance is run_unbalance_response(). This method returns the unbalanced response in the frequency domain for a given magnitide and phase of the unbalance, the node where it’s applied and a frequency range.

ROSS takes the magnitude and phase and converts to a complex force array applied to the given node:

\[\begin{split} force = \left(\begin{array}{cc} F \cdot e^{j\delta}\\ -jF \cdot e^{j\delta}\\ 0\\ 0 \end{array}\right) \end{split}\]

where:

  • \(F\) is the unbalance magnitude;

  • \(\delta\) is the unbalance phase;

  • \(j\) is the complex number notation;

When plotting the results, you can choose to plot the:

  • Bode plot options for a single degree of freedom:

    • amplitude vs frequency: plot_magnitude()

    • phase vs frequency: plot_phase()

    • polar plot of amplitude vs phase: plot_polar_bode()

    • all: plot()

  • Deflected shape plot options:

    • deflected shape 2d: plot_deflected_shape_2d()

    • deflected shape 3d: plot_deflected_shape_3d()

    • bending moment: plot_bending_moment()

    • all: plot_deflected_shape()

run_unbalance_response() is also able to work with clustering points ( see section 7.4.1 ).

1.5.1 Running unbalance response#

To run the Unbalance Response, use the command .unbalance_response()

In this following example, we can obtain the response for a given unbalance and its respective phase in a selected node. Notice that it’s possible to add multiple unbalances instantiating node, magnitude and phase as lists.

The method returns the force response array (complex values), the displacement magnitude (absolute value of the forced response) and the phase of the forced response.

Let’s run an example with 2 unbalances in phase, trying to excite the first and the third natural vibration mode.

Unbalance1: node = 29
            magnitude = 0.03
            phase = 0
Unbalance2: node = 33
            magnitude = 0.02
            phase = 0

n1 = 29
m1 = 0.003
p1 = 0

n2 = 33
m2 = 0.002
p2 = 0

frequency_range=np.linspace(315, 1150, 101)
results2 = rotor3.run_unbalance_response([n1, n2], [m1, m2], [p1, p2], frequency_range)

1.5.2 Plotting results - Bode Plot#

To display the bode plot, use the command .plot(probe)

Where probe is a list of tuples that allows you to choose not only the node where to observe the response, but also the orientation.

Probe orientation equals 0º refers to +X direction (DoFX), and probe orientation equals 90º (or \(\frac{\pi}{2} rad\)) refers to +Y direction (DoFY).

You can insert multiple probes at once.

# probe = (probe_node, probe_orientation)
probe1 = (15, 45) # node 15, orientation 45º
probe2 = (35, 45) # node 35, orientation 45º

results2.plot(probe=[probe1, probe2], probe_units="degrees")

# converting the first plot yaxis to log scale
# plot2 = results2.plot(probe=[probe1, probe2], probe_units="rad")
# plot2.update_yaxes(type="log", row=1, col=1)
# plot2

1.5.3 Plotting results - Deflected shape#

To display the deflected shape configuration, use the command .plot_deflected_shape()

results2.plot_deflected_shape(speed=649)

1.6 Time Response#

ROSS’ method to calculate displacements due a force in time domain is run_time_response(). This function will take a rotor object and plot its time response given a force and a time array.

The force input must be a matrix \(M \times N\), where:

  • \(M\) is the size of the time array;

  • \(N\) is the rotor’s number of DoFs (you can access this value via attribute .ndof).

Each row from the matrix represents a node, and each column represents a time step.

Time Response allows you to plot the response for:

  • a list of probes

  • an orbit for a given node (2d plot)

  • all the nodes orbits (3d plot)

1.6.1 Running time response#

To run the Time Response, use the command .run_time_response().

Building the force matrix is not trivial. We recommend creating a matrix of zeros using numpy.zeros() and then, adding terms to the matrix.

In this examples, let’s create an harmonic force on node 26, in \(x\) and \(y\) directions (remember index notation from Frequency Response (section 7.4.2). We’ll plot results from 0 to 10 seconds of simulation.

speed = 600.0
time_samples = 1001
node = 26
t = np.linspace(0, 16, time_samples)

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

# component on direction x
F[:, 4 * node + 0] = 10 * np.cos(2 * t)
# component on direction y
F[:, 4 * node + 1] = 10 * np.sin(2 * t)

response3 = rotor3.run_time_response(speed, F, t)

1.6.2 Plotting results#

There 3 (three) different options to plot the time response:

  • .plot_1d(): plot time response for given probes.

  • .plot_2d(): plot orbit of a selected node of a rotor system.

  • .plot_3d(): plot orbits for each node on the rotor system in a 3D view.

Ploting time response for list of probes#

probe1 = (3, 0)   # node 3, orientation 0° (X dir.)
probe2 = (3, 90)  # node 3, orientation 90°(Y dir.)

response3.plot_1d(probe=[probe1, probe2], probe_units="degree")

Ploting orbit response for a single node#

node = 26
response3.plot_2d(node=node)

Ploting orbit response for all nodes#

response3.plot_3d()

1.7 Undamped Critical Speed Map (UCS)#

This method will plot the undamped critical speed map for a given range of stiffness values. If the range is not provided, the bearing stiffness at rated speed will be used to create a range.

Whether a synchronous analysis is desired, the method selects only the foward modes and the frequency of the first forward mode will be equal to the speed.

To run the UCS Map, use the command .plot_ucs().

1.7.1 Running and plotting UCS Map#

In this example the UCS Map is calculated for a stiffness range from 10E6 to 10E11 N/m. The other options are left to default.

stiff_range = (6, 11)
ucs_results = rotor3.run_ucs(stiffness_range=stiff_range, num=20, num_modes=16)
ucs_fig = ucs_results.plot()
ucs_fig