Repository: ross-rotordynamics/ross
Pull request #455 (2020-02-25 15:25:21)

Adapt static analysis for CoAxialRotor
Related to issues #443 and #454.

Some improvements to `.run_static()`:
- Supports CoAxialRotor class (multiple shafts)
- More information when returning disks and bearing forces. The user can choose to return a dict with nodal position for each force or with a element tag.
- Remove some repetitive calculus on code.

Some improvements to `StaticResults`:
- `StaticResults` now supports plotting results for multiple shafts.
- Multiple shafts are plotted in the same figure (`.plot_free_body_diagram()` will display multiple figures).
- New: the user can ommit a result by clicking on the respective legend (not valid for `.plot_free_body_diagram()`).


Patch:
diff --git a/ross/api_report.py b/ross/api_report.py
index e97168cd..96be87b6 100644
--- a/ross/api_report.py
+++ b/ross/api_report.py
@@ -211,7 +211,7 @@ class Report:
         """
         # get reaction forces on bearings
         self.rotor.run_static()
-        Fb = self.rotor.bearing_reaction_forces
+        Fb = list(self.rotor.bearing_forces_nodal.values())
         Fb = np.array(Fb) / 9.8065
 
         return Fb
diff --git a/ross/results.py b/ross/results.py
index eabeeffb..ec433e5a 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -1697,15 +1697,15 @@ class StaticResults:
     """
 
     def __init__(
-        self, disp_y, Vx, Bm, df_shaft, df_disks, df_bearings, nodes, nodes_pos, Vx_axis
+        self, disp_y, Vx, Bm, w_shaft, disk_forces, bearing_forces, nodes, nodes_pos, Vx_axis
     ):
 
         self.disp_y = disp_y
         self.Vx = Vx
         self.Bm = Bm
-        self.df_shaft = df_shaft
-        self.df_disks = df_disks
-        self.df_bearings = df_bearings
+        self.w_shaft = w_shaft
+        self.disk_forces = disk_forces
+        self.bearing_forces = bearing_forces
         self.nodes = nodes
         self.nodes_pos = nodes_pos
         self.Vx_axis = Vx_axis
@@ -1724,76 +1724,70 @@ class StaticResults:
         fig : bokeh figure
             Bokeh figure with static deformation plot
         """
-        source = ColumnDataSource(
-            data=dict(x=self.nodes_pos, y0=self.disp_y, y1=[0] * len(self.nodes_pos))
-        )
-
-        TOOLTIPS = [
-            ("Shaft lenght:", "@x"),
-            ("Underformed:", "@y1"),
-            ("Displacement:", "@y0"),
-        ]
-
         # create displacement plot
         fig = figure(
-            tools="pan, wheel_zoom, box_zoom, reset, save, box_select, hover",
-            tooltips=TOOLTIPS,
+            tools="pan, wheel_zoom, box_zoom, reset, save, box_select",
             width=640,
             height=480,
-            title="Static Analysis",
+            title="Deformation",
             x_axis_label="Shaft lenght",
             y_axis_label="Lateral displacement",
         )
-        fig.xaxis.axis_label_text_font_size = "20pt"
-        fig.yaxis.axis_label_text_font_size = "20pt"
+        fig.xaxis.axis_label_text_font_size = "16pt"
+        fig.yaxis.axis_label_text_font_size = "16pt"
         fig.axis.major_label_text_font_size = "16pt"
         fig.title.text_font_size = "14pt"
 
-        interpolated = interpolate.interp1d(
-            source.data["x"], source.data["y0"], kind="cubic"
-        )
-        xnew = np.linspace(
-            source.data["x"][0],
-            source.data["x"][-1],
-            num=len(self.nodes_pos) * 20,
-            endpoint=True,
-        )
+        count = 0
+        for disp_y, Vx, Bm, nodes, nodes_pos, Vx_axis in zip(
+                self.disp_y, self.Vx, self.Bm, self.nodes, self.nodes_pos, self.Vx_axis
+        ):
+            source = ColumnDataSource(
+                data=dict(x=nodes_pos, y=disp_y)
+            )
 
-        ynew = interpolated(xnew)
-        auxsource = ColumnDataSource(data=dict(x=xnew, y0=ynew, y1=[0] * len(xnew)))
+            interpolated = interpolate.interp1d(
+                source.data["x"], source.data["y"], kind="cubic"
+            )
+            xnew = np.linspace(
+                source.data["x"][0],
+                source.data["x"][-1],
+                num=len(nodes_pos) * 20,
+                endpoint=True,
+            )
 
-        fig.line(
-            "x",
-            "y0",
-            source=auxsource,
-            legend="Deformed shaft",
-            line_width=3,
-            line_color=bokeh_colors[9],
-        )
-        fig.circle(
-            "x",
-            "y0",
-            source=source,
-            legend="Deformed shaft",
-            size=8,
-            fill_color=bokeh_colors[9],
-        )
-        fig.line(
-            "x",
-            "y1",
-            source=source,
-            legend="underformed shaft",
-            line_width=3,
-            line_color=bokeh_colors[0],
-        )
-        fig.circle(
-            "x",
-            "y1",
-            source=source,
-            legend="underformed shaft",
-            size=8,
-            fill_color=bokeh_colors[0],
-        )
+            ynew = interpolated(xnew)
+            auxsource = ColumnDataSource(data=dict(x=xnew, y=ynew))
+
+            fig.line(
+                "x",
+                "y",
+                source=auxsource,
+                legend_label="Deformed - shaft " + str(count),
+                line_width=3,
+                line_color=bokeh_colors[9 - count],
+                muted_alpha=0.1,
+                name="def_l",
+            )
+            fig.circle(
+                "x",
+                "y",
+                source=source,
+                legend_label="Deformed - shaft " + str(count),
+                size=8,
+                fill_color=bokeh_colors[9 - count],
+                muted_alpha=0.1,
+                name="def_c",
+            )
+            hover = HoverTool(names=["def_l", "def_c"])
+            hover.tooltips = [
+                ("Shaft lenght:", "@x"),
+                ("Displacement:", "@y"),
+            ]
+            count += 1
+        fig.add_tools(hover)
+        fig.legend.background_fill_alpha = 0.1
+        fig.legend.click_policy = "mute"
 
         return fig
 
@@ -1811,142 +1805,148 @@ class StaticResults:
         fig : bokeh figure
             Bokeh figure with the free-body diagram plot
         """
-        source = ColumnDataSource(
-            data=dict(
-                x=self.nodes_pos, y0=self.disp_y * 1000, y1=[0] * len(self.nodes_pos)
-            )
-        )
-
+        figures = []
+        j = 0
         y_start = 5.0
-        sh_weight = sum(self.df_shaft["m"].values) * 9.8065
-
-        shaft_end = self.nodes_pos[-1]
-        fig = figure(
-            tools="pan, wheel_zoom, box_zoom, reset, save, box_select, hover",
-            width=640,
-            height=480,
-            title="Free-Body Diagram",
-            x_axis_label="shaft lenght",
-            x_range=[-0.1 * shaft_end, 1.1 * shaft_end],
-            y_range=[-3 * y_start, 3 * y_start],
-        )
-        fig.yaxis.visible = False
-        fig.xaxis.axis_label_text_font_size = "20pt"
-        fig.axis.major_label_text_font_size = "16pt"
-        fig.title.text_font_size = "14pt"
-
-        fig.line("x", "y1", source=source, line_width=5, line_color=bokeh_colors[0])
-
-        # fig - plot arrows indicating shaft weight distribution
-        text = str("%.1f" % sh_weight)
-        fig.line(
-            x=self.nodes_pos,
-            y=[y_start] * len(self.nodes_pos),
-            line_width=2,
-            line_color=bokeh_colors[0],
-        )
-
-        ini = self.nodes_pos[0]
-        fin = self.nodes_pos[-1]
-        arrows_list = np.arange(ini, 1.01 * fin, fin / 5.0)
-        for node in arrows_list:
-            fig.add_layout(
-                Arrow(
-                    end=NormalHead(
-                        fill_color=bokeh_colors[2],
-                        fill_alpha=1.0,
-                        size=16,
-                        line_width=2,
-                        line_color=bokeh_colors[0],
-                    ),
-                    x_start=node,
-                    y_start=y_start,
-                    x_end=node,
-                    y_end=0,
-                )
+        for nodes_pos, nodes in zip(self.nodes_pos, self.nodes):
+            fig = figure(
+                tools="pan, wheel_zoom, box_zoom, reset, save, box_select",
+                width=640,
+                height=480,
+                title="Free-Body Diagram - Shaft " + str(j),
+                x_axis_label="Shaft lenght",
+                x_range=[-0.1 * nodes_pos[-1], 1.1 * nodes_pos[-1]],
+                y_range=[-3 * y_start, 3 * y_start],
             )
+            fig.yaxis.visible = False
+            fig.xaxis.axis_label_text_font_size = "16pt"
+            fig.axis.major_label_text_font_size = "16pt"
+            fig.title.text_font_size = "14pt"
 
-        fig.add_layout(
-            Label(
-                x=self.nodes_pos[0],
-                y=y_start,
-                text="W = " + text + "N",
-                text_font_style="bold",
-                text_font_size="10pt",
-                text_baseline="top",
-                text_align="left",
-                y_offset=20,
+            fig.line(
+                nodes_pos,
+                np.zeros(len(nodes_pos)),
+                line_width=5,
+                color=bokeh_colors[0]
             )
-        )
 
-        # fig - calculate the reaction force of bearings and plot arrows
-        for i, node in enumerate(self.df_bearings["n"]):
-            Fb = -self.disp_y[node] * self.df_bearings.loc[i, "kyy"].coefficient[0]
-            text = str("%.1f" % Fb)
-            fig.add_layout(
-                Arrow(
-                    end=NormalHead(
-                        fill_color=bokeh_colors[6],
-                        fill_alpha=1.0,
-                        size=16,
-                        line_width=2,
-                        line_color=bokeh_colors[0],
-                    ),
-                    x_start=self.nodes_pos[node],
-                    y_start=-2 * y_start,
-                    x_end=self.nodes_pos[node],
-                    y_end=0,
-                )
-            )
-            fig.add_layout(
-                Label(
-                    x=self.nodes_pos[node],
-                    y=-2 * y_start,
-                    angle=np.pi / 2,
-                    text="Fb = " + text + "N",
-                    text_font_style="bold",
-                    text_font_size="10pt",
-                    text_baseline="top",
-                    text_align="center",
-                    x_offset=2,
-                )
+            # fig - plot arrows indicating shaft weight distribution
+            text = str("%.1f" % self.w_shaft[j])
+            fig.line(
+                x=nodes_pos,
+                y=[y_start] * len(nodes_pos),
+                line_width=2,
+                line_color=bokeh_colors[0],
             )
 
-        # fig - plot arrows indicating disk weight
-        if len(self.df_disks) != 0:
-            for i, node in enumerate(self.df_disks["n"]):
-                Fd = self.df_disks.loc[i, "m"] * 9.8065
-                text = str("%.1f" % Fd)
+            ini = nodes_pos[0]
+            fin = nodes_pos[-1]
+            arrows_list = np.arange(ini, 1.01 * fin, (fin - ini) / 5.0)
+            for node in arrows_list:
                 fig.add_layout(
                     Arrow(
                         end=NormalHead(
-                            fill_color=bokeh_colors[9],
+                            fill_color=bokeh_colors[2],
                             fill_alpha=1.0,
-                            size=16,
+                            size=14,
                             line_width=2,
                             line_color=bokeh_colors[0],
                         ),
-                        x_start=self.nodes_pos[node],
-                        y_start=2 * y_start,
-                        x_end=self.nodes_pos[node],
+                        x_start=node,
+                        y_start=y_start,
+                        x_end=node,
                         y_end=0,
                     )
                 )
-                fig.add_layout(
-                    Label(
-                        x=self.nodes_pos[node],
-                        y=2 * y_start,
-                        angle=np.pi / 2,
-                        text="Fd = " + text + "N",
-                        text_font_style="bold",
-                        text_font_size="10pt",
-                        text_baseline="top",
-                        text_align="center",
-                        x_offset=2,
-                    )
+
+            fig.add_layout(
+                Label(
+                    x=nodes_pos[0],
+                    y=y_start,
+                    text="Weight = " + text + "N",
+                    text_font_style="bold",
+                    text_font_size="10pt",
+                    text_baseline="top",
+                    text_align="left",
+                    y_offset=20,
                 )
+            )
 
-        return fig
+            # fig - calculate the reaction force of bearings and plot arrows
+            for k, v in self.bearing_forces.items():
+                _, node = k.split("_")
+                node = int(node)
+                if node in nodes:
+                    text = str(v)
+                    var = 1 if v < 0 else - 1
+                    fig.add_layout(
+                        Arrow(
+                            end=NormalHead(
+                                fill_color=bokeh_colors[6],
+                                fill_alpha=1.0,
+                                size=14,
+                                line_width=2,
+                                line_color=bokeh_colors[0],
+                            ),
+                            x_start=nodes_pos[nodes.index(node)],
+                            y_start=var * 2 * y_start,
+                            x_end=nodes_pos[nodes.index(node)],
+                            y_end=0,
+                        )
+                    )
+                    fig.add_layout(
+                        Label(
+                            x=nodes_pos[nodes.index(node)],
+                            y=var * 2 * y_start,
+                            angle=np.pi / 2,
+                            text="Fb = " + text + "N",
+                            text_font_style="bold",
+                            text_font_size="10pt",
+                            text_baseline="top",
+                            text_align="center",
+                            x_offset=2,
+                        )
+                    )
+
+            # fig - plot arrows indicating disk weight
+            for k, v in self.disk_forces.items():
+                _, node = k.split("_")
+                node = int(node)
+                if node in nodes:
+                    text = str(v)
+                    fig.add_layout(
+                        Arrow(
+                            end=NormalHead(
+                                fill_color=bokeh_colors[9],
+                                fill_alpha=1.0,
+                                    size=14,
+                                line_width=2,
+                                line_color=bokeh_colors[0],
+                            ),
+                            x_start=nodes_pos[nodes.index(node)],
+                            y_start=2 * y_start,
+                            x_end=nodes_pos[nodes.index(node)],
+                            y_end=0,
+                        )
+                    )
+                    fig.add_layout(
+                        Label(
+                            x=nodes_pos[nodes.index(node)],
+                            y=2 * y_start,
+                            angle=np.pi / 2,
+                            text="Fd = " + text + "N",
+                            text_font_style="bold",
+                            text_font_size="10pt",
+                            text_baseline="top",
+                            text_align="center",
+                            x_offset=2,
+                        )
+                    )
+            figures.append(fig)
+            j += 1
+        grid_plots = gridplot([figures])
+
+        return grid_plots
 
     def plot_shearing_force(self):
         """Plot the rotor shearing force diagram.
@@ -1962,12 +1962,9 @@ class StaticResults:
         fig : bokeh figure
             Bokeh figure with the shearing force diagram plot
         """
-        shaft_end = self.nodes_pos[-1]
-        source_SF = ColumnDataSource(data=dict(x=self.Vx_axis, y=self.Vx))
-        TOOLTIPS_SF = [("Shearing Force:", "@y")]
+        shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
         fig = figure(
-            tools="pan, wheel_zoom, box_zoom, reset, save, box_select, hover",
-            tooltips=TOOLTIPS_SF,
+            tools="pan, wheel_zoom, box_zoom, reset, save, box_select",
             width=640,
             height=480,
             title="Shearing Force Diagram",
@@ -1975,14 +1972,11 @@ class StaticResults:
             y_axis_label="Force",
             x_range=[-0.1 * shaft_end, 1.1 * shaft_end],
         )
-        fig.xaxis.axis_label_text_font_size = "20pt"
-        fig.yaxis.axis_label_text_font_size = "20pt"
+        fig.xaxis.axis_label_text_font_size = "16pt"
+        fig.yaxis.axis_label_text_font_size = "16pt"
         fig.axis.major_label_text_font_size = "16pt"
         fig.title.text_font_size = "14pt"
 
-        fig.line("x", "y", source=source_SF, line_width=4, line_color=bokeh_colors[0])
-        fig.circle("x", "y", source=source_SF, size=8, fill_color=bokeh_colors[0])
-
         # fig - plot centerline
         fig.line(
             [-0.1 * shaft_end, 1.1 * shaft_end],
@@ -1992,6 +1986,39 @@ class StaticResults:
             line_color=bokeh_colors[0],
         )
 
+        j = 0
+        for Vx, Vx_axis in zip(self.Vx, self.Vx_axis):
+            source_SF = ColumnDataSource(data=dict(x=Vx_axis, y=Vx))
+
+            fig.line(
+                "x",
+                "y",
+                source=source_SF,
+                line_width=4,
+                line_color=bokeh_colors[9 - j],
+                line_alpha=1.0,
+                muted_alpha=0.1,
+                legend_label="Shaft " + str(j),
+                name="line",
+            )
+            fig.circle(
+                "x",
+                "y",
+                source=source_SF,
+                size=8,
+                fill_color=bokeh_colors[9 - j],
+                fill_alpha=1.0,
+                muted_alpha=0.1,
+                legend_label="Shaft " + str(j),
+                name="circle",
+            )
+            hover = HoverTool(names=["line", "circle"])
+            hover.tooltips = [("Shear Force:", "@y")]
+            j += 1
+        fig.add_tools(hover)
+        fig.legend.background_fill_alpha = 0.1
+        fig.legend.click_policy = "mute"
+
         return fig
 
     def plot_bending_moment(self):
@@ -2008,12 +2035,9 @@ class StaticResults:
         fig : bokeh figure
             Bokeh figure with the bending moment diagram plot
         """
-        shaft_end = self.nodes_pos[-1]
-        source_BM = ColumnDataSource(data=dict(x=self.nodes_pos, y=self.Bm))
-        TOOLTIPS_BM = [("Bending Moment:", "@y")]
+        shaft_end = max([sublist[-1] for sublist in self.nodes_pos])
         fig = figure(
-            tools="pan, wheel_zoom, box_zoom, reset, save, box_select, hover",
-            tooltips=TOOLTIPS_BM,
+            tools="pan, wheel_zoom, box_zoom, reset, save, box_select",
             width=640,
             height=480,
             title="Bending Moment Diagram",
@@ -2026,26 +2050,6 @@ class StaticResults:
         fig.axis.major_label_text_font_size = "16pt"
         fig.title.text_font_size = "14pt"
 
-        i = 0
-        while True:
-            if i + 3 > len(self.nodes):
-                break
-
-            interpolated_BM = interpolate.interp1d(
-                self.nodes_pos[i : i + 3], self.Bm[i : i + 3], kind="quadratic"
-            )
-            xnew_BM = np.linspace(
-                self.nodes_pos[i], self.nodes_pos[i + 2], num=42, endpoint=True
-            )
-
-            ynew_BM = interpolated_BM(xnew_BM)
-            auxsource_BM = ColumnDataSource(data=dict(x=xnew_BM, y=ynew_BM))
-            fig.line(
-                "x", "y", source=auxsource_BM, line_width=4, line_color=bokeh_colors[0]
-            )
-            i += 2
-        fig.circle("x", "y", source=source_BM, size=8, fill_color=bokeh_colors[0])
-
         # fig - plot centerline
         fig.line(
             [-0.1 * shaft_end, 1.1 * shaft_end],
@@ -2055,6 +2059,55 @@ class StaticResults:
             line_color=bokeh_colors[0],
         )
 
+        j = 0
+        for Bm, nodes_pos, nodes in zip(self.Bm, self.nodes_pos, self.nodes):
+            source_BM = ColumnDataSource(data=dict(x=nodes_pos, y=Bm))
+            i = 0
+            while True:
+                if i + 3 > len(nodes):
+                    break
+
+                interpolated_BM = interpolate.interp1d(
+                    nodes_pos[i : i + 3], Bm[i : i + 3], kind="quadratic"
+                )
+                xnew_BM = np.linspace(
+                    nodes_pos[i], nodes_pos[i + 2], num=42, endpoint=True
+                )
+
+                ynew_BM = interpolated_BM(xnew_BM)
+                auxsource_BM = ColumnDataSource(data=dict(x=xnew_BM, y=ynew_BM))
+                fig.line(
+                    "x",
+                    "y",
+                    source=auxsource_BM,
+                    line_width=4,
+                    line_color=bokeh_colors[9 - j],
+                    line_alpha=1.0,
+                    muted_alpha=0.1,
+                    legend_label="Shaft " + str(j),
+                    name="line",
+                )
+                i += 2
+
+            fig.circle(
+                "x",
+                "y",
+                source=source_BM,
+                size=8,
+                fill_color=bokeh_colors[9 - j],
+                fill_alpha=1.0,
+                muted_alpha=0.1,
+                legend_label="Shaft " + str(j),
+                name="circle",
+            )
+            hover = HoverTool(names=["line", "circle"])
+            hover.tooltips = [("Beinding Moment:", "@y")]
+            j += 1
+
+        fig.add_tools(hover)
+        fig.legend.background_fill_alpha = 0.1
+        fig.legend.click_policy = "mute"
+
         return fig
 
 
@@ -2119,6 +2172,7 @@ class SummaryResults:
 
         shaft_data = dict(
             tags=self.df_shaft["tag"],
+            sh_number=self.df_shaft["shaft_number"],
             lft_stn=self.df_shaft["n_l"],
             rgt_stn=self.df_shaft["n_r"],
             elem_no=self.df_shaft["_n"],
@@ -2145,6 +2199,7 @@ class SummaryResults:
 
         disk_data = dict(
             tags=self.df_disks["tag"],
+            sh_number=self.df_disks["shaft_number"],
             disk_node=self.df_disks["n"],
             disk_pos=self.nodes_pos[self.df_bearings["n"]],
             disk_mass=self.df_disks["m"],
@@ -2153,9 +2208,10 @@ class SummaryResults:
 
         bearing_data = dict(
             tags=self.df_bearings["tag"],
+            sh_number=self.df_bearings["shaft_number"],
             brg_node=self.df_bearings["n"],
             brg_pos=self.nodes_pos[self.df_bearings["n"]],
-            brg_force=self.brg_forces,
+            brg_force=list(self.brg_forces.values()),
         )
 
         shaft_source = ColumnDataSource(shaft_data)
@@ -2165,42 +2221,45 @@ class SummaryResults:
 
         shaft_titles = [
             "Element Tag",
+            "Shaft Number",
             "Left Station",
             "Right Station",
             "Element Number",
-            "Elem. Left Location (m)",
-            "Elem. Lenght (m)",
-            "Element CG (m)",
-            "Axial CG Location (m)",
-            "Elem. Right Location (m)",
+            "Elem. Left Location",
+            "Elem. Lenght",
+            "Element CG",
+            "Axial CG Location",
+            "Elem. Right Location",
             "Material",
-            "Elem. Mass (kg)",
-            "Inertia (kg.m²)",
+            "Elem. Mass",
+            "Inertia",
         ]
 
         rotor_titles = [
             "Tag",
             "First Station",
             "Last Station",
-            "Starting Pos. (m)",
-            "Total Lenght (m)",
-            "C.G. Locantion (m)",
-            "Total Ip about C.L. (kg.m²)",
+            "Starting Pos.",
+            "Total Lenght",
+            "C.G. Locantion",
+            "Total Ip about C.L.",
         ]
 
         disk_titles = [
             "Tag",
+            "Shaft Number",
             "Disk Station",
-            "C.G. Locantion (m)",
-            "Disk Mass (m)",
-            "Total Ip about C.L. (kg.m²)",
+            "C.G. Locantion",
+            "Disk Mass",
+            "Total Ip about C.L.",
         ]
 
         bearing_titles = [
             "Tag",
+            "Shaft Number",
             "Bearing Station",
-            "Bearing Locantion (m)",
-            "Static Reaction Force (N)",
+            "Bearing Locantion",
+            "Static Reaction Force",
         ]
 
         shaft_formatters = [
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index c47e261a..dbb7fc11 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -2178,141 +2178,270 @@ class Rotor(object):
         -------
         >>> rotor = rotor_example()
         >>> static = rotor.run_static()
-        >>> rotor.bearing_reaction_forces
-        [432.4, 432.4]
+        >>> rotor.bearing_forces_nodal
+        {'node_0': 432.4, 'node_6': 432.4}
+        >>> rotor.bearing_forces_tag
+        {'Bearing 0': 432.4, 'Bearing 1': 432.4}
         """
-        if len(self.df_bearings) == 0:
+        if not len(self.df_bearings):
             raise ValueError("Rotor has no bearings")
 
-        else:
-            # gravity aceleration vector
-            grav = np.zeros((len(self.M()), 1))
+        aux_brg = []
+        for elm in self.bearing_elements:
+            if elm.n not in self.nodes:
+                pass
+            elif elm.n_link in self.nodes:
+                aux_brg.append(
+                    BearingElement(n=elm.n, n_link=elm.n_link, kxx=1e14, cxx=0)
+                )
+            else:
+                aux_brg.append(BearingElement(n=elm.n, kxx=1e14, cxx=0))
 
-            # place gravity effect on shaft and disks nodes
-            for node_y in range(int(len(self.M()) / 4)):
-                grav[4 * node_y + 1] = -9.8065
+        if isinstance(self, CoAxialRotor):
+            aux_rotor = CoAxialRotor(self.shafts, self.disk_elements, aux_brg)
+        else:
+            aux_rotor = Rotor(self.shaft_elements, self.disk_elements, aux_brg)
+        aux_K = aux_rotor.K(0)
+        aux_M = aux_rotor.M()
+
+        for elm in aux_rotor.bearing_elements:
+            if isinstance(elm, SealElement):
+                dofs = elm.dof_global_index
+                try:
+                    aux_K[np.ix_(dofs, dofs)] -= elm.K(0)
+                except TypeError:
+                    aux_K[np.ix_(dofs, dofs)] -= elm.K()
+
+        df_num = aux_rotor.df["shaft_number"].values
+        sh_num = [int(item) for item, count in Counter(df_num).items() if count > 1]
+
+        # gravity aceleration vector
+        g = 9.8065
+        grav = np.zeros(len(aux_K))
+
+        # place gravity effect on shaft and disks nodes
+        for node_y in range(int(len(aux_K) / 4)):
+            grav[4 * node_y + 1] = -g
+
+        # calculates x, for [K]*(x) = [M]*(g)
+        disp = (la.solve(aux_K, aux_M @ grav)).flatten()
+
+        # calculates displacement values in gravity's direction
+        # dof = degree of freedom
+        disp_y = np.array([])
+        for node_dof in range(int(len(disp) / 4)):
+            disp_y = np.append(disp_y, disp[4 * node_dof + 1])
+
+        # Shearing Force
+        BRG = [0] * len(self.nodes_pos)
+        DSK = [0] * len(self.nodes_pos)
+        SCH = [0] * len(self.nodes_pos)
+        BrgForce_nodal = {"node_" + str(i): 0 for i in self.nodes}
+        DskForce_nodal = {"node_" + str(i): 0 for i in self.nodes}
+        BrgForce_tag = {"node_" + str(i): 0 for i in self.nodes}
+        DskForce_tag = {"node_" + str(i): 0 for i in self.nodes}
+
+        # Bearing Forces
+        for i, node in enumerate(aux_rotor.df_bearings["n"]):
+            if not pd.isna(aux_rotor.df_bearings.loc[i, "n_link"]):
+                BRG[node] = (
+                    BRG[node]
+                    + disp_y[node]
+                    * self.df_bearings.loc[
+                        self.df_bearings.tag == aux_rotor.df_bearings.tag, "kyy"
+                    ][0].coefficient[0]
+                )
+                BrgForce_nodal["node_" + str(node)] = np.around(
+                    BrgForce_nodal["node_" + str(node)]
+                    + disp_y[node]
+                    * self.df_bearings.loc[
+                        self.df_bearings.tag == aux_rotor.df_bearings.tag, "kyy"
+                    ][0].coefficient[0],
+                    decimals=1,
+                )
+                BrgForce_tag[aux_rotor.df_bearings.loc[i, "tag"]] = BrgForce_nodal[
+                    "node_" + str(node)
+                ]
 
-            aux_brg = []
-            for n in self.df_bearings["n"]:
-                aux_brg.append(BearingElement(n=n, kxx=1e14, cxx=0))
+                node = int(aux_rotor.df_bearings.loc[i, "n_link"])
+                BRG[node] = (
+                    BRG[node]
+                    - disp_y[node]
+                    * self.df_bearings.loc[self.df_bearings.n_link == node, "kyy"]
+                    .values[0]
+                    .coefficient[0]
+                )
+                BrgForce_nodal["node_" + str(node)] = np.around(
+                    BrgForce_nodal["node_" + str(node)]
+                    - disp_y[node]
+                    * self.df_bearings.loc[self.df_bearings.n_link == node, "kyy"]
+                    .values[0]
+                    .coefficient[0],
+                    decimals=1,
+                )
+                BrgForce_tag[aux_rotor.df_bearings.loc[i, "tag"]] = BrgForce_nodal[
+                    "node_" + str(node)
+                ]
 
-            aux_rotor = Rotor(self.shaft_elements, self.disk_elements, aux_brg)
-            aux_K = np.zeros((aux_rotor.ndof, aux_rotor.ndof))
-
-            for elm in aux_rotor.elements:
-                if elm.__class__.__name__ != "SealElement":
-                    dofs = elm.dof_global_index
-                    n0 = dofs[0]
-                    n1 = dofs[-1] + 1  # +1 to include this dof in the slice
-                    try:
-                        aux_K[n0:n1, n0:n1] += elm.K(0)
-                    except TypeError:
-                        aux_K[n0:n1, n0:n1] += elm.K()
-
-            # calculates x, for [K]*(x) = [M]*(g)
-            disp = (la.solve(aux_K, self.M() @ grav)).flatten()
-
-            # calculates displacement values in gravity's direction
-            # dof = degree of freedom
-            disp_y = np.array([])
-            for node_dof in range(int(len(disp) / 4)):
-                disp_y = np.append(disp_y, disp[4 * node_dof + 1])
-            self.disp_y = disp_y
-
-            # Shearing Force
-            BrgForce = [0] * len(self.nodes_pos)
-            DskForce = [0] * len(self.nodes_pos)
-            SchForce = [0] * len(self.nodes_pos)
-
-            # Bearing Forces
-            BrgForceToReturn = []
-            for i, node in enumerate(self.df_bearings["n"]):
-                BrgForce[node] = (
-                    -disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0]
+            else:
+                BRG[node] = (
+                    BRG[node]
+                    - disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0]
                 )
-                BrgForceToReturn.append(
-                    np.around(
-                        -disp_y[node]
-                        * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0],
-                        decimals=1,
-                    )
+                BrgForce_nodal["node_" + str(node)] = np.around(
+                    BrgForce_nodal["node_" + str(node)]
+                    - disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0],
+                    decimals=1,
                 )
+                BrgForce_tag[aux_rotor.df_bearings.loc[i, "tag"]] = BrgForce_nodal[
+                    "node_" + str(node)
+                ]
 
-            # Disk Forces
-            DskForceToReturn = []
-            if len(self.df_disks) != 0:
-                for i, node in enumerate(self.df_disks["n"]):
-                    DskForce[node] = self.df_disks.loc[i, "m"] * -9.8065
-                    DskForceToReturn.append(
-                        np.around(self.df_disks.loc[i, "m"] * -9.8065, decimals=1)
-                    )
+        # counting nodes with more than 1 bearing attached to
+        node_b = list(aux_rotor.df_bearings["n"])
+        node_b.extend(list(aux_rotor.df_bearings["n_link"]))
+        count = len(node_b) - len(Counter(node_b))
+
+        # Disk Forces
+        if len(self.df_disks):
+            for i, node in enumerate(self.df_disks["n"]):
+                DSK[node] = self.df_disks.loc[i, "m"] * -g
+                DskForce_nodal["node_" + str(node)] = np.around(
+                    self.df_disks.loc[i, "m"] * -g, decimals=1
+                )
+                DskForce_tag[aux_rotor.df_disks.loc[i, "tag"]] = DskForce_nodal[
+                    "node_" + str(node)
+                ]
 
-            # Shaft Weight Forces
-            for i, node in enumerate(self.df_shaft["_n"]):
-                SchForce[node + 1] = self.df_shaft.loc[i, "m"] * -9.8065
+        # Shaft Weight Forces
+        for i, node in enumerate(self.df_shaft["_n"]):
+            SCH[node + 1] = self.df_shaft.loc[i, "m"] * -g
+
+        # Organizing data for each shaft
+        BrgForce = []
+        DskForce = []
+        SchForce = []
+        nodes = []
+        nodes_pos = []
+        displacement = []
+        dsk = []
+        brg = []
+        for i in sh_num:
+            n_min = min(
+                aux_rotor.df_shaft.loc[aux_rotor.df_shaft.shaft_number == i, "n_l"]
+            )
+            n_max = max(
+                aux_rotor.df_shaft.loc[(aux_rotor.df_shaft.shaft_number == i), "n_r"]
+            )
+            BrgForce.append(BRG[n_min : n_max + 1])
+            DskForce.append(DSK[n_min : n_max + 1])
+            SchForce.append(SCH[n_min : n_max + 1])
+            nodes_pos.append(self.nodes_pos[n_min : n_max + 1])
+            displacement.append(disp_y[n_min : n_max + 1])
+            nodes.append(list(range(n_min, n_max + 1)))
+
+            # get bearings and disks for each shaft
+            dsk.append(
+                aux_rotor.df_disks.loc[
+                    aux_rotor.df_disks.shaft_number == i, "tag"
+                ].values
+            )
+            brg.append(
+                aux_rotor.df_bearings.loc[
+                    (aux_rotor.df_bearings.shaft_number == i)
+                    | (
+                        aux_rotor.df_bearings.n_link.isin(list(range(n_min, n_max + 1)))
+                    ),
+                    "tag",
+                ].values
+            )
 
+        Mx = []
+        Vx = []
+        Bm = []
+        Vx_axis = []
+        for j in sh_num:
             # Shearing Force vector
-            Vx = [0] * (len(self.nodes_pos))
-            Vx_axis = []
-            for i in range(int(len(self.nodes))):
-                Vx_axis.append(self.nodes_pos[i])
-                Vx[i] = Vx[i - 1] + BrgForce[i] + DskForce[i] + SchForce[i]
-
-            for i in range(len(Vx) + len(self.df_disks) + len(self.df_bearings)):
-                if DskForce[i] != 0:
-                    Vx.insert(i, Vx[i - 1] + SchForce[i])
-                    DskForce.insert(i + 1, 0)
-                    SchForce.insert(i + 1, 0)
-                    BrgForce.insert(i + 1, 0)
-                    Vx_axis.insert(i, Vx_axis[i])
-
-                if BrgForce[i] != 0:
-                    Vx.insert(i, Vx[i - 1] + SchForce[i])
-                    BrgForce.insert(i + 1, 0)
-                    DskForce.insert(i + 1, 0)
-                    SchForce.insert(i + 1, 0)
-                    Vx_axis.insert(i, Vx_axis[i])
-            self.Vx = [x * -1 for x in Vx]
-            Vx = np.array(self.Vx)
+            aux_Vx = [0] * (len(nodes_pos[j]))
+            aux_Vx_axis = [0] * (len(nodes_pos[j]))
+
+            for i in range(int(len(nodes_pos[j]))):
+                aux_Vx_axis[i] = nodes_pos[j][i]
+                aux_Vx[i] = (
+                    aux_Vx[i - 1] + BrgForce[j][i] + DskForce[j][i] + SchForce[j][i]
+                )
+
+            for i in range(len(aux_Vx_axis) + len(dsk[j]) + len(brg[j]) - count):
+                if DskForce[j][i] != 0:
+                    aux_Vx.insert(i, aux_Vx[i - 1] + SchForce[j][i])
+                    DskForce[j].insert(i + 1, 0)
+                    SchForce[j].insert(i + 1, 0)
+                    BrgForce[j].insert(i + 1, 0)
+                    aux_Vx_axis.insert(i, aux_Vx_axis[i])
+
+                if BrgForce[j][i] != 0:
+                    aux_Vx.insert(i, aux_Vx[i - 1] + SchForce[j][i])
+                    BrgForce[j].insert(i + 1, 0)
+                    DskForce[j].insert(i + 1, 0)
+                    SchForce[j].insert(i + 1, 0)
+                    aux_Vx_axis.insert(i, aux_Vx_axis[i])
+
+            aux_Vx = [x * -1 for x in aux_Vx]
+            Vx.append(np.array(aux_Vx))
+            Vx_axis.append(np.array(aux_Vx_axis))
 
             # Bending Moment vector
-            Mx = []
-            for i in range(len(Vx) - 1):
-                if Vx_axis[i] == Vx_axis[i + 1]:
+            aux_Mx = []
+            for i in range(len(aux_Vx) - 1):
+                if aux_Vx_axis[i] == aux_Vx_axis[i + 1]:
                     pass
                 else:
-                    Mx.append(
+                    aux_Mx.append(
                         (
-                            (Vx_axis[i + 1] * Vx[i + 1])
-                            + (Vx_axis[i + 1] * Vx[i])
-                            - (Vx_axis[i] * Vx[i + 1])
-                            - (Vx_axis[i] * Vx[i])
+                            (aux_Vx_axis[i + 1] * aux_Vx[i + 1])
+                            + (aux_Vx_axis[i + 1] * aux_Vx[i])
+                            - (aux_Vx_axis[i] * aux_Vx[i + 1])
+                            - (aux_Vx_axis[i] * aux_Vx[i])
                         )
                         / 2
                     )
+            Mx.append(aux_Mx)
 
-            Bm = np.zeros(1)
-            for i in range(len(Mx)):
-                Bm = np.append(Bm, Bm[i] + Mx[i])
-            self.Bm = Bm
-
-            sh_weight = np.around(sum(self.df_shaft["m"].values) * 9.8065, decimals=1)
-
-            self.shaft_weight = sh_weight
-            self.disk_forces = DskForceToReturn
-            self.bearing_reaction_forces = BrgForceToReturn
-
-            results = StaticResults(
-                self.disp_y,
-                self.Vx,
-                self.Bm,
-                self.df_shaft,
-                self.df_disks,
-                aux_rotor.df_bearings,
-                self.nodes,
-                self.nodes_pos,
-                Vx_axis,
-            )
+            aux_Bm = np.zeros(1)
+            for i in range(len(Mx[j])):
+                aux_Bm = np.append(aux_Bm, aux_Bm[i] + aux_Mx[i])
+            Bm.append(aux_Bm)
+
+        self.Vx = Vx
+        self.Bm = Bm
+        self.disp_y = displacement
+
+        self.w_shaft = [
+            sum(self.df_shaft.loc[self.df_shaft.shaft_number == i, "m"]) * g
+            for i in sh_num
+        ]
+
+        DskForce_nodal = {k: v for k, v in DskForce_nodal.items() if v != 0}
+        BrgForce_nodal = {k: v for k, v in BrgForce_nodal.items() if v != 0}
+        BrgForce_tag = {k: v for k, v in BrgForce_tag.items() if v != 0}
+        DskForce_tag = {k: v for k, v in DskForce_tag.items() if v != 0}
+
+        self.disk_forces_nodal = DskForce_nodal
+        self.bearing_forces_nodal = BrgForce_nodal
+        self.bearing_forces_tag = BrgForce_tag
+        self.disk_forces_tag = DskForce_tag
+
+        results = StaticResults(
+            self.disp_y,
+            self.Vx,
+            self.Bm,
+            self.w_shaft,
+            self.disk_forces_nodal,
+            self.bearing_forces_nodal,
+            nodes,
+            nodes_pos,
+            Vx_axis,
+        )
 
         return results
 
@@ -2338,7 +2467,7 @@ class Rotor(object):
         >>> # show(table)
         """
         self.run_static()
-        forces = self.bearing_reaction_forces
+        forces = self.bearing_forces_tag
         results = SummaryResults(
             self.df_shaft,
             self.df_disks,


Test patch:
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index 7fb2b343..37ad5a4a 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -723,8 +723,8 @@ def test_static_analysis_rotor3(rotor3):
     rotor3.run_static()
 
     assert_almost_equal(
-        rotor3.disp_y,
-        (
+        rotor3.disp_y[0],
+        np.array(
             [
                 -4.94274533e-12,
                 -4.51249085e-04,
@@ -738,25 +738,29 @@ def test_static_analysis_rotor3(rotor3):
         decimal=6,
     )
     assert_almost_equal(
-        rotor3.Vx,
-        [
-            0,
-            -494.2745,
-            -456.6791,
-            -419.0836,
-            -99.4925,
-            -61.8971,
-            -24.3016,
-            480.9807,
-            518.5762,
-            556.1716,
-            0,
-        ],
+        rotor3.Vx[0],
+        np.array(
+            [
+                0,
+                -494.2745,
+                -456.6791,
+                -419.0836,
+                -99.4925,
+                -61.8971,
+                -24.3016,
+                480.9807,
+                518.5762,
+                556.1716,
+                0,
+            ]
+        ),
         decimal=3,
     )
     assert_almost_equal(
-        rotor3.Bm,
-        [0, -118.8692, -228.3395, -248.5132, -259.2881, -134.3434, 0],
+        rotor3.Bm[0],
+        np.array(
+            [0, -118.8692, -228.3395, -248.5132, -259.2881, -134.3434, 0],
+        ),
         decimal=3,
     )
 
@@ -791,8 +795,8 @@ def test_static_analysis_rotor5(rotor5):
     rotor5.run_static()
 
     assert_almost_equal(
-        rotor5.disp_y,
-        (
+        rotor5.disp_y[0],
+        np.array(
             [
                 8.12651626e-04,
                 4.08939282e-04,
@@ -810,41 +814,45 @@ def test_static_analysis_rotor5(rotor5):
         decimal=6,
     )
     assert_almost_equal(
-        rotor5.Vx,
-        [
-            0,
-            37.5954,
-            75.1908,
-            -494.2745,
-            -456.6791,
-            -419.08368,
-            -99.4925,
-            -61.8971,
-            -24.3016,
-            480.9807,
-            518.5762,
-            556.1716,
-            -75.1908,
-            -37.5954,
-            0,
-        ],
+        rotor5.Vx[0],
+        np.array(
+            [
+                0,
+                37.5954,
+                75.1908,
+                -494.2745,
+                -456.6791,
+                -419.08368,
+                -99.4925,
+                -61.8971,
+                -24.3016,
+                480.9807,
+                518.5762,
+                556.1716,
+                -75.1908,
+                -37.5954,
+                0,
+            ]
+        ),
         decimal=3,
     )
     assert_almost_equal(
-        rotor5.Bm,
-        [
-            0,
-            4.6994,
-            18.7977,
-            -100.0714,
-            -209.5418,
-            -229.7155,
-            -240.4903,
-            -115.5457,
-            18.7977,
-            4.6994,
-            0,
-        ],
+        rotor5.Bm[0],
+        np.array(
+            [
+                0,
+                4.6994,
+                18.7977,
+                -100.0714,
+                -209.5418,
+                -229.7155,
+                -240.4903,
+                -115.5457,
+                18.7977,
+                4.6994,
+                0,
+            ]
+        ),
         decimal=3,
     )
 
@@ -880,8 +888,8 @@ def test_static_analysis_rotor6(rotor6):
     rotor6.run_static()
 
     assert_almost_equal(
-        rotor6.disp_y,
-        (
+        rotor6.disp_y[0],
+        np.array(
             [
                 -1.03951876e-04,
                 -4.93624668e-05,
@@ -893,47 +901,51 @@ def test_static_analysis_rotor6(rotor6):
                 1.72933811e-04,
                 -1.02148266e-11,
                 -3.96409257e-04,
-                -9.20006704e-04,
+                -9.20006704e-04
             ]
         ),
         decimal=6,
     )
     assert_almost_equal(
-        rotor6.Vx,
-        [
-            0,
-            37.5954,
-            75.1908,
-            -104.1543,
-            -66.5589,
-            -28.9635,
-            8.6319,
-            328.2230,
-            365.8184,
-            403.4139,
-            441.0093,
-            -580.4733,
-            -542.8778,
-            -505.2824,
-            0,
-        ],
+        rotor6.Vx[0],
+        np.array(
+            [
+                0,
+                37.5954,
+                75.1908,
+                -104.1543,
+                -66.5589,
+                -28.9635,
+                8.6319,
+                328.2230,
+                365.8184,
+                403.4139,
+                441.0093,
+                -580.4733,
+                -542.8778,
+                -505.2824,
+                0,
+            ]
+        ),
         decimal=3,
     )
     assert_almost_equal(
-        rotor6.Bm,
-        [
-            0,
-            4.6994,
-            18.7977,
-            -2.5414,
-            -14.4817,
-            -17.0232,
-            69.7319,
-            165.8860,
-            271.4389,
-            131.0200,
-            0,
-        ],
+        rotor6.Bm[0],
+        np.array(
+            [
+                0,
+                4.6994,
+                18.7977,
+                -2.5414,
+                -14.4817,
+                -17.0232,
+                69.7319,
+                165.8860,
+                271.4389,
+                131.0200,
+                0,
+            ]
+        ),
         decimal=3,
     )