Stiffness matrix for shaft tapered element is not returning the expected results

[raphaeltimbo]

We are having problems with the stiffness matrix for shaft tapered elements. PR #233 implemented some tests we are skipping for now because of these issues.

@rodrigomoliveira1 , could you describe here the reference for the implementation and what tests did you carry so we can try discuss a solution for this problem?

[rodrigomoliveira1]

The matrices for shaft tapered elements are based in this article:

Genta, G., and Gugliotta, A. (1988). A conical element for finite element rotor dynamics, Journal of Sound and Vibration 120,175-182.

Following these matrices, I have implemented some tests (PR #217 and PR #233). Besides some basic tests to check the correct functioning of each method, I did a simple test that compares a cylindrical element built from class ShaftElement and a a cylindrical element built from ShaftTaperedElement. The expected results were numerically equal matricess. The mass (M), the gyroscopic (G) and the damping (C) matrices are ok, however it's not happening for the stiffness matrix (K).

@pytest.fixture
def tap2():
    #  Timoshenko element
    L = 0.4
    i_d_l = 0.
    i_d_r = 0.
    o_d_l = 0.25
    o_d_r = 0.25

    return ShaftTaperedElement(
        L, i_d_l, i_d_r, o_d_l, o_d_r, steel, shear_effects=True, rotary_inertia=True
    )

@pytest.fixture
def tim2():
    #  Timoshenko element
    le_ = 0.4
    i_d_ = 0
    o_d_ = 0.25
    return ShaftElement(
            le_, i_d_, o_d_, steel, rotary_inertia=True, shear_effects=True
    )

def test_match_mass_matrix(tap2, tim2):
    M_tap = tap2.M()
    M_tim = tim2.M()
    assert_almost_equal(M_tap, M_tim, decimal=5)

@pytest.mark.skip(reason="inconsistencies in stiffness matrix")
def test_match_stiffness_matrix(tap2, tim2):
    K_tap = tap2.K()
    K_tim = tim2.K()
    assert_almost_equal(K_tap, K_tim, decimal=5)

def test_match_gyroscopic_matrix(tap2, tim2):
    G_tap = tap2.G()
    G_tim = tim2.G()
    assert_almost_equal(G_tap, G_tim, decimal=5)

[raphaeltimbo]

I will update here what I am trying to test to solve this issue. I have tried different shear_method_calc ("cowper" or "hutchinson") and the comparison for K are still off, but it looks like "hutchinson" will give a smaller error.

For example, the value for K at the [0, 0] element is 2.40037e+09 for tap2.K() and 4.08082e+09 for tim2.K(). When we change to "hutchinson" we get 3.73722e+09 for tap2.K() and 3.66603e+09 for tim2.K().

[raphaeltimbo]

For reference, this is the stiffness matrix for the element obtained from XLTRC2.

Element properties (ross steel): density: 7810 kg/m3 E: 211e9 N/m2 G: 81.2e9 N/m2

L=0.4, od_left=0.25, id_left=0, od_right=0.25, id_right=0

4.08E+09	0	0	8.16E+08	-4.1E+09	0	0	8.16E+08
0	4.08E+09	-8.2E+08	0	0	-4.1E+09	-8.2E+08	0
0	-8.2E+08	2.64E+08	0	0	8.16E+08	62086082	0
8.16E+08	0	0	2.64E+08	-8.2E+08	0	0	62086082
-4.1E+09	0	0	-8.2E+08	4.08E+09	0	0	-8.2E+08
0	-4.1E+09	8.16E+08	0	0	4.08E+09	8.16E+08	0
0	-8.2E+08	62086082	0	0	8.16E+08	2.64E+08	0
8.16E+08	0	0	62086082	-8.2E+08	0	0	2.64E+08

[rodrigomoliveira1]

I have finally fixed the problems related to the stiffness matrix for shaft tapered elements.

The reference paper has some disagreement with what Friswell presents in his book and in his codes (which are based in the same paper we are using).

• First of all, some terms in the stiffness matrix have their signs missplaced.

• The element stiffness matrix is built from 2 separated matrices (K1 and K2). K1 is proportional to the Young's Modulos (E) and the Inertia (Ie), while K2 should be proportional to the Shear Modulus (G) and the Area (A).

  
      K1 = [
          [  k1,     0,       0,    L*k2,   -k1,    0,       0,    L*k3],
          [   0,    k1,   -L*k2,       0,     0,  -k1,   -L*k3,       0],
          [   0, -L*k2, L**2*k4,       0,     0, L*k2, L**2*k5,       0],
          [L*k2,     0,       0, L**2*k4, -L*k2,    0,       0, L**2*k5],
          [ -k1,     0,       0,   -L*k2,    k1,    0,       0,   -L*k3],
          [   0,   -k1,    L*k2,       0,     0,   k1,    L*k3,       0],
          [   0, -L*k3, L**2*k5,       0,     0, L*k3, L**2*k6,       0],
          [L*k3,     0,       0, L**2*k5, -L*k3,    0,       0, L**2*k6],
      ]
  
      K2 = [
          [  k7,     0,       0,    L*k8,   -k7,     0,       0,    L*k8],
          [   0,    k7,   -L*k8,       0,     0,   -k7,   -L*k8,       0],
          [   0, -L*k8, L**2*k9,       0,     0,  L*k8, L**2*k9,       0],
          [L*k8,     0,       0, L**2*k9, -L*k8,     0,       0, L**2*k9],
          [ -k7,     0,       0,   -L*k8,    k7,     0,       0,   -L*k8],
          [   0,   -k7,    L*k8,       0,     0,    k7,    L*k8,       0],
          [   0, -L*k8, L**2*k9,       0,     0,  L*k8, L**2*k9,       0],
          [L*k8,     0,       0, L**2*k9, -L*k8,     0,       0, L**2*k9],
      ]

K1 = K1 * (E * Ie * K / (105 * L ** 3 * (1 + phi) ** 2)) K2 = K2 * (G * A * phi ** 2 * K2 / (12 * kappa * L * (1 + phi) ** 2)) K = K1 + K2

However, checking again what Friswell has done in his code, the correct formula (and what is working fine now) is:

K = E * Ie / (105 * L ** 3 * (1 + phi) ** 2) * (K1 + 105 * phi * K2)

• The paper also says that the term K1[0][0] is k1 = 630 * a2 + 504 * b2 + 441 * gama + 396 * delta, but it actually lacks an independent value. The correct form is k1 = 1260 + 630 * a2 + 504 * b2 + 441 * gama + 396 * delta. This make sense, since if we set both ends with equal outer and inner radius, then a2, b2, gama and delta will be zero, and K1[0][0] would be zero too without this independent value (1260).

[rodrigomoliveira1]

PR #268 closes this issue