1047 beta limit (#3264)

* Added beta limit from Tholerus * Added rli calculation to iprofile 6 * Updated documentation * Updated documentation * Updated comment in physics_variables.f90 * Fixed typo * Fix typos in plasma_beta.md * Shortened membership testing for iprofile * Refactor iprofile membership testing for improved efficiency --------- Co-authored-by: Christopher Ashe- UKAEA <[email protected]> Co-authored-by: mn3981 <[email protected]>
ukaea · Aug 2, 2024 · 3d3c2eb · 3d3c2eb
1 parent 66f9488
commit 3d3c2eb
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 6 deletions.
diff --git a/documentation/proc-pages/physics-models/plasma_beta.md b/documentation/proc-pages/physics-models/plasma_beta.md
@@ -28,11 +28,12 @@ be calculated.
 | `iprofile` | Description |
 | :-: | - |
 | 0 | `alphaj`, `rli` and `dnbeta` are inputs. |
-| 1 | `alphaj`, `rli` and `dnbeta` are calulcated consistently. `dnbeta` calculated using $g=4l_i$ [^3].  This is only recommended for high aspect ratio tokamaks.|
+| 1 (default) | `alphaj`, `rli` and `dnbeta` are calulcated consistently. `dnbeta` calculated using $g=4l_i$ [^3].  This is only recommended for high aspect ratio tokamaks.|
 | 2 | `alphaj` and `rli` are inputs. `dnbeta` calculated using $g=2.7(1+5\epsilon^{3.5})$ (which gives g = 3.0 for aspect ratio = 3) |
 | 3 | `alphaj` and `rli` are inputs. `dnbeta` calculated using $g=3.12+3.5\epsilon^{1.7}$ [^4]|
 | 4 | `alphaj` and `dnbeta` are inputs. `rli` calculated from elongation [^4]. This is only recommended for spherical tokamaks.|
 | 5 | `alphaj` is an input.  `rli` calculated from elongation and `dnbeta` calculated using $g=3.12+3.5\epsilon^{1.7}$ [^4]. This is only recommended for spherical tokamaks.|
+| 6 | `alphaj` and `c_beta` are inputs.  `rli` calculated from elongation and `dnbeta` calculated using $C_{\beta}=(g-3.7)F_p / 12.5-3.5F_p$, where $F_p$ is the pressure peaking and $C_{\beta}$ is the destabilisation papermeter (default 0.5)[^5]. See Section 2.4 of Tholerus et al. (2024) for a more detailed description.  <u> This is only recommended for spherical tokamaks <u>.|
 
 Further details on the calculation of `alphaj` and `rli` is given in [Plasma Current](./plasma_current.md).
 
@@ -49,4 +50,6 @@ is be set using input parameter `epbetmax`.
 
 [^3]: Tokamaks 4th Edition, Wesson, page 116
 
-[^4]: Menard et al. (2016), Nuclear Fusion, 56, 106023
+[^4]: Menard et al. (2016), Nuclear Fusion, 56, 106023
+
+[^5]: Tholerus et al. (2024), arXiv:2403.09460
diff --git a/process/physics.py b/process/physics.py
@@ -1400,6 +1400,16 @@ def physics(self):
             # Nucl. Fusion, 2016, 44.
             physics_variables.dnbeta = 3.12e0 + 3.5e0 * physics_variables.eps**1.7e0
 
+        if physics_variables.iprofile == 6:
+            # Method used for STEP plasma scoping
+            # Tholerus et al. (2024), arXiv:2403.09460
+            Fp = (physics_variables.ne0 * physics_variables.te0) / (
+                physics_variables.dene * physics_variables.te
+            )
+            physics_variables.dnbeta = 3.7e0 + (
+                (physics_variables.c_beta / Fp) * (12.5e0 - 3.5e0 * Fp)
+            )
+
         # culblm returns the betalim for beta
         physics_variables.betalim = culblm(
             physics_variables.bt,
@@ -2161,7 +2171,7 @@ def culcur(
             alphaj = qstar / q0 - 1.0
             rli = np.log(1.65 + 0.89 * alphaj)
 
-        if iprofile == 4 or iprofile == 5:
+        if iprofile in [4, 5, 6]:
             # Spherical Tokamak relation for internal inductance
             # Menard et al. (2016), Nuclear Fusion, 56, 106023
             rli = 3.4 - kappa

diff --git a/source/fortran/input.f90 b/source/fortran/input.f90
@@ -304,7 +304,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
       ccl0_ma, ccls_ma, ld_ratio_cst
     use physics_variables, only: ipedestal, taumax, i_single_null, fvsbrnni, &
       rhopedt, cvol, fdeut, ffwal, iculbl, itartpf, ilhthresh, &
-      fpdivlim, epbetmax, isc, kappa95, aspect, cwrmax, nesep, csawth, dene, &
+      fpdivlim, epbetmax, isc, kappa95, aspect, cwrmax, nesep, c_beta, csawth, dene, &
       ftar, plasma_res_factor, ssync, rnbeam, beta, neped, hfact, dnbeta, &
       fgwsep, rhopedn, tratio, q0, ishape, fne0, ignite, ftrit, &
       ifalphap, tauee_in, alphaj, alphat, icurr, q, ti, tesep, rli, triang, &
@@ -549,7 +549,10 @@ subroutine parse_input_file(in_file,out_file,show_changes)
        case ('coreradiationfraction')
           call parse_real_variable('coreradiationfraction', coreradiationfraction, 0.0D0, 1.0D0, &
                'Fraction of core radiation subtracted from P_L')
-       case ('csawth')
+            case ('c_beta')
+               call parse_real_variable('c_beta', c_beta, 0.0D0, 1.0D0, &
+                    'Destabalisation parameter for iprofile=6')
+            case ('csawth')
           call parse_real_variable('csawth', csawth, 0.0D0, 10.0D0, &
                'Coefficient for sawteeth effects')
        case ('cvol')
@@ -647,7 +650,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
           call parse_int_variable('ipedestal', ipedestal, 0, 1, &
                'Switch for plasma profile type')
        case ('iprofile')
-          call parse_int_variable('iprofile', iprofile, 0, 5, &
+          call parse_int_variable('iprofile', iprofile, 0, 6, &
                'Switch for current profile consistency')
        case ('ips')
           call parse_int_variable('ips', ips, 0, 1, &

diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90
@@ -89,6 +89,9 @@ module physics_variables
   real(dp) :: bvert
   !! vertical field at plasma (T)
 
+  real(dp) :: c_beta
+  !! Destabalisation parameter for iprofile=6 beta limit
+
   real(dp) :: csawth
   !! coeff. for sawteeth effects on burn V-s requirement
 
@@ -367,6 +370,7 @@ module physics_variables
   !! - =3 use input values for alphaj, rli. Scale dnbeta with aspect ratio (Menard scaling)
   !! - =4 use input values for alphaj, dnbeta. Set rli from elongation (Menard scaling)
   !! - =5 use input value for alphaj.  Set rli and dnbeta from Menard scaling
+  !! - =6 use input values for alphaj, c_beta.  Set rli from Menard and dnbeta from Tholerus
 
   integer :: iradloss
   !! switch for radiation loss term usage in power balance (see User Guide):
@@ -899,6 +903,7 @@ subroutine init_physics_variables
     burnup = 0.0D0
     burnup_in = 0.0D0
     bvert = 0.0D0
+    c_beta = 0.5D0
     csawth = 1.0D0
     cvol = 1.0D0
     cwrmax = 1.35D0