954 bootstrap current fraction scaling formula for steady state machines

documentation/proc-pages/physics-models/plasma_current.md

@@ -68,6 +68,7 @@ existence of pedestals, whereas the Sauter et al. scaling
 | 2 | General scaling -- To use a more general scaling method, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
 | 3 | Numerically fitted scaling [^8] -- To use a numerically fitted scaling method, valid for all aspect ratios, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
 | 4 | Sauter, Angioni and Lin-Liu scaling [^9] [^10] -- Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
+| 5 | Sakai, Fujita and Okamoto scaling [^11] -- Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). The model includes the toroidal diamagnetic current in the calculation due to the dataset, so `idia = 0` can only be used with it
 
 !!! Note "Fixed Bootstrap Current"
     Direct input -- To input the bootstrap current fraction directly, set `bscfmax` 
@@ -114,4 +115,5 @@ Task meeting EU-JA #16, Fusion for Energy, Garching, 24--25 June 2014
 [^7]: Menard et al. (2016), Nuclear Fusion, 56, 106023
 [^8]: H.R. Wilson, Nuclear Fusion **32** (1992) 257
 [^9]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **6** (1999) 2834 
-[^10]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **9** (2002) 5140
+[^10]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **9** (2002) 5140
+[^11]: Ryosuke Sakai, Takaaki Fujita, Atsushi Okamoto, Derivation of bootstrap current fraction scaling formula for 0-D system code analysis, Fusion Engineering and Design, Volume 149, 2019, 111322, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2019.111322.

process/physics.py

@@ -990,6 +990,19 @@ def physics(self):
             current_drive_variables.cboot * self.bootstrap_fraction_sauter()
         )
 
+        current_drive_variables.bscf_sakai = (
+            current_drive_variables.cboot
+            * self.bootstrap_fraction_sakai(
+                betap=physics_variables.betap,
+                q95=physics_variables.q95,
+                q0=physics_variables.q0,
+                alphan=physics_variables.alphan,
+                alphat=physics_variables.alphat,
+                eps=physics_variables.eps,
+                rli=physics_variables.rli,
+            )
+        )
+
         if current_drive_variables.bscfmax < 0.0e0:
             current_drive_variables.bootipf = abs(current_drive_variables.bscfmax)
             current_drive_variables.plasipf = current_drive_variables.bootipf
@@ -1002,6 +1015,10 @@ def physics(self):
                 current_drive_variables.bootipf = current_drive_variables.bscf_wilson
             elif physics_variables.ibss == 4:
                 current_drive_variables.bootipf = current_drive_variables.bscf_sauter
+            elif physics_variables.ibss == 5:
+                # Sakai states that the ACCOME dataset used has the toridal diamagnetic current included in the bootstrap current
+                # So the diamagnetic current calculation should be turned off when using, (idia = 0).
+                current_drive_variables.bootipf = current_drive_variables.bscf_sakai
             else:
                 error_handling.idiags[0] = physics_variables.ibss
                 error_handling.report_error(75)
@@ -4064,6 +4081,14 @@ def outplas(self):
                 current_drive_variables.bscf_wilson,
                 "OP ",
             )
+
+            po.ovarrf(
+                self.outfile,
+                "Bootstrap fraction (Sakai)",
+                "(bscf_sakai)",
+                current_drive_variables.bscf_sakai,
+                "OP ",
+            )
             po.ovarrf(
                 self.outfile,
                 "Diamagnetic fraction (Hender)",
@@ -4115,6 +4140,11 @@ def outplas(self):
                     self.outfile,
                     "  (Sauter et al bootstrap current fraction model used)",
                 )
+            elif physics_variables.ibss == 5:
+                po.ocmmnt(
+                    self.outfile,
+                    "  (Sakai et al bootstrap current fraction model used)",
+                )
 
             if physics_variables.idia == 0:
                 po.ocmmnt(
@@ -4662,6 +4692,55 @@ def bootstrap_fraction_sauter():
 
         return np.sum(da * jboot, axis=0) / physics_variables.plascur
 
+    @staticmethod
+    def bootstrap_fraction_sakai(
+        betap: float,
+        q95: float,
+        q0: float,
+        alphan: float,
+        alphat: float,
+        eps: float,
+        rli: float,
+    ) -> float:
+        """
+        Calculate the bootstrap fraction using the Sakai formula.
+
+        Parameters:
+        betap (float): Plasma poloidal beta.
+        q95 (float): Safety factor at 95% of the plasma radius.
+        q0 (float): Safety factor at the magnetic axis.
+        alphan (float): Density profile index
+        alphat (float): Temperature profile index
+        eps (float): Inverse aspect ratio.
+
+        Returns:
+        float: The calculated bootstrap fraction.
+
+        Notes:
+        The profile assumed for the alphan anf alpat indexes is only a prabolic profile without a pedestal (L-mode).
+        The Root Mean Squared Error for the fitting database of this formula was 0.025
+        Concentrating on the positive shear plasmas using the ACCOME code equilibria with the fully non-inductively driven
+        conditions with neutral beam (NB) injection only are calculated.
+        The electron temperature and the ion temperature were assumed to be equal
+        This can be used for all apsect ratios.
+        The diamagnetic fraction is included in this formula.
+
+        References:
+        Ryosuke Sakai, Takaaki Fujita, Atsushi Okamoto, Derivation of bootstrap current fraction scaling formula for 0-D system code analysis,
+        Fusion Engineering and Design, Volume 149, 2019, 111322, ISSN 0920-3796,
+        https://doi.org/10.1016/j.fusengdes.2019.111322.
+        """
+        # Sakai states that the ACCOME dataset used has the toridal diamagnetic current included in the bootstrap current
+        # So the diamganetic current should not be calculated with this. idia = 0
+        return (
+            10 ** (0.951 * eps - 0.948)
+            * betap ** (1.226 * eps + 1.584)
+            * rli ** (-0.184 - 0.282)
+            * (q95 / q0) ** (-0.042 * eps - 0.02)
+            * alphan ** (0.13 * eps + 0.05)
+            * alphat ** (0.502 * eps - 0.273)
+        )
+
     def fhfac(self, is_):
         """Function to find H-factor for power balance
         author: P J Knight, CCFE, Culham Science Centre

process/utilities/errorlist.json

@@ -8,7 +8,7 @@
   "comment2": [
     "Increment n_errortypes if an error is added to this list"
   ],
-  "n_errortypes": 283,
+  "n_errortypes": 284,
   "errors": [
     {
       "no": 1,
@@ -1424,6 +1424,11 @@
       "no": 283,
       "level": 3,
       "message": "[tfcoil]: Can only have i_tf_turns_integer = 1 with i_tf_wp_geom = 0"
+    },
+    {
+      "no": 284,
+      "level": 3,
+      "message": "CHECK: idia = 0 should be used with the Sakai plasma current scaling."
     }
   ]
 }

source/fortran/current_drive_variables.f90

@@ -42,6 +42,9 @@ module current_drive_variables
   real(dp) :: bscf_wilson
   !! bootstrap current fraction, Wilson et al model
 
+  real(dp) :: bscf_sakai
+  !! Bootstrap current fraction, Sakai et al model
+
   real(dp) :: cboot
   !! bootstrap current fraction multiplier (`ibss=1`)
 

source/fortran/initial.f90

@@ -258,7 +258,7 @@ subroutine check
     use physics_variables, only: aspect, fdeut, fgwped, fhe3, &
         fgwsep, ftrit, ibss, i_single_null, icurr, idivrt, ishape, &
         iradloss, isc, ipedestal, ilhthresh, itart, nesep, rhopedn, rhopedt, &
-        rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar
+        rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar, idia
     use pulse_variables, only: lpulse
     use reinke_variables, only: fzactual, impvardiv
     use tfcoil_variables, only: casthi, casthi_is_fraction, casths, i_tf_sup, &
@@ -828,6 +828,10 @@ subroutine check
         call report_error(283)
     end if
 
+    if ( ibss == 5  .and. idia /= 0 ) then
+        call report_error(284)
+    end if
+
     ! Setting t_cable_tf_is_input to true if t_cable_tf is an input
     if ( abs(t_cable_tf) < epsilon(t_cable_tf) ) then
         t_cable_tf_is_input = .false.

source/fortran/input.f90

@@ -619,7 +619,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
           call parse_real_variable('taumax', taumax, 0.1D0, 100.0D0, &
                'Maximum allowed energy confinement time (s)')
        case ('ibss')
-          call parse_int_variable('ibss', ibss, 1, 4, &
+          call parse_int_variable('ibss', ibss, 1, 5, &
                'Switch for bootstrap scaling')
        case ('iculbl')
           call parse_int_variable('iculbl', iculbl, 0, 3, &

source/fortran/physics_variables.f90

@@ -249,6 +249,7 @@ module physics_variables
   !! - =2 for Nevins et al general scaling
   !! - =3 for Wilson et al numerical scaling
   !! - =4 for Sauter et al scaling
+  !! - =5 for Sakai et al scaling
 
   integer :: iculbl
   !! switch for beta limit scaling (`constraint equation 24`)

tests/unit/test_physics.py

@@ -468,6 +468,90 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph
     assert bfs == pytest.approx(bootstrapfractionsauterparam.expected_bfs)
 
 
+class BootstrapFractionSakaiParam(NamedTuple):
+    betap: Any = None
+
+    q95: Any = None
+
+    q0: Any = None
+
+    alphan: Any = None
+
+    alphat: Any = None
+
+    eps: Any = None
+
+    rli: Any = None
+
+    expected_bfs: Any = None
+
+
+@pytest.mark.parametrize(
+    "bootstrapfractionsakaiparam",
+    (
+        BootstrapFractionSakaiParam(
+            betap=1.3184383457774960,
+            q95=3.5151046634673557,
+            q0=1.0,
+            alphan=1.0,
+            alphat=1.45,
+            eps=1 / 3,
+            rli=1.2098126022585098,
+            expected_bfs=0.34204201506155418,
+        ),
+        BootstrapFractionSakaiParam(
+            betap=1.1701245502231756,
+            q95=5.1746754543339177,
+            q0=2.0,
+            alphan=0.9,
+            alphat=1.3999999999999999,
+            eps=1 / 1.8,
+            rli=0.3,
+            expected_bfs=0.90349498124262029,
+        ),
+    ),
+)
+def test_bootstrap_fraction_sakai(bootstrapfractionsakaiparam, monkeypatch, physics):
+    """
+    Automatically generated Regression Unit Test for bootstrap_fraction_sakai.
+
+    This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT.
+    and tests/regression/input_files/st_regression.IN.DAT.
+
+    :param bootstrapfractionsauterparam: the data used to mock and assert in this test.
+    :type bootstrapfractionsauterparam: bootstrapfractionsauterparam
+
+    :param monkeypatch: pytest fixture used to mock module/class variables
+    :type monkeypatch: _pytest.monkeypatch.monkeypatch
+    """
+
+    monkeypatch.setattr(physics_variables, "betap", bootstrapfractionsakaiparam.betap)
+
+    monkeypatch.setattr(physics_variables, "q95", bootstrapfractionsakaiparam.q95)
+
+    monkeypatch.setattr(physics_variables, "q0", bootstrapfractionsakaiparam.q0)
+
+    monkeypatch.setattr(physics_variables, "alphan", bootstrapfractionsakaiparam.alphan)
+
+    monkeypatch.setattr(physics_variables, "alphat", bootstrapfractionsakaiparam.alphat)
+
+    monkeypatch.setattr(physics_variables, "eps", bootstrapfractionsakaiparam.eps)
+
+    monkeypatch.setattr(physics_variables, "rli", bootstrapfractionsakaiparam.rli)
+
+    bfs = physics.bootstrap_fraction_sakai(
+        betap=bootstrapfractionsakaiparam.betap,
+        q95=bootstrapfractionsakaiparam.q95,
+        q0=bootstrapfractionsakaiparam.q0,
+        alphan=bootstrapfractionsakaiparam.alphan,
+        alphat=bootstrapfractionsakaiparam.alphat,
+        eps=bootstrapfractionsakaiparam.eps,
+        rli=bootstrapfractionsakaiparam.rli,
+    )
+
+    assert bfs == pytest.approx(bootstrapfractionsakaiparam.expected_bfs)
+
+
 class CulcurParam(NamedTuple):
     normalised_total_beta: Any = None