Repository: ross-rotordynamics/ross
Pull request #217 (2019-07-25 18:15:42)

 Add tapered shaft element 
- Add new element:  tapered shaft element
- Add class ShaftTaperedElement
- This class have, basically, the same inputs than "ShaftElement". The difference relies on entering the diameters of each side of the element (left and right).
- Patches methods were adapted to draw the conical shape of the element.

The matrices generated for this element follow the reference below.
> Genta, G., and Gugliotta, A. (1988). A conical element for finite element rotor dynamics, Journal of Sound and Vibration 120,175-182.


Patch:
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 0517fe97..15a497b3 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -245,7 +245,7 @@ class Rotor(object):
         self.df_shaft = df_shaft
 
         # check consistence for disks and bearings location
-        if df.n_l.max() > df[df.type == "ShaftElement"].n_r.max():
+        if df.n_l.max() > df_shaft.n_r.max():
             raise ValueError("Trying to set disk or bearing outside shaft")
 
         self.df = df
@@ -1168,7 +1168,7 @@ class Rotor(object):
         try:
             max_diameter = max([disk.o_d for disk in self.disk_elements])
         except (ValueError, AttributeError):
-            max_diameter = max([shaft.o_d for shaft in self.shaft_elements])
+            max_diameter = max([shaft.o_d_l for shaft in self.shaft_elements])
 
         ax.set_ylim(-1.2 * max_diameter, 1.2 * max_diameter)
         ax.axis("equal")
@@ -1176,8 +1176,6 @@ class Rotor(object):
         ax.set_ylabel("Shaft radius (m)")
 
         # plot nodes
-        text = []
-        x_pos = []
         for node, position in enumerate(self.nodes_pos[::nodes]):
             ax.plot(
                 position,
@@ -1261,11 +1259,6 @@ class Rotor(object):
         #  plot shaft centerline
         shaft_end = self.nodes_pos[-1]
 
-        try:
-            max_diameter = max([disk.o_d for disk in self.disk_elements])
-        except (ValueError, AttributeError):
-            max_diameter = max([shaft.o_d for shaft in self.shaft_elements])
-
         # bokeh plot - create a new plot
         bk_ax = figure(
             tools="pan, wheel_zoom, reset, save",
@@ -1295,7 +1288,7 @@ class Rotor(object):
         x_pos = []
         for node, position in enumerate(self.nodes_pos[::nodes]):
             # bokeh plot
-            text.append(str(node))
+            text.append(str(node*nodes))
             x_pos.append(position)
 
         # bokeh plot - plot nodes
@@ -1322,13 +1315,17 @@ class Rotor(object):
         # plot shaft elements
         for sh_elm in self.shaft_elements:
             position = self.nodes_pos[sh_elm.n]
-            sh_elm.bokeh_patch(position, SR, bk_ax)
+            hover = sh_elm.bokeh_patch(position, SR, bk_ax)
+
+        bk_ax.add_tools(hover)
 
         # plot disk elements
         for disk in self.disk_elements:
             position = (self.nodes_pos[disk.n], self.nodes_o_d[disk.n] / 2)
             length = min(self.nodes_le)
-            disk.bokeh_patch(position, length, bk_ax)
+            hover = disk.bokeh_patch(position, length, bk_ax)
+
+        bk_ax.add_tools(hover)
 
         # plot bearings
         for bearing in self.bearing_seal_elements:
@@ -1340,7 +1337,7 @@ class Rotor(object):
 
         return bk_ax
 
-    def plot_rotor(self, *args, plot_type="matplotlib", **kwargs):
+    def plot_rotor(self, nodes=1, *args, plot_type="matplotlib", **kwargs):
         """Plots a rotor object.
 
         This function will take a rotor object and plot its shaft,
@@ -1372,9 +1369,9 @@ class Rotor(object):
         <matplotlib.axes...
         """
         if plot_type == "matplotlib":
-            return self._plot_rotor_matplotlib(*args, **kwargs)
+            return self._plot_rotor_matplotlib(nodes=nodes, *args, **kwargs)
         elif plot_type == "bokeh":
-            return self._plot_rotor_bokeh(*args, **kwargs)
+            return self._plot_rotor_bokeh(nodes=nodes, *args, **kwargs)
         else:
             raise ValueError(f"{plot_type} is not a valid plot type.")
 
diff --git a/ross/shaft_element.py b/ross/shaft_element.py
index 558b9e54..97e2c141 100644
--- a/ross/shaft_element.py
+++ b/ross/shaft_element.py
@@ -583,11 +583,9 @@ class ShaftElement(Element):
         height = self.o_d / 2 - self.i_d / 2
         if self.n in SR:
             mpl_color = "yellow"
-            bk_color = "yellow"
             legend = "Shaft - Slenderness Ratio < 1.6"
         else:
             mpl_color = self.color
-            bk_color = bokeh_colors[2]
             legend = "Shaft"
 
         #  matplotlib - plot the upper half of the shaft
@@ -645,8 +643,10 @@ class ShaftElement(Element):
                 left=[position],
                 right=[position + self.L],
                 elnum=[self.n],
-                out_d=[self.o_d],
-                in_d=[self.i_d],
+                out_d_l=[self.o_d_l],
+                out_d_r=[self.o_d_r],
+                in_d_l=[self.i_d_l],
+                in_d_r=[self.i_d_r],
                 length=[self.L],
                 mat=[self.material.name],
             )
@@ -659,8 +659,10 @@ class ShaftElement(Element):
                 left=[position],
                 right=[position + self.L],
                 elnum=[self.n],
-                out_d=[self.o_d],
-                in_d=[self.i_d],
+                out_d_l=[self.o_d_l],
+                out_d_r=[self.o_d_r],
+                in_d_l=[self.i_d_l],
+                in_d_r=[self.i_d_r],
                 length=[self.L],
                 mat=[self.material.name],
             )
@@ -695,15 +697,16 @@ class ShaftElement(Element):
         hover = HoverTool(names=["u_shaft", "l_shaft"])
         hover.tooltips = [
             ("Element Number :", "@elnum"),
-            ("Outer Diameter :", "@out_d"),
-            ("Internal Diameter :", "@in_d"),
+            ("Left Outer Diameter :", "@out_d_l"),
+            ("Left Inner Diameter :", "@in_d_l"),
+            ("Right Outer Diameter :", "@out_d_r"),
+            ("Right Inner Diameter :", "@in_d_r"),
             ("Element Length :", "@length"),
             ("Material :", "@mat"),
         ]
         hover.mode = "mouse"
 
-        if len(bk_ax.hover) == 0:
-            bk_ax.add_tools(hover)
+        return hover
 
     @classmethod
     def section(
@@ -777,3 +780,750 @@ class ShaftElement(Element):
         ]
 
         return elements
+
+
+class ShaftTaperedElement(Element):
+    r"""A shaft tapered element.
+    This class will create a shaft element that may take into
+    account shear, rotary inertia and gyroscopic effects.
+    The matrices will be defined considering the following local
+    coordinate vector:
+    .. 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
+        Element length.
+    i_d_l : float
+        Inner diameter of the element at the left position.
+    i_d_r : float
+        Inner diameter of the element at the right position.
+    o_d_l : float
+        Outer diameter of the element at the left position.
+    o_d_r : float
+        Outer diameter of the element at the right position.
+    material : ross.material
+        Shaft material.
+    n : int, optional
+        Element number (coincident with it's first node).
+        If not given, it will be set when the rotor is assembled
+        according to the element's position in the list supplied to
+        the rotor constructor.
+    axial_force : float, optional
+        Axial force.
+    torque : float, optional
+        Torque.
+    shear_effects : bool, optional
+        Determine if shear effects are taken into account.
+        Default is True.
+    rotary_inertia : bool, optional
+        Determine if rotary_inertia effects are taken into account.
+        Default is True.
+    gyroscopic : bool, optional
+        Determine if gyroscopic effects are taken into account.
+        Default is True.
+    shear_method_calc : string, optional
+        Determines which shear calculation method the user will adopt
+        Default is 'cowper'
+    """
+
+    def __init__(
+        self,
+        L,
+        i_d_l,
+        i_d_r,
+        o_d_l,
+        o_d_r,
+        material,
+        n=None,
+        axial_force=0,
+        torque=0,
+        shear_effects=True,
+        rotary_inertia=True,
+        gyroscopic=True,
+        shear_method_calc="cowper",
+    ):
+
+        if type(material) is str:
+            os.chdir(Path(os.path.dirname(ross.__file__)))
+            self.material = Material.use_material(material)
+        else:
+            self.material = material
+
+        self.shear_effects = shear_effects
+        self.rotary_inertia = rotary_inertia
+        self.gyroscopic = gyroscopic
+        self.axial_force = axial_force
+        self.torque = torque
+        self._n = n
+        self.n_l = n
+        self.n_r = None
+        if n is not None:
+            self.n_r = n + 1
+
+        self.shear_method_calc = shear_method_calc
+
+        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)
+        self.o_d_l = float(o_d_l)
+        self.i_d_r = float(i_d_r)
+        self.o_d_r = float(o_d_r)
+        self.color = self.material.color
+
+        # A_l = cross section area from the left side of the element
+        A_l = np.pi * (o_d_l ** 2 - i_d_l ** 2) / 4
+        self.A_l = A_l
+
+        # Second moment of area of the cross section from the left side
+        # of the element
+        Ie_l = np.pi * (o_d_l ** 4 - i_d_l ** 4) / 64
+
+        outer = self.o_d_l ** 2 + self.o_d_l * self.o_d_r + self.o_d_r ** 2
+        inner = self.i_d_l ** 2 + self.i_d_l * self.i_d_r + self.i_d_r ** 2
+        self.volume = np.pi * (self.L / 12) * (outer - inner)
+        self.m = self.material.rho * self.volume
+
+        roj = o_d_l / 2
+        rij = i_d_l / 2
+        rok = o_d_r / 2
+        rik = i_d_r / 2
+
+        # geometrical coefficients
+        delta_ro = rok - roj
+        delta_ri = rik - rij
+        a1 = 2 * np.pi * (roj * delta_ro - rij * delta_ri) / A_l
+        a2 = np.pi * (roj ** 3 * delta_ro - rij ** 3 * delta_ri) / Ie_l
+        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
+
+        self.a1 = a1
+        self.a2 = a2
+        self.b1 = b1
+        self.b2 = b2
+        self.gama = gama
+        self.delta = delta
+
+        # the area is calculated from the cross section located in the middle
+        # of the element
+        self.A = A_l * (1 + a1 * 0.5 + b1 * 0.5 ** 2)
+
+        # Ie is the second moment of area of the cross section - located in
+        # the middle of the element - about the neutral plane
+        Ie = Ie_l * (1 + a2 * 0.5 + b2 * 0.5 ** 2 + gama * 0.5 ** 3 + delta * 0.5 ** 4)
+        self.Ie = Ie
+        self.Ie_l = Ie_l
+
+        phi = 0
+
+        # Slenderness ratio of beam elements (G*A*L^2) / (E*I)
+        sld = (self.material.G_s * self.A * self.L ** 2) / (self.material.E * Ie)
+        self.slenderness_ratio = sld
+
+        # picking a method to calculate the shear coefficient
+        # List of avaible methods:
+        # hutchinson - kappa as per Hutchinson (2001)
+        # cowper - kappa as per Cowper (1996)
+        if shear_effects:
+            r = i_d_l / o_d_l
+            r2 = r * r
+            r12 = (1 + r2) ** 2
+            if shear_method_calc == "hutchinson":
+                # Shear coefficient (phi)
+                # kappa as per Hutchinson (2001)
+                # fmt: off
+                kappa = 6 * r12 * ((1 + self.material.Poisson) /
+                        ((r12 * (7 + 12 * self.material.Poisson + 4 * self.material.Poisson ** 2) +
+                        4 * r2 * (5 + 6 * self.material.Poisson + 2 * self.material.Poisson ** 2))))
+                # fmt: on
+            elif shear_method_calc == "cowper":
+                # kappa as per Cowper (1996)
+                # fmt: off
+                kappa = 6 * r12 * (
+                    (1 + self.material.Poisson)
+                    / (r12 * (7 + 6 * self.material.Poisson) + r2 * (20 + 12 * self.material.Poisson))
+                )
+                # fmt: on
+            else:
+                raise Warning(
+                    "This method of calculating shear coefficients is not implemented. See guide for futher informations."
+                )
+
+            # fmt: off
+            phi = 12 * self.material.E * self.Ie / (self.material.G_s * kappa * self.A * L ** 2)
+            # fmt: on
+            self.kappa = kappa
+        self.phi = phi
+
+    def __eq__(self, other):
+        if self.__dict__ == other.__dict__:
+            return True
+        else:
+            return False
+
+    def save(self, file_name):
+        data = self.load_data(file_name)
+        data["ShaftTaperedElement"][str(self.n)] = {
+            "L": self.L,
+            "i_d_l": self.i_d_l,
+            "i_d_r": self.i_d_r,
+            "o_d_l": self.o_d_l,
+            "o_d_r": self.o_d_r,
+            "material": self.material.name,
+            "n": self.n,
+            "axial_force": self.axial_force,
+            "torque": self.torque,
+            "shear_effects": self.shear_effects,
+            "rotary_inertia": self.rotary_inertia,
+            "gyroscopic": self.gyroscopic,
+            "shear_method_calc": self.shear_method_calc,
+        }
+        self.dump_data(data, file_name)
+
+    @staticmethod
+    def load(file_name="ShaftTaperedElement"):
+        shaft_elements = []
+        with open("ShaftTaperedElement.toml", "r") as f:
+            shaft_elements_dict = toml.load(f)
+            for element in shaft_elements_dict["ShaftTaperedElement"]:
+                shaft_elements.append(
+                    ShaftTaperedElement(
+                        **shaft_elements_dict["ShaftTaperedElement"][element]
+                    )
+                )
+        return shaft_elements
+
+    @property
+    def n(self):
+        return self._n
+
+    @n.setter
+    def n(self, value):
+        self._n = value
+        self.n_l = value
+        if value is not None:
+            self.n_r = value + 1
+
+    def dof_mapping(self):
+        return dict(x0=0, y0=1, alpha0=2, beta0=3, x1=4, y1=5, alpha1=6, beta1=7)
+
+    def __repr__(self):
+        return (
+            f"{self.__class__.__name__}"
+            f"(L={self.L:{0}.{5}}, i_d_l={self.i_d_l:{0}.{5}}, "
+            f"i_d_r={self.i_d_r:{0}.{5}}, o_d_l={self.o_d_l:{0}.{5}},  "
+            f"o_d_r={self.o_d_r:{0}.{5}}, material={self.material.name!r}, "
+            f"n={self.n})"
+        )
+
+    def __str__(self):
+        return (
+            f"\nElem. N:    {self.n}"
+            f"\nLenght:     {self.L:{10}.{5}}"
+            f"\nLeft Int. Diam.: {self.i_d_l:{10}.{5}}"
+            f"\nLeft Out. Diam.: {self.o_d_l:{10}.{5}}"
+            f"\nRight Int. Diam.: {self.i_d_r:{10}.{5}}"
+            f"\nRight Out. Diam.: {self.o_d_r:{10}.{5}}"
+            f'\n{35*"-"}'
+            f"\n{self.material}"
+            f"\n"
+        )
+
+    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(
+        ...                         0.5, 0.05, 0.05, 0.1,
+        ...                         0.15, steel,
+        ...                         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]])
+        """
+        phi = self.phi
+        L = self.L
+        a1 = self.a1
+        a2 = self.a2
+        b1 = self.b1
+        b2 = self.b2
+        delta = self.delta
+        gama = self.gama
+        Ie_l = self.Ie_l
+        A_l = self.A_l
+
+        m1 = (
+            (468 + 882 * phi + 420 * phi ** 2)
+            + a1 * (108 + 210 * phi + 105 * phi ** 2)
+            + b1 * (38 + 78 * phi + 42 * phi ** 2)
+        )
+        m2 = (
+            (66 + 115.5 * phi + 52.5 * phi ** 2)
+            + a1 * (21 + 40.5 * phi + 21 * phi ** 2)
+            + b1 * (8.5 + 18 * phi + 10.5 * phi ** 2)
+        )
+        m3 = (
+            (162 + 378 * phi + 210 * phi ** 2)
+            + a1 * (81 + 189 * phi + 105 * phi ** 2)
+            + b1 * (46 + 111 * phi + 63 * phi ** 2)
+        )
+        m4 = (
+            (39 + 94.5 * phi + 52.5 * phi ** 2)
+            + a1 * (18 + 40.5 * phi + 21 * phi ** 2)
+            + b1 * (9.5 + 21 * phi + 10.5 * phi ** 2)
+        )
+        m5 = (
+            (12 + 21 * phi + 10.5 * phi ** 2)
+            + a1 * (4.5 + 9 * phi + 5.25 * phi ** 2)
+            + b1 * (2 + 4.5 * phi + 3 * phi ** 2)
+        )
+        m6 = (
+            (39 + 94.5 * phi + 52.5 * phi ** 2)
+            + a1 * (21 + 54 * phi + 31.5 * phi ** 2)
+            + b1 * (12.5 + 34.5 * phi + 21 * phi ** 2)
+        )
+        m7 = (
+            (9 + 21 * phi + 10.5 * phi ** 2)
+            + a1 * (4.5 + 10.5 * phi + 5.25 * phi ** 2)
+            + b1 * (2.5 + 6 * phi + 3 * phi ** 2)
+        )
+        m8 = (
+            (468 + 882 * phi + 420 * phi ** 2)
+            + a1 * (360 + 672 * phi + 315 * phi ** 2)
+            + b1 * (290 + 540 * phi + 252 * phi ** 2)
+        )
+        m9 = (
+            (66 + 115.5 * phi + 52.5 * phi ** 2)
+            + a1 * (45 + 75 * phi + 31.5 * phi ** 2)
+            + b1 * (32.5 + 52.5 * phi + 21 * phi ** 2)
+        )
+        m10 = (
+            (12 + 21 * phi + 10.5 * phi ** 2)
+            + a1 * (7.5 + 12 * phi + 5.25 * phi ** 2)
+            + b1 * (5 + 7.5 * phi + 3 * phi ** 2)
+        )
+
+        # fmt: off
+        Mt = np.array([
+                [   m1,     0,        0,    L*m2,     m3,     0,        0,    -L*m4],
+                [    0,    m1,    -L*m2,       0,      0,    m3,     L*m4,        0],
+                [    0, -L*m2,  L**2*m5,       0,      0, -L*m6, -L**2*m7,        0],
+                [ L*m2,     0,        0,  L**2*m5,  L*m6,     0,        0, -L**2*m7],
+                [   m3,     0,        0,     L*m6,    m8,     0,        0,    -L*m9],
+                [    0,    m3,    -L*m6,        0,     0,    m8,     L*m9,        0],
+                [    0,  L*m4, -L**2*m7,        0,     0,  L*m9, L**2*m10,        0],
+                [-L*m4,     0,        0, -L**2*m7, -L*m9,     0,        0, L**2*m10],
+        ])
+        # fmt: on
+        M = self.material.rho * A_l * L * Mt / (1260 * (1 + phi) ** 2)
+
+        if self.rotary_inertia:
+            # fmt: off
+            m11 = 252 + 126 * a2 + 72 * b2 + 45 * gama + 30 * delta
+            m12 = (
+                21 - 105 * phi
+                + a2 * (21 - 42 * phi)
+                + b2 * (15 - 21 * phi)
+                + gama * (10.5 - 12 * phi)
+                + delta * (7.5 - 7.5 * phi)
+            )
+            m13 = (
+                21 - 105 * phi
+                - 63 * a2 * phi
+                - b2 * (6 + 42 * phi)
+                - gama * (7.5 + 30 * phi)
+                - delta * (7.5 + 22.5 * phi)
+            )
+            m14 = (
+                28 + 35 * phi + 70 * phi ** 2
+                + a2 * (7 - 7 * phi + 17.5 * phi ** 2)
+                + b2 * (4 - 7 * phi + 7 * phi ** 2)
+                + gama * (2.75 - 5 * phi + 3.5 * phi ** 2)
+                + delta * (2 - 3.5 * phi + 2 * phi ** 2)
+            )
+            m15 = (
+                7 + 35 * phi - 35 * phi ** 2
+                + a2 * (3.5 + 17.5 * phi - 17.5 * phi ** 2)
+                + b2 * (3 + 10.5 * phi - 10.5 * phi ** 2)
+                + gama * (2.75 + 7 * phi - 7 * phi ** 2)
+                + delta * (2.5 + 5 * phi - 5 * phi ** 2)
+            )
+            m16 = (
+                28 + 35 * phi + 70 * phi ** 2
+                + a2 * (21 + 42 * phi + 52.5 * phi ** 2)
+                + b2 * (18 + 42 * phi + 42 * phi ** 2)
+                + gama * (16.25 + 40 * phi + 35 * phi ** 2)
+                + delta * (15 + 37.5 * phi + 30 * phi ** 2)
+            )
+
+            Mr = np.array([
+                    [  m11,      0,         0,     L*m12,   -m11,     0,         0,     L*m13],
+                    [    0,    m11,    -L*m12,         0,      0,  -m11,    -L*m13,         0],
+                    [    0, -L*m12,  L**2*m14,         0,      0, L*m12, -L**2*m15,         0],
+                    [L*m12,      0,         0,  L**2*m14, -L*m12,     0,         0, -L**2*m15],
+                    [ -m11,      0,         0,    -L*m12,    m12,     0,         0,    -L*m13],
+                    [    0,   -m11,     L*m12,         0,      0,   m12,     L*m13,         0],
+                    [    0, -L*m13, -L**2*m15,         0,      0, L*m13,  L**2*m16,         0],
+                    [L*m13,      0,         0, -L**2*m15, -L*m13,     0,         0,  L**2*m16],
+            ])
+            # fmt: on
+            Mr = self.material.rho * Ie_l * Mr / (210 * L * (1 + phi) ** 2)
+            M = M + Mr
+
+        return M
+
+    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
+        >>> Timoshenko_Element = ShaftTaperedElement(
+        ...                         0.5, 0.05, 0.05, 0.1,
+        ...                         0.15, steel,
+        ...                         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]])
+        """
+        phi = self.phi
+        L = self.L
+        a1 = self.a1
+        a2 = self.a2
+        b1 = self.b1
+        b2 = self.b2
+        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
+        k2 = (
+            630
+            + 210 * a2
+            + 147 * b2
+            + 126 * gama
+            + 114 * delta
+            - phi * (105 * a2 + 105 * b2 + 94.5 * gama + 84 * delta)
+        )
+        k3 = (
+            630
+            + 420 * a2
+            + 357 * b2
+            + 315 * gama
+            + 282 * delta
+            + phi * (105 * a2 + 105 * b2 + 94.5 * gama + 84 * delta)
+        )
+        k4 = (
+            420 + 210 * phi + 105 * phi ** 2
+            + a2 * (105 + 52.5 * phi ** 2)
+            + b2 * (56 - 35 * phi + 35 * phi ** 2)
+            + gama * (42 - 42 * phi + 26.25 * phi ** 2)
+            + delta * (36 - 42 * phi + 21 * phi ** 2)
+        )
+        k5 = (
+            210 - 210 * phi - 105 * phi ** 2
+            + a2 * (105 - 105 * phi - 52.5 * phi ** 2)
+            + b2 * (91 - 70 * phi - 35 * phi ** 2)
+            + gama * (84 - 52.5 * phi - 26.25 * phi ** 2)
+            + delta * (78 - 42 * phi - 21 * phi ** 2)
+        )
+        k6 = (
+            420 + 210 * phi + 105 * phi ** 2
+            + a2 * (315 + 210 * phi + 52.5 * phi ** 2)
+            + b2 * (266 + 175 * phi + 35 * phi ** 2)
+            + gama * (231 + 147 * phi + 26.25 * phi ** 2)
+            + delta * (204 + 126 * phi + 21 * phi ** 2)
+        )
+        k7 = 12 + 6 * a1 + 4 * b1
+        k8 = 6 + 3 * a1 + 2 * b1
+        k9 = 3 + 1.5 * a1 + b1
+
+        K = 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],
+            [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],
+        ])
+
+        K = self.material.E * Ie_l * K / (105 * L ** 3 * (1 + phi) ** 2)
+
+        if self.shear_effects:
+            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],
+            ])
+            K2 = self.material.G_s * A_l * phi ** 2 * K2 / (12 * self.kappa * L * (1 + phi) ** 2)
+            K = K + K2
+        # fmt: on
+
+        return K
+
+    def C(self):
+        """Stiffness matrix for an instance of a shaft tapered element.
+
+        Returns
+        -------
+        C : np.array
+           Damping matrix for the shaft tapered element.
+
+        Examples
+        --------
+        """
+        C = np.zeros((8, 8))
+
+        return C
+
+    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
+        >>> # Timoshenko is the default shaft element
+        >>> Timoshenko_Element = ShaftTaperedElement(
+        ...                         0.5, 0.05, 0.05, 0.1,
+        ...                         0.15, steel,
+        ...                         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.        ]])
+        """
+
+        G = np.zeros((8, 8))
+
+        if self.gyroscopic:
+            phi = self.phi
+            L = self.L
+            a1 = self.a1
+            a2 = self.a2
+            b1 = self.b1
+            b2 = self.b2
+            delta = self.delta
+            gama = self.gama
+            Ie_l = self.Ie_l
+
+            # fmt: off
+            g1 = 252 + 126 * a2 + 72 * b2 + 45 * gama + 30 * delta
+            g2 = (
+                21 - 105 * phi
+                + a2 * (21 - 42 * phi)
+                + b2 * (15 - 21 * phi)
+                + gama * (10.5 - 12 * phi)
+                + delta * (7.5 - 7.5 * phi)
+            )
+            g3 = (
+                21 - 105 * phi
+                - 63 * a2 * phi
+                - b2 * (6 + 42 * phi)
+                - gama * (7.5 + 30 * phi)
+                - delta * (7.5 + 22.5 * phi)
+            )
+            g4 = (
+                28 + 35 * phi + 70 * phi ** 2
+                + a2 * (7 - 7 * phi + 17.5 * phi ** 2)
+                + b2 * (4 - 7 * phi + 7 * phi ** 2)
+                + gama * (2.75 - 5 * phi + 3.5 * phi ** 2)
+                + delta * (2 - 3.5 * phi + 2 * phi ** 2)
+            )
+            g5 = (
+                7 + 35 * phi - 35 * phi ** 2
+                + a2 * (3.5 + 17.5 * phi - 17.5 * phi ** 2)
+                + b2 * (3 + 10.5 * phi - 10.5 * phi ** 2)
+                + gama * (2.75 + 7 * phi - 7 * phi ** 2)
+                + delta * (2.5 + 5 * phi - 5 * phi ** 2)
+            )
+            g6 = (
+                28 + 35 * phi + 70 * phi ** 2
+                + a2 * (21 + 42 * phi + 52.5 * phi ** 2)
+                + b2 * (18 + 42 * phi + 42 * phi ** 2)
+                + gama * (16.25 + 40 * phi + 35 * phi ** 2)
+                + delta * (15 + 37.5 * phi + 30 * phi ** 2)
+            )
+
+            G = np.array([
+                    [    0,    g1,     L*g2,        0,     0,  - g1,     L*g3,        0],
+                    [   g1,     0,        0,    -L*g2,   -g1,     0,        0,    -L*g3],
+                    [-L*g2,     0,        0,  L**2*g4,  L*g2,     0,        0, -L**2*g5],
+                    [    0,  L*g2,  L**2*g4,        0,     0, -L*g2, -L**2*g5,        0],
+                    [    0,   -g1,    -L*g2,        0,     0,    g2,    -L*g3,        0],
+                    [  -g1,     0,        0,     L*g2,    g2,     0,        0,    -L*g3],
+                    [-L*g3,     0,        0, -L**2*g5,  L*g3,     0,        0,  L**2*g6],
+                    [    0,  L*g3, -L**2*g5,        0,     0, -L*g3,  L**2*g6,        0],
+            ])
+            # fmt: on
+            G = self.material.rho * Ie_l * 2 * G / (210 * L * (1 + phi) ** 2)
+
+        return G
+
+    def patch(self, position, SR, ax):
+        """Shaft element patch.
+        Patch that will be used to draw the shaft element.
+        Parameters
+        ----------
+        ax : matplotlib axes, optional
+            Axes in which the plot will be drawn.
+        bk_ax : bokeh plotting axes, optional
+            Axes in which the plot will be drawn.
+        position : float
+            Position in which the patch will be drawn.
+        Returns
+        -------
+        """
+
+        if self.n in SR:
+            mpl_color = "yellow"
+            legend = "Shaft - Slenderness Ratio < 1.6"
+        else:
+            mpl_color = self.color
+            legend = "Shaft"
+
+        shaft_points_u = [
+            [position, self.i_d_l / 2],
+            [position, self.o_d_l / 2],
+            [position + self.L, self.o_d_r / 2],
+            [position + self.L, self.i_d_r / 2],
+            [position, self.i_d_l / 2],
+        ]
+        shaft_points_l = [
+            [position, -self.i_d_l / 2],
+            [position, -self.o_d_l / 2],
+            [position + self.L, -self.o_d_r / 2],
+            [position + self.L, -self.i_d_r / 2],
+            [position, -self.i_d_l / 2],
+        ]
+
+        # matplotlib - plot the upper half of the shaft
+        ax.add_patch(
+            mpatches.Polygon(
+                shaft_points_u,
+                facecolor=mpl_color,
+                linestyle="--",
+                linewidth=0.5,
+                ec="k",
+                alpha=0.8,
+                label=legend,
+            )
+        )
+        # matplotlib - plot the lower half of the shaft
+        ax.add_patch(
+            mpatches.Polygon(
+                shaft_points_l,
+                facecolor=mpl_color,
+                linestyle="--",
+                linewidth=0.5,
+                ec="k",
+                alpha=0.8,
+            )
+        )
+
+    def bokeh_patch(self, position, SR, bk_ax):
+        """Shaft element patch.
+        Patch that will be used to draw the shaft element.
+        Parameters
+        ----------
+        bk_ax : bokeh plotting axes, optional
+            Axes in which the plot will be drawn.
+        position : float
+            Position in which the patch will be drawn.
+        Returns
+        -------
+        """
+        if self.n in SR:
+            bk_color = "yellow"
+            legend = "Shaft - Slenderness Ratio < 1.6"
+        else:
+            bk_color = bokeh_colors[2]
+            legend = "Shaft"
+
+        # bokeh plot - plot the shaft
+        z_upper = [position, position, position + self.L, position + self.L]
+        y_upper = [self.i_d_l / 2, self.o_d_l / 2, self.o_d_r / 2, self.i_d_r / 2]
+        z_lower = [position, position, position + self.L, position + self.L]
+        y_lower = [-self.i_d_l / 2, -self.o_d_l / 2, -self.o_d_r / 2, -self.i_d_r / 2]
+
+        source = ColumnDataSource(
+            dict(
+                z_l=z_lower,
+                y_l=y_lower,
+                z_u=z_upper,
+                y_u=y_upper,
+                elnum=[self.n]*4,
+                out_d_l=[self.o_d_l]*4,
+                out_d_r=[self.o_d_r]*4,
+                in_d_l=[self.i_d_l]*4,
+                in_d_r=[self.i_d_r]*4,
+                length=[self.L]*4,
+                mat=[self.material.name]*4,
+            )
+        )
+
+        bk_ax.patch(
+            x="z_u",
+            y="y_u",
+            source=source,
+            line_color=bokeh_colors[0],
+            line_width=1,
+            fill_alpha=0.5,
+            fill_color=bk_color,
+            legend=legend,
+            name="u_shaft",
+        )
+        bk_ax.patch(
+            x="z_l",
+            y="y_l",
+            source=source,
+            line_color=bokeh_colors[0],
+            line_width=1,
+            fill_alpha=0.5,
+            fill_color=bk_color,
+            legend=legend,
+            name="l_shaft",
+        )
+
+        hover = HoverTool(names=["u_shaft", "l_shaft"])
+        hover.tooltips = [
+            ("Element Number :", "@elnum"),
+            ("Left Outer Diameter :", "@out_d_l"),
+            ("Left Inner Diameter :", "@in_d_l"),
+            ("Right Outer Diameter :", "@out_d_r"),
+            ("Right Inner Diameter :", "@in_d_r"),
+            ("Element Length :", "@length"),
+            ("Material :", "@mat"),
+        ]
+        hover.mode = "mouse"
+
+        return hover


Test patch:
diff --git a/ross/tests/test_shaft_element.py b/ross/tests/test_shaft_element.py
index f9617f11..8f187412 100644
--- a/ross/tests/test_shaft_element.py
+++ b/ross/tests/test_shaft_element.py
@@ -5,7 +5,7 @@ import pytest
 from numpy.testing import assert_almost_equal, assert_allclose
 
 from ross.materials import steel
-from ross.shaft_element import ShaftElement
+from ross.shaft_element import ShaftElement, ShaftTaperedElement
 
 test_dir = os.path.dirname(__file__)
 
@@ -151,3 +151,136 @@ 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_eb():
+    #  Euler-Bernoulli element
+    L = 0.4
+    i_d_l = 0.
+    i_d_r = 0.
+    o_d_l = 0.25
+    o_d_r = 0.10
+
+    return ShaftTaperedElement(
+        L, i_d_l, i_d_r, o_d_l, o_d_r, steel, shear_effects=False, rotary_inertia=False, n=3
+    )
+
+
+def test_tapered_index(tap_eb):
+    assert tap_eb.dof_local_index().x0 == 0
+    assert tap_eb.dof_local_index().y0 == 1
+    assert tap_eb.dof_local_index().alpha0 == 2
+    assert tap_eb.dof_local_index().beta0 == 3
+    assert tap_eb.dof_global_index().x0 == 12
+    assert tap_eb.dof_global_index().y0 == 13
+    assert tap_eb.dof_global_index().alpha0 == 14
+    assert tap_eb.dof_global_index().beta0 == 15
+
+
+def test_parameters_tap(tap_eb):
+    assert tap_eb.phi == 0
+    assert tap_eb.L == 0.4
+    assert tap_eb.i_d == 0.0
+    assert tap_eb.i_d_l == 0.0
+    assert tap_eb.i_d_r == 0.0
+    assert tap_eb.o_d == 0.175
+    assert tap_eb.o_d_l == 0.25
+    assert tap_eb.o_d_r == 0.1
+    assert tap_eb.material.E == 211e9
+    assert tap_eb.material.G_s == 81.2e9
+    assert tap_eb.material.rho == 7810
+    assert_almost_equal(tap_eb.material.Poisson, 0.29926108)
+    assert_almost_equal(tap_eb.A, 0.024052818754)
+    assert_almost_equal(tap_eb.Ie * 1e5, 4.44854429)
+
+
+def test_mass_matrix_tap_eb(tap_eb):
+    # fmt: off
+    M0e_tap = np.array([[42.8500896 ,  0.        ,  0.        ,  2.13520215,  9.90196344,         0.        ,  0.        , -1.01356381],
+                        [ 0.        , 42.8500896 , -2.13520215,  0.        ,  0.        ,         9.90196344,  1.01356381,  0.        ],
+                        [ 0.        , -2.13520215,  0.14254154,  0.        ,  0.        ,        -0.89088462, -0.087628  ,  0.        ],
+                        [ 2.13520215,  0.        ,  0.        ,  0.14254154,  0.89088462,         0.        ,  0.        , -0.087628  ],
+                        [ 9.90196344,  0.        ,  0.        ,  0.89088462, 17.08745904,         0.        ,  0.        , -1.1537686 ],
+                        [ 0.        ,  9.90196344, -0.89088462,  0.        ,  0.        ,        17.08745904,  1.1537686 ,  0.        ],
+                        [ 0.        ,  1.01356381, -0.087628  ,  0.        ,  0.        ,         1.1537686 ,  0.09346986,  0.        ],
+                        [-1.01356381,  0.        ,  0.        , -0.087628  , -1.1537686 ,         0.        ,  0.        ,  0.09346986]])
+    # fmt: on
+    assert_allclose(tap_eb.M(), M0e_tap, rtol=1e-3)
+
+
+def test_stiffness_matrix_tap_eb(tap_eb):
+    # fmt: off
+    K0e_tap = np.array([[-484.29115593,    0.        ,    0.        ,   80.59381675,         484.29115593,    0.        ,    0.        ,   29.13029549],
+                        [   0.        , -484.29115593,  -80.59381675,    0.        ,           0.        ,  484.29115593,  -29.13029549,    0.        ],
+                        [   0.        ,  -80.59381675,   24.33997997,    0.        ,           0.        ,   80.59381675,    7.89754672,    0.        ],
+                        [  80.59381675,    0.        ,    0.        ,   24.33997997,         -80.59381675,    0.        ,    0.        ,    7.89754672],
+                        [ 484.29115593,    0.        ,    0.        ,  -80.59381675,        -484.29115593,    0.        ,    0.        ,  -29.13029549],
+                        [   0.        ,  484.29115593,   80.59381675,    0.        ,           0.        , -484.29115593,   29.13029549,    0.        ],
+                        [   0.        ,  -29.13029549,    7.89754672,    0.        ,           0.        ,   29.13029549,    3.75457147,    0.        ],
+                        [  29.13029549,    0.        ,    0.        ,    7.89754672,         -29.13029549,    0.        ,    0.        ,    3.75457147]])
+    # fmt: on
+    assert_almost_equal(tap_eb.K() / 1e7, K0e_tap, decimal=5)
+
+
+@pytest.fixture
+def tap_tim():
+    #  Timoshenko element
+    L = 0.4
+    i_d_l = 0.
+    i_d_r = 0.
+    o_d_l = 0.25
+    o_d_r = 0.10
+
+    return ShaftTaperedElement(
+        L, i_d_l, i_d_r, o_d_l, o_d_r, steel, shear_effects=True, rotary_inertia=True
+    )
+
+
+def test_parameters_tap_tim(tap_tim):
+    assert_almost_equal(tap_tim.phi, 0.40668276)
+    assert_almost_equal(tap_tim.material.Poisson, 0.29926108)
+    assert_almost_equal(tap_tim.A, 0.02405282)
+    assert_almost_equal(tap_tim.Ie * 1e5, 4.44854429057887)
+
+
+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, 16.54087906,         0.        ,  0.        , -1.14296052],
+                        [ 0.        , 10.27411395, -1.04155143,  0.        ,  0.        ,        16.54087906,  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]])
+    # fmt: on
+    assert_almost_equal(tap_tim.M(), M0e_tim, decimal=5)
+
+
+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]])
+    # 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.        ,        -0.39661891, -0.05013491,  0.        ],
+                        [-1.19360236,  0.        ,  0.        , -0.15864756, -0.39661891,         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.        ]])
+    # fmt: on
+    assert_almost_equal(tap_tim.G(), G0e_tim, decimal=5)