🚧 3388 update the beta limit section of the docs to be up to date and include all models

documentation/proc-pages/development/standards.md

@@ -306,15 +306,15 @@ subroutine constraint_eqn_001(args)
   !! - \( T_i \) -- density weighted average ion temperature [keV]
   !! - \( B_{tot} \) -- total toroidal + poloidal field [T]
 
-  use physics_variables, only: betaft, beta_beam, dene, ten, dnitot, tin, btot, beta
+  use physics_variables, only: beta_fast_alpha, beta_beam, dene, ten, dnitot, tin, btot, beta
   use constants, only: electron_charge,rmu0
 
   implicit none
 
   type(constraint_args_type), intent(out) :: args
   !! constraint derived type
 
-    args%cc = 1.0D0 - (betaft + beta_beam + &
+    args%cc = 1.0D0 - (beta_fast_alpha + beta_beam + &
       2.0D3*rmu0*electron_charge * (dene*ten + dnitot*tin)/btot**2 )/beta
     args%con = beta * (1.0D0 - args%cc)
     args%err = beta * args%cc

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

@@ -0,0 +1,195 @@
+# Plasma Beta
+
+## Overview
+
+The efficiency of confinement of plasma pressure by the magnetic field is represented by the ratio:
+
+$$
+\beta = \frac{2\mu_0p}{B^2}
+$$
+
+There are several different measures of this type, arising from different choices of definition and from the need to quantify different equilibrium properties.
+
+In its expanded form of magnetic field components:
+
+$$
+\beta = \frac{2\mu_0p}{B^2_{\text{T}}+B^2_{\text{P}}}
+$$
+
+This is often broken down into separate related quantities knows as the toroidal $\beta_{\text{t}}$ and the poloidal beta $\beta_{\text{p}}$ whose definitions are:
+
+$$
+\beta_{\text{t}} = \frac{2\mu_0p}{B^2_{\text{T}}} \\
+\beta_{\text{p}} = \frac{2\mu_0p}{B^2_{\text{P}}}
+$$
+
+The relationship between these quantities can be written as:
+
+$$
+\frac{1}{\beta} = \frac{1}{\beta_{\text{t}}} + \frac{1}{\beta_{\text{p}}}
+$$
+
+implying that the total $\beta$ is dominated by the smaller of the two contributions.Observe that by definition $\beta \le 1$. However,either $\beta_{\text{t}}$ or $\beta_{\text{p}}$, but not both, can be greater than unity.
+
+The above in its simplest form is known as the total thermal beta $\beta_{\text{th}}$ as it only takes into account the pressure produced by the ions and electrons such that:
+
+$$
+\beta_{\text{th}} = \frac{2\mu_0 \langle n_{\text{e}}T_{\text{e}}+n_{\text{i}}T_{\text{i}}\rangle}{B^2}
+$$
+
+In future reactors there will be a notable pressure contribution due to the productions of fast alpha particles and that from plasma-beam interactions if used. So the true total beta can be defined as:
+
+$$
+\beta_{\text{tot}} = \frac{2\mu_0 \langle n_{\text{e}}T_{\text{e}}+n_{\text{i}}T_{\text{i}} + p_{\alpha}+ p_{\text{beam}}\rangle}{B^2} \\
+\beta_{\text{tot}} = \beta_{\text{th}}+\beta_{\alpha}+ \beta_{\text{beam}}
+$$
+
+Models for the fast alpha particle pressure contribution can be found [here](plasma_alpha_beta_contribution.md).
+
+----------------
+
+## Derivation of plasma beta parameter
+
+From the lowest order momentum balance equation:
+
+$$
+\nabla p = \mathbf{j} \times \mathbf{B}
+$$
+
+as it pertains to plasma equilibrium. The current and the field must also satisfy Maxwell’s equations
+
+$$
+\mu_0\mathbf{j} = \nabla \times \mathbf{B},  \quad \nabla \cdot \mathbf{B} = 0
+$$
+
+A consequence of the above is that the current and  magnetic field lie on isobaric surfaces. This follows directly from the observation that $\mathbf{j} \cdot \nabla p = \mathbf{B} \cdot \nabla p = 0$ and $\nabla p$ is everywhere normal to the surface $p = \text{const}$. The force is everywhere normal to the isobaric surface and just balances the pressure gradient force, $-\nabla p$. Although the current and the magnetic field lie in a common flux surface, they can only be parallel in regions where the pressure gradient vanishes. Since currents that are parallel to the field do not contribute to the $\mathbf{j} \times \mathbf{B}$ force, they are known as “force-free currents.”
+
+Taking the divergence of $\mu_0\mathbf{j} = \nabla \times \mathbf{B}$  yields $\nabla \cdot \mathbf{j} = 0$ which is consistent with the quasineutrality assumption. The vector product of $\mathbf{B}$ with the momentum balance equation leads to an expression for the current perpendicular to $\mathbf{B}$,
+
+$$
+\mathbf{j}_{\perp} = \frac{\mathbf{B} \times \nabla p}{B^2}
+$$
+
+We see that there can be no perpendicular current in the absence of a pressure gradient. This leads to a momentum balance expressed in terms of the pressure and magnetic field.
+
+$$
+\nabla \left(p+\frac{B^2}{2\mu_0}\right) = \frac{1}{\mu_0}\left(\mathbf{B} \cdot \nabla \right)\mathbf{B}
+$$
+
+The right hand side vanishes when the field lines are straight and parallel, in which case above reduces to a simple statement that the total (kinetic plus magnetic) pressure is constant everywhere within a confined plasma, 
+
+$$
+p + \frac{B^2}{2\mu_0} \equiv \frac{B^2}{2\mu_0}\left(1+\beta \right) = \text{const}
+$$
+
+where we have introduced the quantity
+
+$$
+\beta \equiv \frac{2\mu_0p}{B^2}
+$$
+
+------------------------
+
+## Beta Limit
+
+The plasma beta limit[^1] is given by 
+
+$$\begin{aligned}
+\beta < 0.01\, g \, \frac{I(\mbox{MA})}{a(\mbox{m}) \, B_0(\mbox{T})}
+\end{aligned}$$
+
+where $B_0$ is the axial vacuum toroidal field. The beta
+coefficient $g$ is set using input parameter `dnbeta`. To apply the beta limit, 
+constraint equation 24 should be turned on with iteration variable 36
+(`fbetatry`). 
+
+By default, $\beta$ is defined with respect to the total equilibrium B-field [^2]. 
+
+| `i_beta_component` | Description |
+| :-: | - |
+| 0 (default) | Apply the $\beta$ limit to the total plasma beta (including the contribution from fast ions) |
+| 1 | Apply the $\beta$ limit to only the thermal component of beta |
+| 2 | Apply the $\beta$ limit to only the thermal plus neutral beam contributions to beta |
+| 3 | Apply the $\beta$ limit to the total beta (including the contribution from fast ions), calculated using only the toroidal field |
+
+### Setting the Beta $g$ Coefficient
+
+Switch `iprofile` determines how the beta $g$ coefficient `dnbeta` should 
+be calculated.
+
+| `iprofile` | Description |
+| :-: | - |
+| 0 | `alphaj`, `rli` and `dnbeta` are inputs. |
+| 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).
+
+----------------------
+
+## Key Constraints
+
+### Beta consistency
+
+This constraint can be activated by stating `icc = 1` in the input file.
+
+Ensures the relationship between $\beta$, density, temperature and total magnetic field is withheld by checking the fixed input or iteration variable $\mathtt{beta}$ is consistent in value with the rest of the physics parameters
+
+$$
+\mathtt{beta} \equiv \frac{2\mu_0 \langle n_{\text{e}}T_{\text{e}}+n_{\text{i}}T_{\text{i}}\rangle}{B^2} + \beta_{\alpha} + \beta_{\text{beam}}
+$$
+
+**It is highly recommended to always have this constraint on as it is a global consistency checker**
+
+----------------
+
+### Poloidal beta and inverse aspect upper limit
+
+This constraint can be activated by stating `icc = 6` in the input file [^6].
+
+To apply a limit to the value of $\epsilon\beta_p$, where $\epsilon = a/R$ is
+the inverse aspect ratio and $\beta_p$ is the poloidal $\beta$, constraint equation no. 6 should be 
+turned on with iteration variable no. 8 (`fbeta`). The limiting value of $\epsilon\beta_p$ 
+is be set using input parameter `epbetmax`.
+
+--------------------
+
+### Beta upper limit
+
+This constraint can be activated by stating `icc = 24` in the input file.
+
+It is the general setting of the $\beta$ limit depending on the $\beta_{\text{N}}$ value calculated in the [beta limit](#beta-limit) calculations.
+
+The upper limit value of beta is calculated by `calculate_beta_limit()`. The scaling value `fbetatry` can be varied also.
+
+**It is recommended to have this constraint on as it is a plasma stability model**
+
+--------------------
+
+### Poloidal upper limit
+
+This constraint can be activated by stating `icc = 48` in the input file.
+
+The value of `beta_poloidal_max` can be set to the desired maximum poloidal beta. The scaling value `fbeta_poloidal` can be varied also.
+
+-------------------
+
+### Beta lower limit
+
+This constraint can be activated by stating `icc = 84` in the input file.
+
+[^1]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989',
+
+[^2]: D.J. Ward, 'PROCESS Fast Alpha Pressure', Work File Note F/PL/PJK/PROCESS/CODE/050
+
+[^3]: Tokamaks 4th Edition, Wesson, page 116
+
+[^4]: Menard et al. (2016), Nuclear Fusion, 56, 106023
+
+[^5]: Tholerus et al. (2024), arXiv:2403.09460
+
+[^6]: M. E. Mauel et al., “Operation at the tokamak equilibrium poloidal beta-limit in TFTR,” Nuclear Fusion, vol. 32, no. 8, pp. 1468–1473, Aug. 1992. doi:10.1088/0029-5515/32/8/i14 

examples/data/csv_output_large_tokamak_MFILE.DAT

@@ -343,15 +343,15 @@
  Safety_factor_at_95%_flux_surface_______________________________________ (q95)_________________________      3.4251E+00 ITV
  Cylindrical_safety_factor_(qcyl)________________________________________ (qstar)_______________________      2.8387E+00    
  Total_plasma_beta_______________________________________________________ (beta)________________________      3.6055E-02 ITV
- Total_poloidal_beta_____________________________________________________ (betap)_______________________      1.2531E+00 OP 
+ Total_poloidal_beta_____________________________________________________ (beta_poloidal)_______________________      1.2531E+00 OP 
  Total_toroidal_beta_____________________________________________________ ______________________________      3.7123E-02 OP 
  Fast_alpha_beta_________________________________________________________ (betaft)______________________      4.4814E-03 OP 
  Beam_ion_beta___________________________________________________________ (betanb)______________________      0.0000E+00 OP 
  (Fast_alpha_+_beam_beta)/(thermal_beta)_________________________________ (gammaft)_____________________      1.4194E-01 OP 
  Thermal_beta____________________________________________________________ ______________________________      3.1574E-02 OP 
  Thermal_poloidal_beta___________________________________________________ ______________________________      1.0974E+00 OP 
  Thermal_toroidal_beta_(=_beta-exp)______________________________________ ______________________________      3.2509E-02 OP 
- 2nd_stability_beta_:_beta_p_/_(R/a)_____________________________________ (eps*betap)___________________      4.1770E-01    
+ 2nd_stability_beta_:_beta_p_/_(R/a)_____________________________________ (eps*beta_poloidal)___________________      4.1770E-01    
  2nd_stability_beta_upper_limit__________________________________________ (epbetmax)____________________      1.3800E+00    
  Beta_g_coefficient______________________________________________________ (dnbeta)______________________      4.7593E+00    
  Normalised_thermal_beta_________________________________________________ ______________________________      2.5527E+00    
@@ -1571,7 +1571,7 @@ gamma = 0.3
 i_bootstrap_current = 4
 
 * Switch for beta limit scaling
-iculbl = 1
+i_beta_component = 1
 
 * Switch for plasma current scaling
 i_plasma_current    = 4

examples/data/large_tokamak_1_MFILE.DAT

@@ -344,15 +344,15 @@
  Safety_factor_at_95%_flux_surface_______________________________________ (q95)_________________________      3.5961E+00 ITV
  Cylindrical_safety_factor_(qcyl)________________________________________ (qstar)_______________________      2.9805E+00    
  Total_plasma_beta_______________________________________________________ (beta)________________________      3.3648E-02 ITV
- Total_poloidal_beta_____________________________________________________ (betap)_______________________      1.2857E+00 OP 
+ Total_poloidal_beta_____________________________________________________ (beta_poloidal)_______________________      1.2857E+00 OP 
  Total_toroidal_beta_____________________________________________________ ______________________________      3.4553E-02 OP 
  Fast_alpha_beta_________________________________________________________ (betaft)______________________      4.2236E-03 OP 
  Beam_ion_beta___________________________________________________________ (betanb)______________________      0.0000E+00 OP 
  (Fast_alpha_+_beam_physics_variables.beta)/(thermal_physics_variables.be (gammaft)_____________________      1.4354E-01 OP 
  Thermal_beta____________________________________________________________ ______________________________      2.9425E-02 OP 
  Thermal_poloidal_beta___________________________________________________ ______________________________      1.1243E+00 OP 
  Thermal_toroidal_physics_variables.beta_(=_beta-exp)____________________ ______________________________      3.0215E-02 OP 
- 2nd_stability_physics_variables.beta_:_beta_p_/_(R/a)___________________ (eps*betap)___________________      4.2857E-01    
+ 2nd_stability_physics_variables.beta_:_beta_p_/_(R/a)___________________ (eps*beta_poloidal)___________________      4.2857E-01    
  2nd_stability_physics_variables.beta_upper_limit________________________ (epbetmax)____________________      1.3800E+00    
  Beta_g_coefficient______________________________________________________ (dnbeta)______________________      4.9099E+00    
  Normalised_thermal_beta_________________________________________________ ______________________________      2.4977E+00    
@@ -1565,7 +1565,7 @@ gamma = 0.3
 i_bootstrap_current = 4
 
 * Switch for beta limit scaling
-iculbl = 1
+i_beta_component = 1
 
 * Switch for plasma current scaling
 i_plasma_current    = 4