Links are not working for more than one support

[raphaeltimbo]

The current test where we have links has only one support. We should add another test with supports on both bearings to evaluate if the implementation is correct.

[rodrigomoliveira1]

I have tried this simple example below and it doesn't work:

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

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

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

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

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

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

rotor = rotor_example()
rotor.M().shape

Traceback (most recent call last):

  File "<ipython-input-41-d8ae27a43e7d>", line 1, in <module>
    rotor.M().shape

  File "C:/Users/USER/Google Drive/Backup notebook/Lavi/ross-master/ross/rotor_assembly.py", line 678, in M
    M0[np.ix_(dofs, dofs)] += elm.M()

IndexError: index 32 is out of bounds for axis 0 with size 32

The size of the mass matrix is correct (32 x 32), however when we come to observe the elements' indexes we can notice a bug on the last line (shown in the print below). The Global Index should be GlobalIndex(x_8=30, y_8=31)

GlobalIndex(x_0=0, y_0=1, alpha_0=2, beta_0=3, x_1=4, y_1=5, alpha_1=6, beta_1=7)
GlobalIndex(x_1=4, y_1=5, alpha_1=6, beta_1=7, x_2=8, y_2=9, alpha_2=10, beta_2=11)
GlobalIndex(x_2=8, y_2=9, alpha_2=10, beta_2=11, x_3=12, y_3=13, alpha_3=14, beta_3=15)
GlobalIndex(x_3=12, y_3=13, alpha_3=14, beta_3=15, x_4=16, y_4=17, alpha_4=18, beta_4=19)
GlobalIndex(x_4=16, y_4=17, alpha_4=18, beta_4=19, x_5=20, y_5=21, alpha_5=22, beta_5=23)
GlobalIndex(x_5=20, y_5=21, alpha_5=22, beta_5=23, x_6=24, y_6=25, alpha_6=26, beta_6=27)
GlobalIndex(x_2=8, y_2=9, alpha_2=10, beta_2=11)
GlobalIndex(x_4=16, y_4=17, alpha_4=18, beta_4=19)
GlobalIndex(x_0=0, y_0=1, x_7=28, y_7=29)
GlobalIndex(x_7=28, y_7=29)
GlobalIndex(x_8=32, y_8=33)

[rodrigomoliveira1]

I'm debugging some methods like .dof_global_index() and .dof_local_index() from class Element and I've noticed something strange.

I added a simple print in both methods, just to check whenever they're called and the code starting running their routines respectively. Look what I did:

    def dof_local_index(self):

        print("running dof_local_index")

        dof_mapping = self.dof_mapping()
        dof_tuple = namedtuple("LocalIndex", dof_mapping)
        local_index = dof_tuple(**dof_mapping)

        return local_index

and

    def dof_global_index(self):

        print("running dof_global_index")

       dof_mapping = self.dof_mapping()
        global_dof_mapping = {}
        for k, v in dof_mapping.items():
            dof_letter, dof_number = k.split("_")
            global_dof_mapping[dof_letter + "_" + str(int(dof_number) + self.n)] = v
        dof_tuple = namedtuple("GlobalIndex", global_dof_mapping)

        for k, v in global_dof_mapping.items():
            global_dof_mapping[k] = 4 * self.n + v

        global_index = dof_tuple(**global_dof_mapping)

        return global_index

When I run an example for a single shaft element, calling both method this is what happens:

>>> shaft = ShaftElement(1, 0, 1, steel, 2)
>>> shaft.dof_global_index()
running dof_global_index
Out[1]: GlobalIndex(x_2=8, y_2=9, alpha_2=10, beta_2=11, x_3=12, y_3=13, alpha_3=14, beta_3=15)

>>> shaft.dof_local_index()
running dof_local_index
Out[2]: LocalIndex(x_0=0, y_0=1, alpha_0=2, beta_0=3, x_1=4, y_1=5, alpha_1=6, beta_1=7)

So, we can see that everything is running fine, as expected.

But now, let's run and example for point mass element (the same behaviour occurs for disk, bearings and seals):

>>> pointmass = point_mass_example()
>>> pointmass.dof_global_index()
Out[3]: GlobalIndex(x_0=0, y_0=1)

>>> pointmass.dof_local_index()
Out[4]: LocalIndex(x_0=0, y_0=1)

Although the outputs are correct, the print tests are not being printed when running elements but shaft elements.

let's see another example to get a wider view: Using the .rotor_example(), we'll return the rotor's mass matrix, checking every element degree of freedom (notice the two prints added inside .M())

    def M(self):

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

        for elm in self.elements:
            print("element type: ", elm.type)
            dofs = elm.dof_global_index()
            print("dof ", dofs)
            M0[np.ix_(dofs, dofs)] += elm.M()

        return M0

>>> rotor1 = rotor_example()
>>> rotor1.M().shape
element type:  ShaftElement
running dof_global_index
dof  GlobalIndex(x_0=0, y_0=1, alpha_0=2, beta_0=3, x_1=4, y_1=5, alpha_1=6, beta_1=7)
element type:  ShaftElement
running dof_global_index
dof  GlobalIndex(x_1=4, y_1=5, alpha_1=6, beta_1=7, x_2=8, y_2=9, alpha_2=10, beta_2=11)
element type:  ShaftElement
running dof_global_index
dof  GlobalIndex(x_2=8, y_2=9, alpha_2=10, beta_2=11, x_3=12, y_3=13, alpha_3=14, beta_3=15)
element type:  ShaftElement
running dof_global_index
dof  GlobalIndex(x_3=12, y_3=13, alpha_3=14, beta_3=15, x_4=16, y_4=17, alpha_4=18, beta_4=19)
element type:  ShaftElement
running dof_global_index
dof  GlobalIndex(x_4=16, y_4=17, alpha_4=18, beta_4=19, x_5=20, y_5=21, alpha_5=22, beta_5=23)
element type:  ShaftElement
running dof_global_index
dof  GlobalIndex(x_5=20, y_5=21, alpha_5=22, beta_5=23, x_6=24, y_6=25, alpha_6=26, beta_6=27)
element type:  DiskElement
dof  GlobalIndex(x_2=8, y_2=9, alpha_2=10, beta_2=11)
element type:  DiskElement
dof  GlobalIndex(x_4=16, y_4=17, alpha_4=18, beta_4=19)
element type:  BearingElement
dof  GlobalIndex(x_0=0, y_0=1)
element type:  BearingElement
dof  GlobalIndex(x_6=24, y_6=25)
Out[5]: (28, 28)

As we can see, when the loop reaches elements different than shaft elements, the prints do not return. I don't know if this kind of behaviour was supposed to happen, but this hinders the changes to fix the link bug.

[rodrigomoliveira1]

Nevermind, I've found the problem. Some elements were importing .element.py from another file.

[rodrigomoliveira1]

I found a solution for this problem, but I don't know if it's the best way to solve it.

Here, we have the attribute that indicates de number of degrees of freedom.

        # number of dofs
        self.ndof = (
            4 * max([el.n for el in shaft_elements])
            + 8
            + 2 * len([el for el in point_mass_elements])
        )

This equation won't match the .dof_global_index() method to return the global dof's for each element when instantiating more than one support, for example. Then, the code will break at some point when running the mass / stiffness / damping matrices.

We could change it to + 4 * len([el for el in point_mass_elements]) and turn the PointMass mass matrix a 4x4, considering inertias equal zero. This creates another problem: we would have singular matrices when calculating global matrices. So I solved this by adding these two lines to .M(), .C(), .G() and .K() methods in rotor_assembly.py:

    def M(self):
        M0 = np.zeros((self.ndof, self.ndof))

        for elm in self.elements:e)
            dofs = elm.dof_global_index()
            M0[np.ix_(dofs, dofs)] += elm.M()

        M0 = M0[~np.all(M0 == 0, axis=1)]
        M0 = M0[:, ~np.all(M0 == 0, axis=0)]

        return M0

These 2 lines remove rows and columns filled with zeros in a given matrix avoiding, then, getting a singular matrix.

        M0 = M0[~np.all(M0 == 0, axis=1)]
        M0 = M0[:, ~np.all(M0 == 0, axis=0)]

The problem is: the attribute ndof won't match the actual matrices shape. We could assign a new value for ndof when running the matrices methods. @raphaeltimbo, @GabrielBachmann any ideas?

[raphaeltimbo]

I believe we have a problem on how we are calculating the global index for the PointMass. For example:

>>> for p in (rotor.point_mass_elements):
>>>    print(f'{p.tag}: {p.dof_global_index()}')
Point Mass 0: GlobalIndex(x_7=28, y_7=29)
Point Mass 1: GlobalIndex(x_8=32, y_8=33)

I think the correct result here should be: Point Mass 1: GlobalIndex(x_8=30, y_8=31)

[rodrigomoliveira1]

Totally agree with you, because the .dof_global_index() method always multiplies the elements' n value by 4. But for point masses, it should be multiplied by 2.

        for k, v in global_dof_mapping.items():
            global_dof_mapping[k] = 4 * self.n + v

[ross-bott]

Hi there! I have marked this issue as stale because it has not had activity for 45 days. Consider the following options:

• If the issue refers to a large task, break it in smaller issues that can be solved in less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.

[raphaeltimbo]

We have a PR open to solve this (#407).

[rodrigomoliveira1]

PR #407 has been merged and closes this issue.