Repository: ross-rotordynamics/ross
Pull request #268 (2019-09-03 14:48:20)

Fix inconsistencies in the stiffness matrix
- Fix some typos in geometry coefficients;
- Fix stiffness matrix formula for shaft tapered elements;
- Fix docstrings for M(), K(), G() in ShaftTaperedElement;
- Add better docstrings to ShaftTaperedElement.

Fix tests to match the changes done in shaft_element.py:
- Fix stiffness matrix test;
- Fix element attribute tests ;
- Fix comparation tests.


Patch:
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 6d45226d..8cd6519c 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -695,6 +695,7 @@ class ShaftTaperedElement(Element):
     .. math:: [x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]^T
     Where :math:`\alpha_0` and :math:`\alpha_1` are the bending on the yz plane
     and :math:`\beta_0` and :math:`\beta_1` are the bending on the xz plane.
+
     Parameters
     ----------
     L : float
@@ -730,6 +731,44 @@ class ShaftTaperedElement(Element):
     shear_method_calc : string, optional
         Determines which shear calculation method the user will adopt
         Default is 'cowper'
+
+    Returns
+    -------
+    A shaft element
+
+    Attributes
+    ----------
+    Poisson : float
+        Poisson coefficient for the element.
+    A : float
+        Element mid-section area.
+    Ie : float
+        Ie is the second moment of area of the mid cross section about
+        the neutral plane Ie = pi*r**2/4
+    phi : float
+        Constant that is used according to [1] to consider rotary
+        inertia and shear effects. If these are not considered phi=0.
+    References
+    ----------
+    .. [1] 'Dynamics of Rotating Machinery' by MI Friswell, JET Penny, SD Garvey
+       & AW Lees, published by Cambridge University Press, 2010 pp. 158-166.
+    Examples
+    --------
+    >>> from ross.materials import steel
+    >>> Euler_Bernoulli_Element = ShaftTaperedElement(
+    ...                         L=0.5, i_d_l=0.05, i_d_r=0.05, o_d_l=0.1,
+    ...                         o_d_r=0.15, material=steel,
+    ...                         rotary_inertia=False,
+    ...                         shear_effects=False)
+    >>> Euler_Bernoulli_Element.phi
+    0
+    >>> Timoshenko_Element = ShaftTaperedElement(
+    ...                         L=0.5, i_d_l=0.05, i_d_r=0.05, o_d_l=0.1,
+    ...                         o_d_r=0.15, material=steel,
+    ...                         rotary_inertia=True,
+    ...                         shear_effects=True)
+    >>> Timoshenko_Element.phi
+    0.21026260191948362
     """
 
     def __init__(
@@ -770,7 +809,6 @@ class ShaftTaperedElement(Element):
 
         self.L = float(L)
 
-        # diameters
         self.o_d = (float(o_d_l) + float(o_d_r)) / 2
         self.i_d = (float(i_d_l) + float(i_d_r)) / 2
         self.i_d_l = float(i_d_l)
@@ -805,7 +843,7 @@ class ShaftTaperedElement(Element):
         b1 = np.pi * (delta_ro ** 2 - delta_ri ** 2) / A_l
         b2 = 3 * np.pi * (roj ** 2 * delta_ro ** 2 - rij ** 2 * delta_ri ** 2) / (2 * Ie_l)
         gama = np.pi * (roj * delta_ro ** 3 - rij * delta_ri ** 3) / Ie_l
-        delta = np.pi * (delta_ro ** 4 - delta_ri ** 4) / 4 * Ie_l
+        delta = np.pi * (delta_ro ** 4 - delta_ri ** 4) / (4 * Ie_l)
 
         self.a1 = a1
         self.a2 = a2
@@ -835,7 +873,7 @@ class ShaftTaperedElement(Element):
         # hutchinson - kappa as per Hutchinson (2001)
         # cowper - kappa as per Cowper (1996)
         if shear_effects:
-            r = i_d_l / o_d_l
+            r = ((i_d_l + i_d_r) / 2) / ((o_d_l + o_d_r) / 2)
             r2 = r * r
             r12 = (1 + r2) ** 2
             if shear_method_calc == "hutchinson":
@@ -941,9 +979,11 @@ class ShaftTaperedElement(Element):
 
     def M(self):
         r"""Mass matrix for an instance of a shaft tapered element.
+
         Returns
         -------
         Mass matrix for the shaft tapered element.
+
         Examples
         --------
         >>> Timoshenko_Element = ShaftTaperedElement(
@@ -952,10 +992,10 @@ class ShaftTaperedElement(Element):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.M()[:4, :4]
-        array([[11.35295267,  0.        ,  0.        ,  0.85960486],
-               [ 0.        , 11.35295267, -0.85960486,  0.        ],
-               [ 0.        , -0.85960486,  0.08653475,  0.        ],
-               [ 0.85960486,  0.        ,  0.        ,  0.08653475]])
+        array([[11.36986417,  0.        ,  0.        ,  0.86197637],
+               [ 0.        , 11.36986417, -0.86197637,  0.        ],
+               [ 0.        , -0.86197637,  0.08667495,  0.        ],
+               [ 0.86197637,  0.        ,  0.        ,  0.08667495]])
         """
         phi = self.phi
         L = self.L
@@ -1090,9 +1130,11 @@ class ShaftTaperedElement(Element):
 
     def K(self):
         r"""Stiffness matrix for an instance of a shaft tapered element.
+
         Returns
         -------
         Stiffness matrix for the shaft tapered element.
+
         Examples
         --------
         >>> from ross.materials import steel
@@ -1102,13 +1144,13 @@ class ShaftTaperedElement(Element):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.K()[:4, :4]/1e6
-        array([[211.30832227,   0.        ,   0.        ,  54.96237663],
-               [  0.        , 211.30832227, -54.96237663,   0.        ],
-               [  0.        , -54.96237663,  15.71572869,   0.        ],
-               [ 54.96237663,   0.        ,   0.        ,  15.71572869]])
+        array([[209.25641985,   0.        ,   0.        ,  38.62129051],
+               [  0.        , 209.25641985, -38.62129051,   0.        ],
+               [  0.        , -38.62129051,  11.56619973,   0.        ],
+               [ 38.62129051,   0.        ,   0.        ,  11.56619973]])
         """
-        phi = self.phi
         L = self.L
+        phi = self.phi
         a1 = self.a1
         a2 = self.a2
         b1 = self.b1
@@ -1116,10 +1158,9 @@ class ShaftTaperedElement(Element):
         delta = self.delta
         gama = self.gama
         Ie_l = self.Ie_l
-        A_l = self.A_l
 
         # fmt: off
-        k1 = 630 * a2 + 504 * b2 + 441 * gama + 396 * delta
+        k1 = 1260 + 630 * a2 + 504 * b2 + 441 * gama + 396 * delta
         k2 = (
             630
             + 210 * a2
@@ -1161,7 +1202,7 @@ class ShaftTaperedElement(Element):
         k8 = 6 + 3 * a1 + 2 * b1
         k9 = 3 + 1.5 * a1 + b1
 
-        K = np.array([
+        K1 = np.array([
             [  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],
@@ -1172,20 +1213,18 @@ class ShaftTaperedElement(Element):
             [L*k3,     0,       0, L**2*k5, -L*k3,    0,       0, L**2*k6],
         ])
 
-        K = self.material.E * Ie_l * K / (105 * L ** 3 * (1 + phi) ** 2)
-
         K2 = np.array([
             [  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],
+            [ -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],
         ])
-        K2 = self.material.G_s * A_l * phi ** 2 * K2 / (12 * self.kappa * L * (1 + phi) ** 2)
-        K = K + K2
+
+        K = self.material.E * Ie_l / (105 * L ** 3 * (1 + phi) ** 2) * (K1 + 105 * phi * K2)
         # fmt: on
 
         return K
@@ -1200,6 +1239,17 @@ class ShaftTaperedElement(Element):
 
         Examples
         --------
+        >>> from ross.materials import steel
+        >>> Timoshenko_Element = ShaftTaperedElement(
+        ...                         0.5, 0.05, 0.05, 0.1,
+        ...                         0.15, steel,
+        ...                         rotary_inertia=True,
+        ...                         shear_effects=True)
+        >>> Timoshenko_Element.C()[:4, :4]
+        array([[0., 0., 0., 0.],
+               [0., 0., 0., 0.],
+               [0., 0., 0., 0.],
+               [0., 0., 0., 0.]])
         """
         C = np.zeros((8, 8))
 
@@ -1207,9 +1257,11 @@ class ShaftTaperedElement(Element):
 
     def G(self):
         """Gyroscopic matrix for an instance of a shaft tapred element.
+
         Returns
         -------
         Gyroscopic matrix for the shaft tapered element.
+ 
         Examples
         --------
         >>> from ross.materials import steel
@@ -1220,10 +1272,10 @@ class ShaftTaperedElement(Element):
         ...                         rotary_inertia=True,
         ...                         shear_effects=True)
         >>> Timoshenko_Element.G()[:4, :4]
-        array([[ 0.        ,  0.2998902 , -0.00955058,  0.        ],
-               [-0.2998902 ,  0.        ,  0.        , -0.00955058],
-               [ 0.00955058,  0.        ,  0.        ,  0.00663842],
-               [ 0.        ,  0.00955058, -0.00663842,  0.        ]])
+        array([[ 0.        ,  0.30940809, -0.01085902,  0.        ],
+               [-0.30940809,  0.        ,  0.        , -0.01085902],
+               [ 0.01085902,  0.        ,  0.        ,  0.0067206 ],
+               [ 0.        ,  0.01085902, -0.0067206 ,  0.        ]])
         """
 
         G = np.zeros((8, 8))


Test patch:
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index ee1731d3..76183a51 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -153,6 +153,8 @@ def test_from_table():
         assert shaft[9].n == 8
         assert_allclose(shaft[8].material.E, 206842718795.05, rtol=3e-06)
         assert_allclose(shaft[9].material.E, 6894.75, atol=0.008)
+
+
 # Shaft Tapered Element tests
 @pytest.fixture
 def tap_tim():
@@ -192,49 +194,51 @@ def test_parameters_tap_tim(tap_tim):
     assert tap_tim.material.rho == 7810
     assert_almost_equal(tap_tim.material.Poisson, 0.29926108)
     assert_almost_equal(tap_tim.A, 0.024052818754)
-    assert_almost_equal(tap_tim.Ie * 1e5, 4.44854429)
-    assert_almost_equal(tap_tim.phi, 0.40668276)
+    assert_almost_equal(tap_tim.Ie * 1e5,  4.60385984)
+    assert_almost_equal(tap_tim.phi, 0.4208816002)
 
 
 def test_mass_matrix_tap_tim(tap_tim):
     # fmt: off
-    M0e_tim = np.array([[41.85725794,  0.        ,  0.        ,  1.8603991 , 10.27411395,  0.        ,  0.        , -1.14850813],
-                        [ 0.        , 41.85725794, -1.8603991 ,  0.        ,  0.        , 10.27411395,  1.14850813,  0.        ],
-                        [ 0.        , -1.8603991 ,  0.18670638,  0.        ,  0.        , -1.04155143, -0.10238917,  0.        ],
-                        [ 1.8603991 ,  0.        ,  0.        ,  0.18670638,  1.04155143,  0.        ,  0.        , -0.10238917],
-                        [10.27411395,  0.        ,  0.        ,  1.04155143, 17.33598969,  0.        ,  0.        , -1.14296052],
-                        [ 0.        , 10.27411395, -1.04155143,  0.        ,  0.        , 17.33598969,  1.14296052,  0.        ],
-                        [ 0.        ,  1.14850813, -0.10238917,  0.        ,  0.        ,  1.14296052,  0.09198448,  0.        ],
-                        [-1.14850813,  0.        ,  0.        , -0.10238917, -1.14296052,  0.        ,  0.        ,  0.09198448]])
+    M0e_tim = np.array([
+            [41.84146162,  0.        ,  0.        ,  1.85735573, 10.27465893,  0.        ,  0.        , -1.16147981],
+            [ 0.        , 41.84146162, -1.85735573,  0.        ,  0.        , 10.27465893,  1.16147981,  0.        ],
+            [ 0.        , -1.85735573,  0.18686024,  0.        ,  0.        , -1.04154454, -0.10311193,  0.        ],
+            [ 1.85735573,  0.        ,  0.        ,  0.18686024,  1.04154454,  0.        ,  0.        , -0.10311193],
+            [10.27465893,  0.        ,  0.        ,  1.04154454, 17.35069604,  0.        ,  0.        , -1.1330391 ],
+            [ 0.        , 10.27465893, -1.04154454,  0.        ,  0.        , 17.35069604,  1.1330391 ,  0.        ],
+            [ 0.        ,  1.16147981, -0.10311193,  0.        ,  0.        ,  1.1330391 ,  0.09851808,  0.        ],
+            [-1.16147981,  0.        ,  0.        , -0.10311193, -1.1330391 ,  0.        ,  0.        ,  0.09851808]])
     # fmt: on
     assert_almost_equal(tap_tim.M(), M0e_tim, decimal=5)
 
 
-@pytest.mark.skip(reason="inconsistencies in stiffness matrix")
 def test_stiffness_matrix_tap_tim(tap_tim):
     # fmt: off
-    K0e_tim = np.array([[-19.5879458 ,   0.        ,   0.        ,   5.57911856,  19.5879458 ,   0.        ,   0.        ,   1.92061626],
-                        [  0.        , -19.5879458 ,  -5.57911856,   0.        ,   0.        ,  19.5879458 ,  -1.92061626,   0.        ],
-                        [  0.        ,  -5.57911856,   1.78916039,   0.        ,   0.        ,   5.57911856,   0.44248703,   0.        ],
-                        [  5.57911856,   0.        ,   0.        ,   1.78916039,  -5.57911856,   0.        ,   0.        ,   0.44248703],
-                        [ 19.5879458 ,   0.        ,   0.        ,  -5.57911856, -19.5879458 ,   0.        ,   0.        ,   0.03401395],
-                        [  0.        ,  19.5879458 ,   5.57911856,   0.        ,   0.        , -19.5879458 ,  -0.03401395,   0.        ],
-                        [  0.        ,  -1.92061626,   0.44248703,   0.        ,   0.        ,  -0.03401395,   0.32575947,   0.        ],
-                        [  1.92061626,   0.        ,   0.        ,   0.44248703,   0.03401395,   0.        ,   0.        ,   0.32575947]])
+    K0e_tim = np.array([
+            [ 23.3411494 ,   0.        ,   0.        ,   6.2946873 , -23.3411494 ,   0.        ,   0.        ,   3.04177246],
+            [  0.        ,  23.3411494 ,  -6.2946873 ,   0.        ,   0.        , -23.3411494 ,  -3.04177246,   0.        ],
+            [  0.        ,  -6.2946873 ,   1.91793266,   0.        ,   0.        ,   6.2946873 ,   0.59994226,   0.        ],
+            [  6.2946873 ,   0.        ,   0.        ,   1.91793266,  -6.2946873 ,   0.        ,   0.        ,   0.59994226],
+            [-23.3411494 ,   0.        ,   0.        ,  -6.2946873 ,  23.3411494 ,   0.        ,   0.        ,  -3.04177246],
+            [  0.        , -23.3411494 ,   6.2946873 ,   0.        ,   0.        ,  23.3411494 ,   3.04177246,   0.        ],
+            [  0.        ,  -3.04177246,   0.59994226,   0.        ,   0.        ,   3.04177246,   0.61676672,   0.        ],
+            [  3.04177246,   0.        ,   0.        ,   0.59994226,  -3.04177246,   0.        ,   0.        ,   0.61676672]])
     # fmt: on
     assert_almost_equal(tap_tim.K() / 1e8, K0e_tim, decimal=5)
 
 
 def test_gyroscopic_matrix_tap_tim(tap_tim):
     # fmt: off
-    G0e_tim = np.array([[ 0.        ,  1.19360236,  0.15864756,  0.        ,  0.        , -1.19360236, -0.05013491,  0.        ],
-                        [-1.19360236,  0.        ,  0.        ,  0.15864756,  1.19360236,  0.        ,  0.        , -0.05013491],
-                        [-0.15864756,  0.        ,  0.        ,  0.11754246,  0.15864756,  0.        ,  0.        , -0.01449657],
-                        [ 0.        , -0.15864756, -0.11754246,  0.        ,  0.        ,  0.15864756,  0.01449657,  0.        ],
-                        [ 0.        , -1.19360236, -0.15864756,  0.        ,  0.        ,  1.19360236,  0.05013491,  0.        ],
-                        [ 1.19360236,  0.        ,  0.        , -0.15864756, -1.19360236,  0.        ,  0.        ,  0.05013491],
-                        [ 0.05013491,  0.        ,  0.        , -0.01449657, -0.05013491,  0.        ,  0.        , -0.00213197],
-                        [ 0.        ,  0.05013491,  0.01449657,  0.        ,  0.        , -0.05013491,  0.00213197,  0.        ]])
+    G0e_tim = np.array([
+            [ 0.        ,  1.23853234,  0.15544731,  0.        ,  0.        , -1.23853234, -0.03173553,  0.        ],
+            [-1.23853234,  0.        ,  0.        ,  0.15544731,  1.23853234,  0.        ,  0.        , -0.03173553],
+            [-0.15544731,  0.        ,  0.        ,  0.11850408,  0.15544731,  0.        ,  0.        , -0.01563679],
+            [ 0.        , -0.15544731, -0.11850408,  0.        ,  0.        ,  0.15544731,  0.01563679,  0.        ],
+            [ 0.        , -1.23853234, -0.15544731,  0.        ,  0.        ,  1.23853234,  0.03173553,  0.        ],
+            [ 1.23853234,  0.        ,  0.        , -0.15544731, -1.23853234,  0.        ,  0.        ,  0.03173553],
+            [ 0.03173553,  0.        ,  0.        , -0.01563679, -0.03173553,  0.        ,  0.        ,  0.01089192],
+            [ 0.        ,  0.03173553,  0.01563679,  0.        ,  0.        , -0.03173553, -0.01089192,  0.        ]])
     # fmt: on
     assert_almost_equal(tap_tim.G(), G0e_tim, decimal=5)
 
@@ -270,7 +274,6 @@ def test_match_mass_matrix(tap2, tim2):
     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()