From 3d3c2ebf692647d42225c10242fd0cfbd0899e70 Mon Sep 17 00:00:00 2001 From: stuartmuldrew <58300539+stuartmuldrew@users.noreply.github.com> Date: Fri, 2 Aug 2024 13:44:49 +0100 Subject: [PATCH] 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 <91618944+chris-ashe@users.noreply.github.com> Co-authored-by: mn3981 --- .../proc-pages/physics-models/plasma_beta.md | 7 +++++-- process/physics.py | 12 +++++++++++- source/fortran/input.f90 | 9 ++++++--- source/fortran/physics_variables.f90 | 5 +++++ 4 files changed, 27 insertions(+), 6 deletions(-) diff --git a/documentation/proc-pages/physics-models/plasma_beta.md b/documentation/proc-pages/physics-models/plasma_beta.md index 05f295483..ee64be300 100644 --- 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. This is only recommended for spherical tokamaks .| 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 \ No newline at end of file +[^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 index 47f19602f..ac895a08c 100644 --- 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 index c399e4e4e..9654d2e26 100644 --- 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 index b3ccad7a1..995df178f 100644 --- 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