Implement co-axial rotors

[raphaeltimbo]

I am opening this issue so we can discuss the implementation of co-axial rotors. Discussion for this started on issue #245 , since it involves linking rotors together. Some points regarding the implementation:

• We will need to discuss the best command interface to the user regarding this. I was thinking on the following options:
  • Implement a new CoAxialRotor class that would inherit from Rotor. To create the co-axial rotor, the user would give the different rotors already assembled and the correspondent links.
  • Change the current Rotor class to receive more than one shaft. To implement this, we would need to have a way to differentiate the shafts (XLTRC2 has a ‘shaft’ number for that).

Right now I am leaning towards the first option. It looks like it will be easier to implement and we won’t need to change code that is already working.

As a reference, Friswell's example 6.6.1 uses a co-axial rotor.

[rodrigomoliveira1]

I'm trying to work with a new class like you said. However, that's a lot of data to adjust to make things work. For example:

If I instantiate 2 rotors:

rotor1 = rotor_example1() rotor2 = rotor_example2()

and I use both as input to CoAxialRotor, then, for at least one of these rotors I must modify the nodes, nodes_pos, dof_global_index and all the others attributes related to their nodes, including all DataFrames (df, df_shaft, df_disks, ...). And we'll need to modify the methods that calculate the rotor matrices.

On the other hand, changing the current Rotor class to receive more than one shaft would be simpler. We could adapt it to receive the shafts and the starting point position for each one. The links between them could be done through a BearingElement with a n_link. Thus, we wouldn't need to change the other methods, once everything else is calculated based on the geometry (nodes and nodes positions, ...).

[rodrigomoliveira1]

I've managed to implement the coaxial rotor code modifying the Rotor class. Take a look at the class structure now:

class Rotor(object):
    r"""A rotor object.

    This class will create a rotor with the shaft,
    disk, bearing and seal elements provided.

    Parameters
    ----------
    axial_shaft : list
        List with the shaft elements (main shaft)
    disk_elements : list
        List with the disk elements
    bearing_seal_elements : list
        List with the bearing elements
    point_mass_elements: list
        List with the point mass elements
    coaxial_shaft : list
        List with the shaft elements of a co-axial shaft.
    shaft_start_pos : list
        List indicating the initial node position for each shaft.
        Default is zero for each shaft created.
    sparse : bool, optional
        If sparse, eigenvalues will be calculated with arpack.
        Default is True.
    n_eigen : int, optional
        Number of eigenvalues calculated by arpack.
        Default is 12.
    tag : str
        A tag for the rotor

    Returns
    -------
    A rotor object.

    Attributes
    ----------
    evalues : array
        Rotor's eigenvalues.
    evectors : array
        Rotor's eigenvectors.
    wn : array
        Rotor's natural frequencies in rad/s.
    wd : array
        Rotor's damped natural frequencies in rad/s.

    Examples
    --------
    >>> #  Rotor without damping with 2 shaft elements 1 disk and 2 bearings
    >>> import ross as rs
    >>> steel = rs.materials.steel
    >>> z = 0
    >>> le = 0.25
    >>> i_d = 0
    >>> o_d = 0.05
    >>> tim0 = rs.ShaftElement(le, i_d, o_d,
    ...                        material=steel,
    ...                        shear_effects=True,
    ...                        rotary_inertia=True,
    ...                        gyroscopic=True)
    >>> tim1 = rs.ShaftElement(le, i_d, o_d,
    ...                        material=steel,
    ...                        shear_effects=True,
    ...                        rotary_inertia=True,
    ...                        gyroscopic=True)
    >>> shaft_elm = [tim0, tim1]
    >>> disk0 = rs.DiskElement.from_geometry(1, steel, 0.07, 0.05, 0.28)
    >>> stf = 1e6
    >>> bearing0 = rs.BearingElement(0, kxx=stf, cxx=0)
    >>> bearing1 = rs.BearingElement(2, kxx=stf, cxx=0)
    >>> rotor = rs.Rotor(shaft_elm, [disk0], [bearing0, bearing1])
    >>> modal = rotor.run_modal(speed=0)
    >>> modal.wd[0] # doctest: +ELLIPSIS
    215.3707...
    """

    def __init__(
        self,
        axial_shaft,
        disk_elements=None,
        bearing_seal_elements=None,
        point_mass_elements=None,
        coaxial_shaft=None,
        shaft_start_pos=[0.0, 0.0],
        sparse=True,
        n_eigen=12,
        min_w=None,
        max_w=None,
        rated_w=None,
        tag=None,
    ):

I did for 2 shafts only. I don't know if there is a system with more co-axial shafts. Anyway, it can be adapted to more than 2 shafts if needed. Basically, I've changed the __init__() to recognize 2 different shafts: an rotor.axial_shaft and a rotor.coaxial_shaft. There's no need to point out the nodes range. The only thing here is that the user must input the axis starting points (shaft_start_pos) to define the shafts' positions.

You can still call for rotor.shaft_elements, which will display all the elements from both shafts. So, now we have

• rotor.shaft_elements - returns the list of shaft elements from both shafts;
• rotor.axial_shaft - returns the list of shaft elements from axial_shaft input;
• rotor.coaxial_shaft - returns the list of shaft elements from coaxial_shaft input;

Here is an example:

def rotor_example():
    i_d = 0
    o_d = 0.05
    n = 10
    L = [0.25 for _ in range(n)]

    axial_shaft = [
        ShaftElement(
            l, i_d, o_d, material=steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
        )
        for l in L
    ]

    i_d = 0.15
    o_d = 0.20
    n = 6
    L = [0.25 for _ in range(n)]

    coaxial_shaft = [
        ShaftElement(
            l, i_d, o_d, material=steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
        )
        for l in L
    ]

    disk0 = DiskElement.from_geometry(
        n=1, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk1 = DiskElement.from_geometry(
        n=9, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk2 = DiskElement.from_geometry(
        n=13, material=steel, width=0.07, i_d=0.20, o_d=0.48
    )
    disk3 = DiskElement.from_geometry(
        n=15, material=steel, width=0.07, i_d=0.20, o_d=0.48
    )

    stfx = 1e6
    stfy = 0.8e6
    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=0)
    bearing1 = BearingElement(10, kxx=stfx, kyy=stfy, cxx=0)
    bearing2 = BearingElement(11, kxx=stfx, kyy=stfy, cxx=0)
    bearing3 = BearingElement(8, n_link=17, kxx=stfx, kyy=stfy, cxx=0)

    return Rotor(axial_shaft, [disk0, disk1, disk2, disk3], [bearing0, bearing1, bearing2, bearing3],
                 coaxial_shaft=coaxial_shaft, shaft_start_pos=[0.0, 0.5])

>>> rotor = rotor_example()
>>> show(rotor.plot_rotor())

Of course the representation of the bearing that connects the 2 shafts are not good. We have to think in another way to draw this specific bearing. Another problem is the node numbers representation. We need to change how it's displayed to avoid overlapping.

[raphaeltimbo]

As discussed, I believe we should create a separate class for co-axial rotors to avoid having to many arguments on our main rotor class.

[rodrigomoliveira1]

@raphaeltimbo, i've changed what i've done before to create a class for co-axial rotors. I tried something more general allowing the user to add as many shafts he/she wants.

The instantiation is very similar to class Rotor, by the way. There're 2 main differences when using CoAxialRotor.

• The user must input a list with lists of ShaftElements, so the code can differentiate a shaft from another. This way we don't need to identify every element to each shaft.

• The user must input the shafts' positions, where the fisrt node of each shaft is located. This new attribute is called shaft_start_pos. It's basically a list.

However, I had to recreate the .__init__() method to CoAxialRotor.

This is an example of how to instantiate CoAxialRotor class:

    i_d = 0
    o_d = 0.05
    n = 10
    L = [0.25 for _ in range(n)]

    axial_shaft = [ShaftElement(l, i_d, o_d, material=steel) for l in L]

    i_d = 0.15
    o_d = 0.20
    n = 6
    L = [0.25 for _ in range(n)]

    coaxial_shaft = [ShaftElement(l, i_d, o_d, material=steel) for l in L]

    shaft = [axial_shaft, coaxial_shaft]

    disk0 = DiskElement.from_geometry(
        n=1, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk1 = DiskElement.from_geometry(
        n=9, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk2 = DiskElement.from_geometry(
        n=13, material=steel, width=0.07, i_d=0.20, o_d=0.48
    )
    disk3 = DiskElement.from_geometry(
        n=15, material=steel, width=0.07, i_d=0.20, o_d=0.48
    )

    stfx = 1e6
    stfy = 0.8e6
    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=0)
    bearing1 = BearingElement(10, kxx=stfx, kyy=stfy, cxx=0)
    bearing2 = BearingElement(11, kxx=stfx, kyy=stfy, cxx=0)
    bearing3 = BearingElement(7, n_link=16, kxx=stfx, kyy=stfy, cxx=0)

    return CoAxialRotor(shaft, [disk0, disk1, disk2, disk3],
                         [bearing0, bearing1, bearing2, bearing3],
                         shaft_start_pos=[0, 0.5])

[rodrigomoliveira1]

Now i have to work on the coupling between properties of nodes of different shafts. The way I see, the most automatic way to do it is checking the n_link from bearings. If a bearing n_link corresponds to another shaft node, then, the mass, stiffness, etc, properties from one shaft node (where the bearing is located - attribute n) will couple with the other shaft node (where the bearing links - attribute n_link)

[rodrigomoliveira1]

I didn't open a PR for this yet because I found a bug while plotting co-axial rotors with bearings in the same axial position.

The bearings should be drawn at node 2 and 11 respectively. But the second bearing appears at the top of the first bearing, like they were linked.

[rodrigomoliveira1]

The problem is solved. Now bearing icons are plotted respecting the diameters from both shafts they connect.

(generic example below)

[raphaeltimbo]

Nice! When you have the tests + PR completed we can merge this.

[raphaeltimbo]

@rodrigomoliveira1 , check if your code is working for the following example:

import ross as rs
from bokeh.plotting import show
steel = rs.materials.steel

def rotor_example():
    i_d = 0
    o_d = 0.05
    n = 6
    L = [0.25 for _ in range(n)]

    shaft_elem = [
        rs.ShaftElement(
            l, i_d, o_d, material=steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
        )
        for l in L
    ]

    disk0 = rs.DiskElement.from_geometry(
        n=2, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk1 = rs.DiskElement.from_geometry(
        n=4, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )

    stfx = 1e6
    stfy = 0.8e6
    bearing0 = rs.BearingElement(0, n_link=7, kxx=stfx, kyy=stfy, cxx=0, scale_factor=1)
    support0 = rs.BearingElement(7, kxx=stfx, kyy=stfy, cxx=0, tag="Support", scale_factor=1)
    bearing1 = rs.BearingElement(6, n_link=8, kxx=stfx, kyy=stfy, cxx=0, scale_factor=1)
    support1 = rs.BearingElement(8, kxx=stfx, kyy=stfy, cxx=0, tag="Support", scale_factor=1)

    point_mass0 = rs.PointMass(7, m=1.0)
    point_mass1 = rs.PointMass(8, m=1.0)

    return rs.Rotor(shaft_elem, [disk0, disk1], [bearing0, support0, support1, bearing1],
                 [point_mass0, point_mass1])

rotor = rotor_example()
show(rotor.plot_rotor())

For now I am getting the following plot:

[rodrigomoliveira1]

It's odd. I tried the same code than you, but did a minor change:

import ross as rs
from bokeh.plotting import show
steel = rs.materials.steel

def rotor_example():
    i_d = 0
    o_d = 0.05
    n = 6
    L = [0.25 for _ in range(n)]

    shaft_elem = [
        rs.ShaftElement(
            l, i_d, o_d, material=steel, shear_effects=True, rotary_inertia=True, gyroscopic=True
        )
        for l in L
    ]

    disk0 = rs.DiskElement.from_geometry(
        n=2, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )
    disk1 = rs.DiskElement.from_geometry(
        n=4, material=steel, width=0.07, i_d=0.05, o_d=0.28
    )

    stfx = 1e6
    stfy = 0.8e6
    bearing0 = rs.BearingElement(0, n_link=7, kxx=stfx, kyy=stfy, cxx=0, scale_factor=1)
    support0 = rs.BearingElement(7, kxx=stfx, kyy=stfy, cxx=0, tag="Support", scale_factor=1)
    bearing1 = rs.BearingElement(6, n_link=8, kxx=stfx, kyy=stfy, cxx=0, scale_factor=1)
    support1 = rs.BearingElement(8, kxx=stfx, kyy=stfy, cxx=0, tag="Support", scale_factor=1)

    point_mass0 = rs.PointMass(7, m=1.0)
    point_mass1 = rs.PointMass(8, m=1.0)

    # here I did [bearing0, bearing1, support0, support1] instead of 
    # [bearing0, support0, support1, bearing1]
    return rs.Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1, support0, support1],
                 [point_mass0, point_mass1])

rotor = rotor_example()
show(rotor.plot_rotor()) 

and got this:

I'm investigatig why the arguments position causes this issue

[rodrigomoliveira1]

I've solved this problem.

I've just had to sort the bearing_seal_elements by the element node, just like it's done to shaft_elements And while cycling through the bearing without z_pos, I replaced bearing_seal_elements by self.bearing_seal_elements

self.bearing_seal_elements = sorted(bearing_seal_elements, key=lambda el: el.n)

I'm just waiting PRs #432 and #433 to be merged and I upload the modifications. If I open anything new now, it's getting a conflict.

[rodrigomoliveira1]

While it is still unviable to create a PR to add this class, I working on adapting the other methods to the CoAxialRotor.

I thought it would be simpler to do it. I'm trying to make it possible to plot a result for each shaft. But since the results are stored in simple lists, the data for each shaft are mixed, which makes the visualization of the graphs confusing and the results presented wrong.

So, I decided to open a new column on the pd.DataFrame (df["shaft_number"]). It defines which shaft each element belongs to. However some BearingElements would belong to more the one shaft, due the connections between them, and this can be a little tricky to the DataFrame.

Look the example below:

For now df["shaft_number"] returns a Series of floats. (check the last column below)

>>> rotor.df[["type", "n_l", "n_link", "shaft_number"]]

              type   n_l  n_link  shaft_number
0     ShaftElement   0.0     NaN           0.0
1   BearingElement   0.0    18.0           0.0
2     ShaftElement   1.0     NaN           0.0
3      DiskElement   1.0     NaN           0.0
4     ShaftElement   2.0     NaN           0.0
5     ShaftElement   3.0     NaN           0.0
6     ShaftElement   4.0     NaN           0.0
7     ShaftElement   5.0     NaN           0.0
8     ShaftElement   6.0     NaN           0.0
9     ShaftElement   7.0     NaN           0.0
10  BearingElement   8.0    17.0           0.0
11    ShaftElement   8.0     NaN           0.0
12    ShaftElement   9.0     NaN           0.0
13     DiskElement   9.0     NaN           0.0
14  BearingElement  10.0     NaN           0.0
15  BearingElement  11.0    19.0           1.0
16    ShaftElement  11.0     NaN           1.0
17    ShaftElement  12.0     NaN           1.0
18    ShaftElement  13.0     NaN           1.0
19     DiskElement  13.0     NaN           1.0
20    ShaftElement  14.0     NaN           1.0
21     DiskElement  15.0     NaN           1.0
22    ShaftElement  15.0     NaN           1.0
23    ShaftElement  16.0     NaN           1.0
24  BearingElement  18.0     NaN           0.0
25  BearingElement  19.0     NaN           1.0
26       PointMass   NaN     NaN           0.0
27       PointMass   NaN     NaN           1.0

Notice that the BearingElement in the position 10 belongs to the shaft number 0 because BearingElement.n == 8 (node 8 belongs to the inner shaft). But it also should belong to the shaft number 1, because BearingElement.n_link == 17 (node 17 belongs to the outer shaft).

I tried making it a list, so in the index 10 of the DataFrame we would get this:

10 BearingElement 8.0 17.0 list(0.0, 1.0)

But this way, it's difficult to iterate over this column because we would have something like: [0.0, ... , 0.0, list(0.0, 1.0), 1.0, ... , 1.0]

If you guys have any idea to overcome this situation, I'm open.

[raphaeltimbo]

I thought of maybe having two different rotor objects within a coaxial rotor object, this way we could have all the methods that we already have for each rotor, and we would also need some code to link the rotor when needed (getting global matrices for example).

[rodrigomoliveira1]

I thought of maybe having two different rotor objects within a coaxial rotor object, this way we could have all the methods that we already have for each rotor, and we would also need some code to link the rotor when needed (getting global matrices for example).

So, I'm thinking of the instances for the class.

• a list of rotors e.g: [rotor1, rotor2]

• a list of nodes to connect the rotors e.g: [[rotor1.nodes[2], rotor2.nodes[3]], [rotor1.nodes[5], rotor2.nodes[6]]]

• a list to indicating the initial position of each rotor e.g: [0.0, 0.5]

Maybe, we have to renumerate the nodes, for following rotors, to build the global matrices and set the DoF's in the right location.

class CoAxialRotor(Rotor):
    r"""A co-axial rotor object.

    This class will create a system of co-axial rotors with a system of rotors provided.

    Parameters
    ----------
    rotors: list
        List of rotors
    nodes_connection: list
        List indicating the rotor nodes to be connected
    shaft_start_pos : list
        List indicating the initial node position for each rotor.
    tag : str
        A tag for the co-axial rotor
    """

    def __init__(
        self,
        rotors,
        nodes_connection,
        shaft_start_pos,
        tag=None,
    ):

[raphaeltimbo]

class CoAxialRotor(Rotor):
…
    nodes_connection: list
        List indicating the rotor nodes to be connected
    shaft_start_pos : list
        List indicating the initial node position for each rotor.
…

Are these two parameters needed? If the user is providing bearings with their corresponding links, these would already have the “nodes_connection” and from that I think we could also get the “shaft_start_pos” (although it might be complicated to do it). Here the user would have to provide shaft lists with nodes already defined.

[rodrigomoliveira1]

Actually, I just removed shaft_start_pos. I was trying to do like you said before but wasn't successful. But I did it now.

It'll work better for a single connection. For multiple connections the user has to make sure the nodal positions match (just by plotting it's possible to check), or the code will place the coaxial rotor based on the bearing with the greatest n value

If the user is providing bearings with their corresponding links, these would already have the “nodes_connection”

I agree, but for multiple connections (e.g. 2 rotors with 2 n_link each), maybe the code would need more information. Or we could assume an increasing order (based on bearing.n) to link the rotors.

[rodrigomoliveira1]

I would go with the option two for the coaxial object. If the user wants a separated plot for each rotor he should be able to get that by accessing each rotor. Something like: results = coaxial_rotor.rotors[0].run_static()

@raphaeltimbo, I've thinking of what you said in the Issue #443, and I don't think it's going to work. Running a simulation for each rotor separately would result in loss of information, because the rotors would no longer be coupled. Besides, a single rotor would have mismatched matrices shapes due the bearing.n_link (.K() would return a matrix larger than .M())

In my opinion, we could instantiate the CoAxialRotor class with rotors object, but without the option to run them individually