Merge pull request #1280 from kiracofe8/rotor_Id_summary

Rotor id summary

Author: raphaeltimbo

docs/user_guide/tutorial_part_5.ipynb

@@ -4592,10 +4592,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "A = [[0, 1],\n",
-    "     [-2, -0.5]]\n",
-    "B = [[0],\n",
-    "     [1]]\n",
+    "A = [[0, 1], [-2, -0.5]]\n",
+    "B = [[0], [1]]\n",
     "C = [[1, 0]]\n",
     "\n",
     "Q_lqr = np.diag([10, 1])\n",
@@ -4631,7 +4629,7 @@
    "outputs": [],
    "source": [
     "C_lead = lead_lag(tau=0.02, alpha=0.1, k=2.0)\n",
-    "C_lag  = lead_lag(tau=0.5, alpha=5.0, k=1.0)"
+    "C_lag = lead_lag(tau=0.5, alpha=5.0, k=1.0)"
    ]
   },
   {
@@ -4723,7 +4721,7 @@
     "    pid(2, 1, 0.02, n_f=200),\n",
     "    lead_lag(0.03, 0.1),\n",
     "    notch_filter(120, 0.02, 0.2),\n",
-    "    low_pass_filter(400)\n",
+    "    low_pass_filter(400),\n",
     ")"
    ]
   },
@@ -4784,12 +4782,14 @@
    ],
    "source": [
     "plot_frequency_response(\n",
-    "    C, N, C * N,\n",
+    "    C,\n",
+    "    N,\n",
+    "    C * N,\n",
     "    legends=[\"Controller\", \"Filter\", \"Combined\"],\n",
     "    w_min=1e-1,\n",
     "    w_max=1e5,\n",
     "    n_points=1500,\n",
-    "    title=\"Bode (Controller / Filter / Combined)\"\n",
+    "    title=\"Bode (Controller / Filter / Combined)\",\n",
     ")"
    ]
   },

ross/results.py

@@ -5290,13 +5290,14 @@ class SummaryResults(Results):
         The figure object with the tables plot.
     """
 
-    def __init__(self, df_shaft, df_disks, df_bearings, brg_forces, CG, Ip, tag):
+    def __init__(self, df_shaft, df_disks, df_bearings, brg_forces, CG, Ip, It, tag):
         self.df_shaft = df_shaft
         self.df_disks = df_disks
         self.df_bearings = df_bearings
         self.brg_forces = brg_forces
         self.CG = CG
         self.Ip = Ip
+        self.It = It
         self.tag = tag
 
     def plot(self):
@@ -5334,6 +5335,7 @@ def plot(self):
             "Total lenght": [self.df_shaft["nodes_pos_r"].iloc[-1]],
             "CG": ["{:.3f}".format(self.CG)],
             "Ip": ["{:.3e}".format(self.Ip)],
+            "It": ["{:.3e}".format(self.It)],
             "Rotor Mass": [
                 "{:.3f}".format(np.sum(self.df_shaft["m"]) + np.sum(self.df_disks["m"]))
             ],

ross/rotor_assembly.py

@@ -567,6 +567,24 @@ def flatten(l):
         self.G0 = G0
         self.Ksdt0 = Ksdt0
 
+        # Calculation of overall rotor transverse (diametral) inertia (includes only DOFs located at the shaft element DOF, excludes point masses that are outside the shaft).
+        # This is only calculating Iyy. Assuming Ixx is the same.
+        # First, set up a vector corresponding to rigid body rotation of the entire rotor
+        v = np.zeros([self.ndof])
+        for i, elm in enumerate(self.shaft_elements):
+            dofs = list(elm.dof_global_index.values())
+            y0 = elm.dof_mapping()["y_0"]
+            a0 = elm.dof_mapping()["alpha_0"]
+            y1 = elm.dof_mapping()["y_1"]
+            a1 = elm.dof_mapping()["alpha_1"]
+
+            v[dofs[y0]] = -(nodes_pos_l[i] - self.CG)  # y
+            v[dofs[y1]] = -(nodes_pos_r[i] - self.CG)  # y
+            v[dofs[a0]] = 1  # alpha
+            v[dofs[a1]] = 1  # alpha
+        # Then, use the vector to compute diametral aka transverse inertia of the entire rotor.
+        self.It = v @ (self.M0 @ v.T)
+
     def _set_tag(self, tag):
         """Set the tag for the current rotor."""
         self.tag = tag or "Rotor 0"
@@ -4334,6 +4352,7 @@ def summary(self):
             forces,
             self.CG,
             self.Ip,
+            self.It,
             self.tag,
         )
         return results

ross/tests/test_rotor_assembly.py

@@ -2907,3 +2907,153 @@ def test_run_amb_sensitivity():
         results.max_abs_sensitivities["Magnetic Bearing 0"]["x"],
         results_custom_freq.max_abs_sensitivities["Magnetic Bearing 0"]["x"],
     )
+
+
+@pytest.fixture
+def rotor3a():
+    #  same as rotor 3 but no disks
+    i_d = 0
+    o_d = 0.05
+    n = 6
+    L = [0.25 for _ in range(n)]
+
+    shaft_elem = [
+        ShaftElement(
+            l,
+            i_d,
+            o_d,
+            material=steel,
+            shear_effects=True,
+            rotary_inertia=True,
+            gyroscopic=True,
+        )
+        for l in L
+    ]
+
+    stfx = 1e6
+    stfy = 0.8e6
+    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=0)
+    bearing1 = BearingElement(6, kxx=stfx, kyy=stfy, cxx=0)
+
+    return Rotor(shaft_elem, [], [bearing0, bearing1])
+
+
+@pytest.fixture
+def rotor3b():
+    #  #  same as rotor 3, disks at center.
+    i_d = 0
+    o_d = 0.05
+    n = 6
+    L = [0.25 for _ in range(n)]
+
+    shaft_elem = [
+        ShaftElement(
+            l,
+            i_d,
+            o_d,
+            material=steel,
+            shear_effects=True,
+            rotary_inertia=True,
+            gyroscopic=True,
+        )
+        for l in L
+    ]
+
+    disk0 = DiskElement.from_geometry(3, steel, 0.07, 0.05, 0.28)
+    disk1 = DiskElement.from_geometry(3, steel, 0.07, 0.05, 0.35)
+
+    stfx = 1e6
+    stfy = 0.8e6
+    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=0)
+    bearing1 = BearingElement(6, kxx=stfx, kyy=stfy, cxx=0)
+
+    return Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+
+
+@pytest.fixture
+def rotor3c():
+    #  no density on shaft elements.  test ONLY the It of the disks
+    i_d = 0
+    o_d = 0.05
+    n = 6
+    L = [0.25 for _ in range(n)]
+
+    mat0 = Material(name="Steel0", rho=0, E=211e9, Poisson=0.3)
+
+    shaft_elem = [
+        ShaftElement(
+            l,
+            i_d,
+            o_d,
+            material=mat0,
+            shear_effects=True,
+            rotary_inertia=True,
+            gyroscopic=True,
+        )
+        for l in L
+    ]
+
+    disk0 = DiskElement.from_geometry(0, steel, 0.07, 0.05, 0.28)
+    disk1 = DiskElement.from_geometry(6, steel, 0.07, 0.05, 0.35)
+
+    stfx = 1e6
+    stfy = 0.8e6
+    bearing0 = BearingElement(0, kxx=stfx, kyy=stfy, cxx=0)
+    bearing1 = BearingElement(6, kxx=stfx, kyy=stfy, cxx=0)
+
+    return Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1])
+
+
+# this test checks if Rotor.m calculates mass the same way as directly summing the mass matrix
+def test_rotor_m(rotor3):
+    mm = rotor3.M()[::6, ::6].sum()
+
+    assert rotor3.m == pytest.approx(mm, rel=1e-8)
+
+
+# this test checks if Rotor.CG calculates the CG same way as a direct vector calcution from the mass matrix
+def test_rotor_CG(rotor3):
+    mm = rotor3.M()[::6, ::6]
+    mm2 = np.sum(mm, axis=0)
+    CG2 = np.matmul(mm2, rotor3.nodes_pos) / rotor3.m
+
+    assert rotor3.CG == pytest.approx(CG2, rel=1e-8)
+
+
+def test_rotor_it(rotor1, rotor3, rotor3a, rotor3b, rotor3c):
+    L = 0.5
+    r = 0.05 / 2
+    It1 = (rotor1.m / 12) * (3 * r**2 + L**2)
+    tol = 2e-8
+
+    assert rotor1.It == pytest.approx(It1, rel=tol)  # check with no disks
+
+    L3 = 0.25 * 6
+    It3a = (rotor3a.m / 12) * (3 * r**2 + L3**2)
+    assert rotor3a.It == pytest.approx(It3a, rel=tol)  # second check with no disks
+
+    # disk0 It=0.17808928257067666, m=32.58972765304033
+    # disk1 It=0.42358058405433824, m=51.5252611115262
+
+    It3b = It3a + 0.17808928257067666 + 0.42358058405433824
+    assert rotor3b.It == pytest.approx(It3b, rel=tol)  # check with disks at CG
+
+    It3c = (
+        0.17808928257067666
+        + 32.58972765304033 * (0 - rotor3c.CG) ** 2
+        + 0.42358058405433824
+        + 51.5252611115262 * (1.5 - rotor3c.CG) ** 2
+    )
+    assert rotor3c.It == pytest.approx(
+        It3c, rel=tol
+    )  # check with ONLY disks off center
+
+    # check with off center disks and shaft elements
+    It3 = (
+        (It3a + (rotor3.CG - 0.75) ** 2 * rotor3a.m)
+        + 0.17808928257067666
+        + 32.58972765304033 * (0.5 - rotor3.CG) ** 2
+        + 0.42358058405433824
+        + 51.5252611115262 * (1.0 - rotor3.CG) ** 2
+    )
+    assert rotor3.It == pytest.approx(It3, rel=tol)