Magnet solenoid tools (#2304)

* add solenoid tools file

* first pass max field

* first pass hoop stress

* flake8

* max flux

* rename

* docstrings

* axial stress start

* replace estimate with calculation

* remove temp scaffolding

* flake8

* add test

* flkae8

* add livingstone test

bluemira/magnets/solenoid_tools.py

@@ -0,0 +1,153 @@
+# bluemira is an integrated inter-disciplinary design tool for future fusion
+# reactors. It incorporates several modules, some of which rely on other
+# codes, to carry out a range of typical conceptual fusion reactor design
+# activities.
+#
+# Copyright (C) 2021-2023 M. Coleman, J. Cook, F. Franza, I.A. Maione, S. McIntosh,
+#                         J. Morris, D. Short
+#
+# bluemira is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# bluemira is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with bluemira; if not, see <https://www.gnu.org/licenses/>.
+
+"""
+Tools for simple solenoid calculations.
+"""
+
+from typing import Tuple
+
+import numpy as np
+
+from bluemira.magnetostatics.semianalytic_2d import semianalytic_Bx, semianalytic_Bz
+
+
+def calculate_B_max(
+    rho_j: float, r_inner: float, r_outer: float, height: float, z_0: float = 0.0
+) -> float:
+    """
+    Calculate the maximum self-field in a solenoid. This is always located
+    at (r_inner, z_0)
+
+    Parameters
+    ----------
+    rho_j:
+        Current density across the solenoid winding pack [A/m^2]
+    r_inner:
+        Solenoid inner radius [m]
+    r_outer:
+        Solenoid outer radius [m]
+    height:
+        Solenoid vertical extent [m]
+
+    Returns
+    -------
+    Maximum field in a solenoid [T]
+
+    Notes
+    -----
+    Cross-checked graphically with k data from Boom and Livingstone, "Superconducting
+    solenoids", 1962, Fig. 6
+    """
+    dxc = 0.5 * (r_outer - r_inner)
+    xc = r_inner + dxc
+    dzc = 0.5 * height
+    x_bmax = r_inner
+    current = rho_j * (height * (r_outer - r_inner))
+    Bx_max = current * semianalytic_Bx(xc, z_0, x_bmax, z_0, dxc, dzc)
+    Bz_max = current * semianalytic_Bz(xc, z_0, x_bmax, z_0, dxc, dzc)
+    return np.hypot(Bx_max, Bz_max)
+
+
+def calculate_hoop_radial_stress(
+    B_in: float,
+    B_out: float,
+    rho_j: float,
+    r_inner: float,
+    r_outer: float,
+    r: float,
+    poisson_ratio: float = 0.3,
+) -> Tuple[float]:
+    """
+    Calculate the hoop and radial stress at a radial location in a solenoid
+
+    Parameters
+    ----------
+    B_in:
+        Field at the inside edge of the solenoid [T]
+    B_out:
+        Field at the outside edge of the solenoid [T]
+    rho_j:
+        Current density across the solenoid winding pack [A/m^2]
+    r_inner:
+        Solenoid inner radius [m]
+    r_outer:
+        Solenoid outer radius [m]
+    r:
+        Radius at which to calculate [m]
+    poisson_ratio:
+        Poisson ratio of the material
+
+    Returns
+    -------
+    hoop_stress:
+        Hoop stress at the radial location [Pa]
+    radial_stress:
+        Radial stress at the radial location [Pa]
+
+    Notes
+    -----
+    Wilson, Superconducting Magnets, 1982, equations 4.10 and 4.11
+    Must still factor in the fraction of load-bearing material
+    """
+    alpha = r_outer / r_inner
+    eps = r / r_inner
+    nu = poisson_ratio
+    alpha2 = alpha**2
+    eps2 = eps**2
+    ratio2 = alpha2 / eps2
+
+    K = (alpha * B_in - B_out) * rho_j * r_inner / (alpha - 1)  # noqa: N806
+    M = (B_in - B_out) * rho_j * r_inner / (alpha - 1)  # noqa: N806
+    a = K * (2 + nu) / (3 * (alpha + 1))
+    c = M * (3 + nu) / 8
+
+    b = alpha2 + alpha + 1
+    b1 = ratio2 - eps * (1 + 2 * nu) * (alpha + 1) / (2 + nu)
+    b2 = -ratio2 - eps * (alpha + 1)
+
+    d = alpha2 + 1 + ratio2 - eps2 * (1 + 3 * nu) / (3 + nu)
+    e = alpha2 + 1 - ratio2 - eps2
+
+    hoop_stress = a * (b + b1) - c * d
+    radial_stress = a * (b + b2) - c * e
+
+    return hoop_stress, radial_stress
+
+
+def calculate_flux_max(B_max: float, r_inner: float, r_outer: float) -> float:
+    """
+    Calculate the maximum flux achievable from a solenoid
+
+    Parameters
+    ----------
+    B_max:
+        Maximum field in the solenoid [T]
+    r_inner:
+        Solenoid inner radius [m]
+    r_outer:
+        Solenoid outer radius [m]
+
+    Returns
+    -------
+    Maximum flux achievable from a solenoid [V.s]
+    """
+    return np.pi / 3 * B_max * (r_outer**2 + r_inner**2 + r_outer * r_inner)

tests/magnets/test_solenoid_tools.py

@@ -0,0 +1,77 @@
+# bluemira is an integrated inter-disciplinary design tool for future fusion
+# reactors. It incorporates several modules, some of which rely on other
+# codes, to carry out a range of typical conceptual fusion reactor design
+# activities.
+#
+# Copyright (C) 2021-2023 M. Coleman, J. Cook, F. Franza, I.A. Maione, S. McIntosh,
+#                         J. Morris, D. Short
+#
+# bluemira is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# bluemira is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with bluemira; if not, see <https://www.gnu.org/licenses/>.
+
+
+import numpy as np
+import pytest
+
+from bluemira.base.constants import MU_0
+from bluemira.magnetostatics.semianalytic_2d import semianalytic_Bx, semianalytic_Bz
+from bluemira.magnets.solenoid_tools import calculate_B_max, calculate_hoop_radial_stress
+
+
+class TestHoopRadialStress:
+    def test_radial_boundary_conditions(self):
+        x, z = 4, 0
+        dx, dz = 0.25, 0.8
+        rho_j = 1e6
+        nu = 0.3
+        current = rho_j * (4 * dx * dz)
+        Bx_in = current * semianalytic_Bx(x - dx, z, x, z, dx, dz)
+        Bz_in = current * semianalytic_Bz(x - dx, z, x, z, dx, dz)
+        B_in = np.hypot(Bx_in, Bz_in)
+        Bx_out = current * semianalytic_Bx(x + dx, z, x, z, dx, dz)
+        Bz_out = current * semianalytic_Bz(x + dx, z, x, z, dx, dz)
+        B_out = np.hypot(Bx_out, Bz_out)
+        sigma_theta_in, sigma_r_in = calculate_hoop_radial_stress(
+            B_in, B_out, rho_j, x - dx, x + dx, x - dx, nu
+        )
+        np.testing.assert_almost_equal(sigma_r_in, 0.0)
+        sigma_theta_out, sigma_r_out = calculate_hoop_radial_stress(
+            B_in, B_out, rho_j, x - dx, x + dx, x + dx, nu
+        )
+        np.testing.assert_almost_equal(sigma_r_out, 0.0)
+        assert sigma_theta_in > sigma_theta_out
+
+
+class TestCalculateBmax:
+    @pytest.mark.parametrize(
+        "alpha,beta,k_expected",
+        [[1.375, 1.2, 1.12], [1.322, 1.5, 1.07], [1.287, 2.0, 1.03]],
+    )
+    def test_Bmax(self, alpha, beta, k_expected):
+        """
+        Expected values from Boom and Livingstone, 1962, Example 4 table (unlabelled)
+        """
+        r_inner = 3
+        r_outer = alpha * r_inner
+        dz = beta * r_inner
+        rho_j = 1e6
+        B_max = calculate_B_max(rho_j, r_inner, r_outer, 2 * dz, 0)
+        B_0 = (
+            rho_j
+            * r_inner
+            * MU_0
+            * beta
+            * np.log((alpha + np.hypot(alpha, beta)) / (1 + np.hypot(1, beta)))
+        )
+        k_calc = B_max / B_0
+        np.testing.assert_almost_equal(k_calc, k_expected, decimal=2)