Tutorial - MultiRotor System

Tutorial - MultiRotor System#

======================

This is a basic tutorial on how to use the MultiRotor class for rotordynamics analysis. Before starting this tutorial, be sure you’re already familiar with ROSS library. What changes here is basically the part of building the model, since running the analyses will be practically the same as seen in the other tutorials.

When two shafts are joined together by gears, coupling can occur between the lateral and torsional vibration. This interaction is modeled in ROSS based on [Rao et al., 1998] and [Kaplan et al., 2013]. In this work, a typical spur gear pair is modeled as a pair of rigid disks connected by a spring and a damping, considering the pressure angle (\(\alpha\)) and the orientation angle (\(\varphi\)) as shown in the following:

Gear mesh
Figure 1: Global coordinate system of a spur gear pair (Yang et al., 2016).

As you can see, an element not yet shown is needed to model the multi-rotor system, the GearElement or GearElementTVMS. This new element resembles the DiskElement with added attributes, which will be shown next.

The gear mesh stiffness can be constant or varying on time, according to {cite} Ma2014. The GearElement is associated with the constant mesh stiffness and the GearElementTVMS must be created for the time varying mesh stiffness.

Section 1: Constant Mesh Stiffness#

1.1 GearElement Class#

ROSS offers two ways to create a gear element:

  1. Inputting mass and inertia data (GearElement)

  2. Inputting geometrical and material data (GearElementTVMS)

The class GearElement allows you to create gear elements. It is a subclass of DiskElement that requires information related to its pitch diameter and pressure angle. To use constant gear mesh stiffness it is necessary to use the GearElementclass.

import ross as rs
import numpy as np
import pandas as pd

from ross.units import Q_

# Make sure the default renderer is set to 'notebook' for inline plots in Jupyter
import plotly.io as pio
import plotly.graph_objects as go

pio.renderers.default = "notebook"

gear = rs.GearElement(
    n=0,
    m=5,
    Id=0.002,
    Ip=0.004,
    n_teeth=20,
    pitch_diameter=0.5,
    pr_angle=Q_(22.5, "deg"),
    tag="Gear",
)
gear

GearElement(Id=0.002, Ip=0.004, m=5.0, color='Goldenrod', n=0, scale_factor=1.0, tag='Gear')

1.2 MultiRotor Class#

MultiRotor is a subclass of Rotor class. It takes two rotors (driving and driven) as arguments and couple them with their gears. The object created has several methods that can be used to evaluate the dynamics of the model (they all start with the prefix .run_).

To use this class, you must input the already instantiated rotors and each one need at least one gear element.

The shaft elements are renumbered starting with the elements of the driving rotor.

To assemble the matrices, the driving and driven matrices are joined. For the stiffness matrix, the coupling is considered at the nodes of the gears in contact.

1.2.1 Creating a spur geared two-shaft rotor system#

Let’s create a simple model with two rotors connected by a pair of spur gears and supported by flexible bearings, as shown in Fig. 2. For more details on the description of the model, see the work of [Rao et al., 1998].

MultiRotor
Figure 2: Rotors connected by a pair of spur gears (Friswell et al., 2010).

In Figure 2, the first node is 1 (one), but we must remember that in ROSS the node count starts at 0 (zero).

1.2.1.1 Creating material#

# Creating material
material = rs.Material(name="mat_steel", rho=7800, E=207e9, G_s=79.5e9)

1.2.1.2 Creating the driving rotor#

# Rotor 1
L1 = [0.1, 4.24, 1.16, 0.3]
d1 = [0.3, 0.3, 0.22, 0.22]
shaft1 = [
    rs.ShaftElement(
        L=L1[i],
        idl=0.0,
        odl=d1[i],
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L1))
]
generator = rs.DiskElement(
    n=1,
    m=525.7,
    Id=16.1,
    Ip=32.2,
)
disk = rs.DiskElement(
    n=2,
    m=116.04,
    Id=3.115,
    Ip=6.23,
)

gear1 = rs.GearElement(
    n=4,
    m=726.4,
    Id=56.95,
    Ip=113.9,
    n_teeth=328,
    base_diameter=0.5086 * 2,
    pr_angle=Q_(22.5, "deg"),
    helix_angle=0,
)

bearing1 = rs.BearingElement(n=0, kxx=183.9e6, kyy=200.4e6, cxx=3e3)
bearing2 = rs.BearingElement(n=3, kxx=183.9e6, kyy=200.4e6, cxx=3e3)

rotor1 = rs.Rotor(
    shaft1,
    [generator, disk, gear1],
    [bearing1, bearing2],
)

rotor1.plot_rotor()

1.2.1.3 Creating the driven rotor#

# Rotor 2

L2 = [0.3, 5, 0.1]
d2 = [0.15, 0.15, 0.15]
shaft2 = [
    rs.ShaftElement(
        L=L2[i],
        idl=0.0,
        odl=d2[i],
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L2))
]

base_radius = rs.Q_(0.03567, "m")
pressure_angle = rs.Q_(22.5, "deg")
pitch_diameter = 2 * base_radius / np.cos(pressure_angle)

gear2 = rs.GearElement(
    n=0,
    m=5,
    Id=0.002,
    Ip=0.004,
    n_teeth=23,
    pitch_diameter=pitch_diameter,
    pr_angle=pressure_angle,
    helix_angle=0,
)

turbine = rs.DiskElement(n=2, m=7.45, Id=0.0745, Ip=0.149)

bearing3 = rs.BearingElement(n=1, kxx=10.1e6, kyy=41.6e6, cxx=3e3)
bearing4 = rs.BearingElement(n=3, kxx=10.1e6, kyy=41.6e6, cxx=3e3)

rotor2 = rs.Rotor(
    shaft2,
    [gear2, turbine],
    [bearing3, bearing4],
)

rotor2.plot_rotor()

1.2.1.4 Connecting rotors#

To build the multi-rotor model, we need to inform, in the following order:

  • the driving rotor,

  • the driven rotor,

  • the tuple with the pair of coupled nodes (first number corresponds to the gear node of the driving rotor, and the second of the driven rotor),

  • the gear ratio, and

  • the gear mesh stiffness.

Finally, we can inform:

  • the orientation angle (if not defined, zero is adopted as the default),

  • the position of the driven rotor in relation to the driving rotor only for visualization in the plot (“above” or “below”), and

  • a tag.

# Creating multi-rotor
multirotor = rs.MultiRotor(
    rotor1,
    rotor2,
    coupled_nodes=(4, 0),
    gear_mesh_stiffness=1e8,
    orientation_angle=0,
    position="below",
)

multirotor.plot_rotor()

1.2.2 Running analyses#

We will run some analyses for the multi-rotor in this section and even compare results from the literature.

1.2.2.2 Campbell diagram#

To obtain the Campbell diagram we can proceed in the same way as seen for a single rotor. Remember that the reference speeds / frequencies are in relation to the driving rotor.

In the Campbell diagram below, the dashed lines show the shaft rotation speeds corresponding to the generator (blue, node 1) and turbine (yellow, node 7).

frequency_range = Q_(np.arange(0, 5000, 100), "RPM")

gear_ratio = multirotor.mesh.gear_ratio

campbell = multirotor.run_campbell(frequency_range, frequencies=13)
campbell.plot(frequency_units="Hz", harmonics=[1, round(gear_ratio, 3)]).show()

nodes = [2, 7]
unb_mag = [35.505e-3, 0.449e-3]
unb_phase = [0, 0]

dt = 1e-3
t = np.arange(0, 1200, dt)
speed1 = Q_(5000, "RPM").to_base_units().m  # Generator rotor speed

# Unbalance force
F = multirotor.unbalance_force_over_time(nodes, unb_mag, unb_phase, speed1, t)

# Time response
time_resp = multirotor.run_time_response(speed1, F.T, t)
amp_resp = time_resp.yout

Due to computational cost limitations, the entire analysis performed for comparison is not included in this tutorial. However, by running the run_time_response function across the complete range of speeds, it is possible to obtain the following responses at nodes 2 and 7 as shown in the figures below:

MultiRotor MultiRotor
Figure 3: Comparison of ROSS results with Yang et al. (2016) for the spur geared two-shaft rotor system.

1.2.3 Creating a spur geared multi-shaft rotor system#

Let’s create a more complex model with three connected rotors. More details of this example can be found in [Yang et al., 2016].

# Rotor 1
L1 = [300, 92, 200, 200, 92, 300]
r1 = [61.5, 75, 75, 75, 75, 61.5]
shaft1_coupled = [
    rs.ShaftElement(
        L=L1[i] * 1e-3,
        idl=0.0,
        odl=r1[i] * 2e-3,
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L1))
]

D1 = rs.DiskElement(
    n=0,
    m=66.63,
    Id=0.431,
    Ip=0.735,
)
D2 = rs.DiskElement(
    n=6,
    m=69.83,
    Id=0.542,
    Ip=0.884,
)

cxx = 3e3
B1 = rs.BearingElement(n=2, kxx=5.5e8, kyy=6.7e8, cxx=cxx)
B2 = rs.BearingElement(n=4, kxx=5.5e8, kyy=6.7e8, cxx=cxx)

pressure_angle = Q_(22.5, "deg")

G1 = rs.GearElement(
    n=3,
    m=14.37,
    Id=0.068,
    Ip=0.136,
    n_teeth=37,
    base_diameter=0.19,
    pr_angle=pressure_angle,
)

rotor1 = rs.Rotor(
    shaft1_coupled,
    [D1, D2, G1],
    [B1, B2],
)

# Rotor 2
L2 = [80, 200, 200, 640]
r2 = [160.5, 160.5, 130.5, 130.5]
shaft2 = [
    rs.ShaftElement(
        L=L2[i] * 1e-3,
        idl=0.0,
        odl=r2[i] * 2e-3,
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L2))
]

B3 = rs.BearingElement(n=1, kxx=3.2e9, kyy=4.6e9, cxx=cxx)
B4 = rs.BearingElement(n=3, kxx=3.2e9, kyy=4.6e9, cxx=cxx)

G2 = rs.GearElement(
    n=2,
    m=813.79,
    Id=52.36,
    Ip=104.72,
    n_teeth=244,
    base_diameter=1.23,
    pr_angle=pressure_angle,
)

rotor2 = rs.Rotor(
    shaft2,
    [G2],
    [B3, B4],
)

# Rotor 3
L3 = [300, 110, 200, 200, 110, 300]
r3 = [64.5, 76.5, 76.5, 76.5, 76.5, 64.5]
shaft3 = [
    rs.ShaftElement(
        L=L3[i] * 1e-3,
        idl=0.0,
        odl=r3[i] * 2e-3,
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L1))
]

D3 = rs.DiskElement(
    n=0,
    m=95.06,
    Id=1.097,
    Ip=1.532,
)
D4 = rs.DiskElement(
    n=6,
    m=96.22,
    Id=1.110,
    Ip=1.620,
)

G3 = rs.GearElement(
    n=3,
    m=19.52,
    Id=0.098,
    Ip=0.195,
    n_teeth=47,
    base_diameter=0.24,
    pr_angle=pressure_angle,
)

B5 = rs.BearingElement(n=2, kxx=7.2e8, kyy=8.4e8, cxx=cxx)
B6 = rs.BearingElement(n=4, kxx=7.2e8, kyy=8.4e8, cxx=cxx)

rotor3 = rs.Rotor(
    shaft3,
    [D3, D4, G3],
    [B5, B6],
)

# Connect rotor 1 with rotor 2 (driving rotor)
multi_rotor1 = rs.MultiRotor(
    rotor2,
    rotor1,
    coupled_nodes=(2, 3),
    gear_mesh_stiffness=2.55e8,
    orientation_angle=Q_(270, "deg"),
    position="above",
)

# Connect rotor 3 with rotor 2 in multi rotor
psi = 90
final_system = rs.MultiRotor(
    multi_rotor1,
    rotor3,
    coupled_nodes=(2, 3),
    gear_mesh_stiffness=2.6e8,
    orientation_angle=Q_(270 - psi, "deg"),
    position="below",
)

final_system.plot_rotor()

If we apply two unbalances with 92 g.mm in phase opposition at two ends of rotor 1 (nodes 5 and 11), while one unbalance with 294 g.mm at the middle of rotor 3 (node 15), we can obtain the following responses at node 15:

MultiRotor
Figure 4: Comparison of ROSS results with Yang et al. (2016) for the spur geared multi-shaft rotor system.

We can also vary the orientation angle between the rotors:

MultiRotor
Figure 5: Unbalance responses at node 15 with three orientation angles.

Section 2: Time Variant Mesh Stiffness#

The method for Time varying Mesh Stiffness is used when it is necessary to consider the gear mesh stiffness variation across the time. This behavior comes when considering the complete tooth profile and the contact ratio of the gear mesh.

2.1 GearElementTVMS Class#

The class GearElementTVMS allows you to create gear elements which includes the calculation for the time varying mesh stiffness.

To create a GearElementTVMS object some design parameters are needed.

# Rebuilding Rotor 1

L1 = [0.1, 4.24, 1.16, 0.3]
d1 = [0.3, 0.3, 0.22, 0.22]
shaft1 = [
    rs.ShaftElement(
        L=L1[i],
        idl=0.0,
        odl=d1[i],
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L1))
]
generator = rs.DiskElement(
    n=1,
    m=525.7,
    Id=16.1,
    Ip=32.2,
)
disk = rs.DiskElement(
    n=2,
    m=116.04,
    Id=3.115,
    Ip=6.23,
)

bearing1 = rs.BearingElement(n=0, kxx=183.9e6, kyy=200.4e6, cxx=3e3)
bearing2 = rs.BearingElement(n=3, kxx=183.9e6, kyy=200.4e6, cxx=3e3)

# Building Gear 1 as GearElement TVMS


# Common Gear Mesh Variables

steel = rs.Material(name="Steel", rho=Q_(7850, "kg/m**3"), E=2e11, Poisson=0.3)
helix_angle = Q_(0, "degree")
pressure_angle = Q_(22.5, "degree")


# Calculation for the particular gear variables
# this calculations were done only for adapt the paper variables to the ones needed in the GearElementTVMS

bore_diameter = Q_(0.22, "m")  # rotor outer diameter = internal diameter of the gear
m = Q_(726.4, "kg")  # mass
Ip = Q_(113.9, "kg*m**2")  # Polar Inertia

pitch_diameter = np.sqrt(
    (8 * Ip) / (m) - bore_diameter**2
)  # represents the gear outer diameter

rho = Q_(7850, "kg/m**3")

width = (4 * m) / (np.pi * rho * (pitch_diameter**2 - bore_diameter**2))  # gear width

n_teeth = 328
module = pitch_diameter / n_teeth

node = 4

# Recreating Gear Element with GearElementTVMS
gear1 = rs.GearElementTVMS(
    n=node,
    material=steel,
    width=width,
    bore_diameter=bore_diameter,
    module=module,
    n_teeth=n_teeth,
    pr_angle=pressure_angle,
    helix_angle=helix_angle,
)

# recreating rotor 1

rotor1 = rs.Rotor(
    shaft1,
    [generator, disk, gear1],
    [bearing1, bearing2],
)

rotor1.plot_rotor()

# Rebuilding Rotor 2

L2 = [0.3, 5, 0.1]
d2 = [0.15, 0.15, 0.15]
shaft2 = [
    rs.ShaftElement(
        L=L2[i],
        idl=0.0,
        odl=d2[i],
        material=material,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(len(L2))
]

turbine = rs.DiskElement(n=2, m=7.45, Id=0.0745, Ip=0.149)

bearing3 = rs.BearingElement(n=1, kxx=10.1e6, kyy=41.6e6, cxx=3e3)
bearing4 = rs.BearingElement(n=3, kxx=10.1e6, kyy=41.6e6, cxx=3e3)

# Building Gear 2 as GearElement TVMS

# Common Gear Mesh Variables

steel = rs.Material(name="Steel", rho=7850, E=2e11, Poisson=0.3)
helix_angle = Q_(0, "degree")
pressure_angle = Q_(22.5, "degree")

# Calculation for the particular gear variables
# this calculations were done only for adapt the paper variables to the ones needed in the GearElementTVMS

bore_diameter = Q_(0.015, "m")  # internal diameter of the gear
m = 5  # mass
base_radius = Q_(0.03567, "m")  # gear base radius

pitch_diameter = (
    2 * base_radius / np.cos(pressure_angle)
)  # represents the gear outer diameter

n_teeth = 23
module = pitch_diameter / n_teeth

node = 0

# Recreating Gear Element with GearElementTVMS
gear2 = rs.GearElementTVMS(
    n=node,
    material=steel,
    width=width,
    bore_diameter=bore_diameter,
    module=module,
    n_teeth=n_teeth,
    pr_angle=pressure_angle,
    helix_angle=helix_angle,
)

# recreating rotor 2

rotor2 = rs.Rotor(
    shaft2,
    [gear2, turbine],
    [bearing3, bearing4],
)

rotor2.plot_rotor()

2.2 Gear Mesh Stiffness - Default Profile#

To obtain the default gear mesh stiffness profile it is necessary to use the GearElementTVMS and in the MultiRotor flag update_mesh_stiffness must be TRUE.

multirotor = rs.MultiRotor(
    rotor1,
    rotor2,
    coupled_nodes=(4, 0),
    update_mesh_stiffness=True,
    position="below",
    orientation_angle=0,
)

multirotor.plot_rotor()

C:\Users\Murillo\OneDrive - Universidade Federal de Uberlândia\Área de Trabalho\Mestrado\ENGRENAMENTO\Implementacao\ross_dev_backlash\ross\ross\gear_element.py:707: UserWarning:

Extrapolating gear body coefficients described by Sainsot et. al. (2014). Be careful when post-processing the results.

To plot the stiffness profile just use .mesh.plot_stiffness_profile function presenting the number of mesh periods to plot.

multirotor.mesh.plot_stiffness_profile(n_mesh_period=2)

C:\Users\Murillo\OneDrive - Universidade Federal de Uberlândia\Área de Trabalho\Mestrado\ENGRENAMENTO\Implementacao\ross_dev_backlash\ross\ross\gear_element.py:707: UserWarning:

Extrapolating gear body coefficients described by Sainsot et. al. (2014). Be careful when post-processing the results.

It is possible to see the values of average gear mesh stiffness and gear mesh contact ratio using the functions .mesh.stiffness and .mesh.contact_ratio

stiff_avg = multirotor.mesh.stiffness
print(f"Média de Rigidez = {stiff_avg:.2e} N/m")

contact_ratio = multirotor.mesh.contact_ratio
print(f"Razão de Contato = {contact_ratio:.2f}")

Média de Rigidez = 1.85e+09 N/m
Razão de Contato = 1.64

2.3 Gear Mesh Stiffness - Rectangular Profile#

To obtain the rectangular gear mesh stiffness profile it is necessary to use the GearElementTVMS, the flag update_mesh_stiffness must be TRUE and it is important to use two more argument.

square_varying_stiffness argument must be TRUE and square_stiffness_amplitude_ratiorepresents the amplitude ratio considering the mean value.

multirotor = rs.MultiRotor(
    rotor1,
    rotor2,
    coupled_nodes=(4, 0),
    update_mesh_stiffness=True,
    square_varying_stiffness=True,
    square_stiffness_amplitude_ratio=0.225,
    position="below",
    orientation_angle=0,
)

To plot the stiffness profile just use .mesh.plot_stiffness_profile function presenting the number of mesh periods to plot.

multirotor.mesh.plot_stiffness_profile(n_mesh_period=2)

It is possible to see the values of average gear mesh stiffness and gear mesh contact ratio using the functions .mesh.stiffness and .mesh.contact_ratio

stiff_avg = multirotor.mesh.stiffness
print(f"Média de Rigidez = {stiff_avg:.2e} N/m")

contact_ratio = multirotor.mesh.contact_ratio
print(f"Razão de Contato = {contact_ratio:.2f}")

Média de Rigidez = 1.85e+09 N/m
Razão de Contato = 1.64

2.4 Analysis and multi-shaft rotor system#

By using the class GearElementTVMS to calculte the time varying mesh stiffness there is not any changes in the way to run analyses and build multi-shaft rotor system, as shwon in the section 1.

References#

[Fri10]

Michael I Friswell. Dynamics of rotating machines. Cambridge University Press, 2010.

[KDA+13]

Jason Kaplan, sa Dousti, Paul Allaire, Bradley Nichols, Timothy Dimond, and Alexandrina Untaroiu. Rotor dynamic modeling of gears and geared systems. Proceedings of the ASME Turbo Expo, 7:, 06 2013. doi:10.1115/GT2013-94654.

[RSC98] (1,2)

JS Rao, TN Shiau, and JR Chang. Theoretical analysis of lateral response due to torsional excitation of geared rotors. Mechanism and Machine Theory, 33(6):761–783, 1998. doi:https://doi.org/10.1016/S0094-114X(97)00056-6.

[YWWD16]

Yi Yang, Jiaying Wang, Xurong Wang, and Yiping Dai. A general method to predict unbalance responses of geared rotor systems. Journal of Sound and Vibration, 381:246–263, 2016. doi:https://doi.org/10.1016/j.jsv.2016.06.031.

Contents