Error in tapered element stiffness formula

This error has been raised by HaukeWittich in https://github.com/ross-rotordynamics/ross/discussions/878

^{Originally posted by **HaukeWittich** March 17, 2022} Hello, I have the impression that there is a small bug in the shaft element part, more specifically in the second part of the stiffness matrix Kh2, which reflects the shear deformation. This part is depending on the surface A. For a tapered element the surface is variable described via a quadratic polynomial function with a constant A_l, which is the surface of the left hand side of the element. Further there is a transformation with phi, which is defined with the mean surface A in the middle of the element. The quotient between this two surfaces is not respected in my eyes. Thank you.

To give more details on this error, our reference to implement the tapered element was Genta, G., and A. Gugliotta. "A conical element for finite element rotor dynamics." Journal of Sound and Vibration 120.1 (1988): 175-182, which we can read it here

In the appendix A of this reference, we can find the expression for the stiffness matrix:

We have implemented this equation in the code with some transformation, like this:

K = self.material.E * Ie_l / (105 * L ** 3 * (1 + phi) ** 2) * (K1 + 105 * phi * K2)

Where K1 is the first matrix (left) of the equation and K2 is the last matrix.

So, to confirm whether there was an error in the value multiplying K2 I have performed the following calculations:

import ross as rs

steel = rs.materials.steel

sh = rs.ShaftElement(
    L=0.4,
    idl=0.0,
    idr=0.0,
    odl=0.25,
    odr=0.10,
    n=0,
    material=steel
)
print(
    "Transformed equation: ",
    ((sh.material.E * sh.Ie_l) / (105 * sh.L**3 * (1 + sh.phi) ** 2)) * sh.phi * 105,
)

chi = 1 / sh.kappa
print(
    "Paper equation: ",
    sh.material.G_s * sh.A_l * sh.phi**2 / (12 * chi * sh.L * (1 + sh.phi) ** 2),
)
-----------------------
Transformed equation:  131788214.46286184
Paper equation:  64576225.08680231

The ratio between the 'transformed' value and 'paper' value is exactly equal to the ratio between the area at the left of the element and at the center of the element.

To double check, I performed an analysis using the cylindrical code that we had validated on ROSS v0.2.0 with 1000 elements for a tapered shaft and compared with one element of the code with "Transformed equation" and "Paper equation":

Ccode using cylindrical elements with ROSS v0.2.0 and 1000 elements.

import ross as rs

steel = rs.materials.steel
L_total = 0.5
odl = 0.1
odr = 0.01
shaft_elements = []

N = 4
delta_d = (odl - odr) / N

odl_el = odl
odr_el = odl - delta_d

for i in range(N):
    L = L_total / N

    sh_el = rs.ShaftElement(
        n=i,
        idl=0.0,
        odl=odl_el,
        idr=0.0,
        odr=odr_el,
        material=steel,
        L=L,
    )

    if i == 0:
        sh1 = sh_el
    shaft_elements.append(sh_el)

    odl_el -= delta_d
    odr_el -= delta_d

bearing0 = rs.BearingElement(n=0, kxx=1e20, kyy=1e20, cxx=0)
bearing1 = rs.BearingElement(n=N, kxx=1e20, kyy=1e20, cxx=0)

rotor = rs.Rotor(shaft_elements=shaft_elements, bearing_elements=[bearing0, bearing1])

modal = rotor.run_modal(speed=0)
modal.wn

---> array([ 1554.93375489,  1554.93375587,  9499.70295089,  9499.70296119,
       19790.70688825, 19790.70771746])

Code with 4 tapered elements:

import ross as rs

steel = rs.materials.steel
L_total = 0.5
odl = 0.1
odr = 0.01
shaft_elements = []

N = 4
delta_d = (odl - odr) / N

odl_el = odl
odr_el = odl - delta_d

for i in range(N):
    L = L_total / N

    sh_el = rs.ShaftElement(
        n=i,
        idl=0.0,
        odl=odl_el,
        idr=0.0,
        odr=odr_el,
        material=steel,
        L=L,
    )
    shaft_elements.append(sh_el)

    odl_el -= delta_d
    odr_el -= delta_d

bearing0 = rs.BearingElement(n=0, kxx=1e20, kyy=1e20, cxx=0)
bearing1 = rs.BearingElement(n=N, kxx=1e20, kyy=1e20, cxx=0)

rotor = rs.Rotor(shaft_elements=shaft_elements, bearing_elements=[bearing0, bearing1])

modal = rotor.run_modal(speed=0)
modal.wn

---> results with 'Transformed equation': array([ 1643.7248728 ,  1643.72487315, 10014.71779303, 10014.71812789,
       21599.29763002, 21599.29868402])

---> results with 'Paper equation': array([ 1630.75091817,  1630.75091846,  9899.54380766,  9899.54396821,
       21074.45305736, 21074.45655191])

As we can see, we are closer with the 'Paper equation'. I am going to make a PR to fix this issue.

petrobras / ross

Error in tapered element stiffness formula #880