diff --git a/documentation/proc-pages/eng-models/heating-and-current-drive.md b/documentation/proc-pages/eng-models/heating-and-current-drive.md deleted file mode 100644 index cea6e5c553..0000000000 --- a/documentation/proc-pages/eng-models/heating-and-current-drive.md +++ /dev/null @@ -1,69 +0,0 @@ -# Auxiliary Power Systems - -## Current Drive - -The use of inductive current drive leads to pulsed plant operation because of the limited flux swing that can be achieved using the central solenoid. This poses problems due to the fact that fatigue failures may result, and there would be a need for thermal storage to maintain output of electricity between pulses, and supply power for starting a new pulse.However, the plasma current can also be produced and maintained (partially or wholly) using non-inductive means which, in principle, removes this restriction. `PROCESS` contains a number of auxiliary current drive schemes, including various RF methods (Lower Hybrid, Electron Cyclotron,Electron Bernstein Wave, and Ion Cyclotron (Fast Wave) current drives) and also Neutral Beam current drive systems. The code calculates the efficiency and the resulting power requirements of the chosen system. - -The fraction of the required plasma current to be produced by non-inductive means, `fvsbrnni`, should be set, and flag `irfcd` should be set to 0 for purely inductive scenarios, or 1 otherwise. The current drive efficiency model to be used in this latter case is defined by the value of switch `iefrf`: - -- `iefrf` = 1: Fenstermacher Lower Hybrid model -- `iefrf` = 2: Ion cyclotron model[^1], -- `iefrf` = 3: Fernstermacher electron cyclotron resonance model -- `iefrf` = 4: Ehst Lower Hybrid model -- `iefrf` = 5: ITER neutral beam model[^1],[^2], -- `iefrf` = 6: Culham Lower Hybrid model[^2], -- `iefrf` = 7: Culham electron cyclotron model[^2], -- `iefrf` = 8: Culham neutral beam model[^2], -- `iefrf` = 9: Oscillating Field current drive (RFPs only - OBSOLETE-REMOVED), -- `iefrf` = 10: ECRH user input gamma, -- `iefrf` = 11: ECRH "HARE" model [^3], -- `iefrf` = 12: EBW user scaling input[^4], -- `iefrf` = 13: ECRH O-mode cutoff with Zeff and Te [^4], - - -(Note that, at present, the neutral beam models do not include the effect of an edge transport barrier (pedestal) in the plasma profile.) - -It is sometimes useful to adjust artificially the current drive efficiency values produced by these routines. This can be achieved by setting teh scaling coefficients `feffcd`. The wall plug to plasma efficiencies can also be adjusted, by changing the relevant variable (`etaech`, `etalh`, `etanbi` or `etaof`). - -## Plasma heating - -In addition to current drive, some auxiliary power can be used to heat the plasma. The value of input parameters `pheat` determines the amount of auxiliary *heating* power (in Watts) to be applied to the plasma. This variable may be used as an iteration variable (no. 11). - -## Neutral beam access - -If present, a neutral beam injection system needs sufficient space between the TF coils to be able to intercept the plasma tangentially. The major radius `rtanbeam` at which the centre-line of the beam is tangential to the toroidal direction is user-defined using input parameter `frbeam`, which is the ratio of `rtanbeam` to the plasma major radius `rmajor`. The maximum possible tangency radius `rtanmax` is determined by the geometry of the TF coils - see Figure 1, and this can be enforced using constraint equation no. 20 with iteration variable no. 33 (`fportsz`). The thickness of the beam duct walls may be set using input parameter `nbshield`. - -
-
- NBI port -

-
Figure 1: Top-down schematic of the neutral beam access geometry. The beam with the maximum possible tangency radius is shown here. -
-
-
-
- -## Neutral beam losses - -Input parameter `forbitloss` can be used to specify the fraction of the net injected neutral beam power that is lost between the beam particles' ionisation and thermalisation (known as the first orbit loss). This quantity cannot easily be calculated as it depends on the field ripple and other three-dimensional effects. The power lost is assumed to be absorbed by the first wall. - -The power in the beam atoms that are not ionised as they pass through the plasma (shine-through) is calculated by the code. There are two constraint equations that can be used to control the beam penetration and deposition, as follows: - -- It is necessary to use a beam energy that simultaneously gives adequate penetration of the beam to the centre of the plasma and tolerable shine-through of the beam on the wall after the beam has traversed the plasma. The number of exponential decay lengths, $\tau$, for the beam power to fall before it reaches the plasma centre should be in the region of ~ 4-6[^2],. Constraint equation no. 14 may be used to force $\tau$ to be equal to the value given by input parameter `tbeamin`, and is therefore in effect a beam energy consistency equation. -- Alternatively, constraint equation no. 59 with iteration variable no. 105 (`fnbshineef`) may be used to ensure that the beam power fraction emerging from the plasma is no more than the value given by input parameter `nbshinefmax`. - -It is recommended that only one of these two constraint equations is used during a run. - -## Ignited plasma - -Switch `ignite` can be used to denote whether the plasma is ignited, i.e. fully self-sustaining without the need for any injected auxiliary power during the burn. If `ignite` = 1, the calculated injected power does not contribute to the plasma power balance, although the cost of the auxiliary power system is taken into account (the system is then assumed to be required to provide heating etc during the plasma start-up phase only - use `pheat` to indicate the power requirement). If `ignite` = 0, the plasma is not ignited, and the auxiliary power is taken into account in the plasma power balance during the burn phase. Also, constraint equation no. 28 can be turned on to enforce the fusion gain *Q* to be at least `bigqmin`. - -[^1]: N. A. Uckan and ITER Physics Group, *"ITER Physics Design Guidelines: 1989"*, ITER Documentation Series, No. 10, IAEA/ITER/DS/10 (1990) - -[^2]: T. C. Hender, M. K. Bevir, M. Cox, R. J. Hastie, P. J. Knight, C. N. Lashmore-Davies, B. Lloyd, G. P. Maddison, A. W. Morris, M. R. O'Brien, M.F. Turner abd H. R. Wilson, *"Physics Assessment for the European Reactor Study"*, AEA Fusion Report AEA FUS 172 (1992) - -[^3]: E. Poli, M. Müller, H. Zohm, M. Kovari, *"Fast evaluation of the current driven by electron cyclotron waves for reactor studies"*, Physics of Plasmas 1 December 2018; 25 (12): 122501 - -[^4]: Laqua, H & Maassberg, H & Marushchenko, Nikolai & Volpe, Francesco & Weller, A & Kasparek, W. (2003). *"Electron-Bernstein-Wave Current Drive in an Overdense Plasma at the Wendelstein 7-AS Stellarator"*, Physical review letters. 90. 075003. 10.1103/PhysRevLett.90.075003. \ No newline at end of file diff --git a/documentation/proc-pages/fusion-devices/spherical-tokamak.md b/documentation/proc-pages/fusion-devices/spherical-tokamak.md index b628d1c6ee..1a48d8ca6f 100644 --- a/documentation/proc-pages/fusion-devices/spherical-tokamak.md +++ b/documentation/proc-pages/fusion-devices/spherical-tokamak.md @@ -43,8 +43,8 @@ Switch `itart` provides overall control of the ST switches within the code, and | Switch | Conventional aspect ratio (`itart = 0`) | Low aspect ratio (`itart = 1`) | | --- | --- | --- | | `ishape` | 0, 2, 3, 4 | 1, 5, 6, 7, 8 | -| `ibss` | 1, 2, 3 | 2, 3 | -| `icurr` | 1, 3, 4, 5, 6, 7 | 2, 9 | +| `i_bootstrap_current` | 1, 2, 3 | 2, 3 | +| `i_plasma_current` | 1, 3, 4, 5, 6, 7 | 2, 9 | | `itfsup` | 0, 1 | 0 |
Table 1: Summary of the switch values in 'PROCESS' that relate to conventional aspect ratio and low aspect ratio machines.
diff --git a/documentation/proc-pages/physics-models/error.txt b/documentation/proc-pages/physics-models/error.txt index 2f9529e681..f7d77d834e 100644 --- a/documentation/proc-pages/physics-models/error.txt +++ b/documentation/proc-pages/physics-models/error.txt @@ -1042,15 +1042,15 @@ Scalings}\label{bootstrap-diamagnetic-and-pfirsch-schluxfcter-current-scalings} The fraction of the plasma current provided by the so-called bootstrap effect can be either input into the code directly, or calculated using one of four methods, as summarised here. Note that methods -\texttt{ibss\ =\ 1-3} use fits to parabolic density and temperature +\texttt{i_bootstrap_current\ =\ 1-3} use fits to parabolic density and temperature profiles, and do not take into account the existence of pedestals (\texttt{ipedestal\ =\ 1}), whereas the Sauter et al. scaling -(\texttt{ibss\ =\ 4}) allows general profiles to be used. +(\texttt{i_bootstrap_current\ =\ 4}) allows general profiles to be used. \begin{longtable}[]{@{}cl@{}} \toprule \begin{minipage}[b]{0.05\columnwidth}\centering\strut -\texttt{ibss}\strut +\texttt{i_bootstrap_current}\strut \end{minipage} & \begin{minipage}[b]{0.03\columnwidth}\raggedright\strut Description\strut \end{minipage}\tabularnewline @@ -1060,7 +1060,7 @@ Description\strut 1\strut \end{minipage} & \begin{minipage}[t]{0.03\columnwidth}\raggedright\strut ITER scaling -- To use the ITER scaling method for the bootstrap current -fraction. Set \texttt{bscfmax} to the maximum required bootstrap current +fraction. Set \texttt{bootstrap_current_fraction_max} to the maximum required bootstrap current fraction (\(\leq 1\)). This method is valid at high aspect ratio only.\strut \end{minipage}\tabularnewline @@ -1068,7 +1068,7 @@ only.\strut 2\strut \end{minipage} & \begin{minipage}[t]{0.03\columnwidth}\raggedright\strut General scaling\footnotemark{} -- To use a more general scaling method, -set \texttt{bscfmax} to the maximum required bootstrap current fraction +set \texttt{bootstrap_current_fraction_max} to the maximum required bootstrap current fraction (\(\leq 1\)).\strut \end{minipage} \footnotetext{W.M. Nevins, `Summary Report: ITER Specialists' Meeting on @@ -1078,7 +1078,7 @@ set \texttt{bscfmax} to the maximum required bootstrap current fraction 3\strut \end{minipage} & \begin{minipage}[t]{0.03\columnwidth}\raggedright\strut Numerically fitted scaling\footnotemark{} -- To use a numerically fitted -scaling method, valid for all aspect ratios, set \texttt{bscfmax} to the +scaling method, valid for all aspect ratios, set \texttt{bootstrap_current_fraction_max} to the maximum required bootstrap current fraction (\(\leq 1\)).\strut \end{minipage} \footnotetext{H.R. Wilson, Nuclear Fusion \textbf{32} (1992) 257}\tabularnewline @@ -1086,7 +1086,7 @@ maximum required bootstrap current fraction (\(\leq 1\)).\strut 4\strut \end{minipage} & \begin{minipage}[t]{0.03\columnwidth}\raggedright\strut Sauter, Angioni and Lin-Liu scaling\footnotemark{} \footnotemark{} -- -Set \texttt{bscfmax} to the maximum required bootstrap current fraction +Set \texttt{bootstrap_current_fraction_max} to the maximum required bootstrap current fraction (\(\leq 1\)).\strut \end{minipage} \addtocounter{footnote}{-1} @@ -1099,7 +1099,7 @@ Set \texttt{bscfmax} to the maximum required bootstrap current fraction \end{longtable} !!! Note ``Fixed Bootstrap Current'' Direct input -- To input the -bootstrap current fraction directly, set \texttt{bscfmax} to \((-1)\) +bootstrap current fraction directly, set \texttt{bootstrap_current_fraction_max} to \((-1)\) times the required value (e.g. -0.73 sets the bootstrap faction to 0.73). @@ -2547,7 +2547,7 @@ Overfull \hbox (3.56793pt too wide) in paragraph at lines 1038--1043 \T1/lmr/m/n/10 (-20) tion. Overfull \hbox (22.65782pt too wide) in paragraph at lines 1038--1043 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1038--1043 \T1/lmr/m/n/10 (-20) max- @@ -2619,7 +2619,7 @@ Overfull \hbox (19.63986pt too wide) in paragraph at lines 1046--1050 \T1/lmr/m/n/10 (-20) method, Overfull \hbox (22.65782pt too wide) in paragraph at lines 1046--1050 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1046--1050 \T1/lmr/m/n/10 (-20) max- @@ -2679,7 +2679,7 @@ Overfull \hbox (2.53915pt too wide) in paragraph at lines 1056--1060 \T1/lmr/m/n/10 (-20) tios, Overfull \hbox (22.65782pt too wide) in paragraph at lines 1056--1060 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1056--1060 \T1/lmr/m/n/10 (-20) max- @@ -2727,7 +2727,7 @@ Overfull \hbox (7.44649pt too wide) in paragraph at lines 1064--1068 \T1/lmr/m/n/10 (-20) ing[][][] Overfull \hbox (22.65782pt too wide) in paragraph at lines 1064--1068 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1064--1068 \T1/lmr/m/n/10 (-20) max- @@ -4677,7 +4677,7 @@ Overfull \hbox (3.56793pt too wide) in paragraph at lines 1038--1043 \T1/lmr/m/n/10 (-20) tion. Overfull \hbox (22.65782pt too wide) in paragraph at lines 1038--1043 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1038--1043 \T1/lmr/m/n/10 (-20) max- @@ -4749,7 +4749,7 @@ Overfull \hbox (19.63986pt too wide) in paragraph at lines 1046--1050 \T1/lmr/m/n/10 (-20) method, Overfull \hbox (22.65782pt too wide) in paragraph at lines 1046--1050 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1046--1050 \T1/lmr/m/n/10 (-20) max- @@ -4809,7 +4809,7 @@ Overfull \hbox (2.53915pt too wide) in paragraph at lines 1056--1060 \T1/lmr/m/n/10 (-20) tios, Overfull \hbox (22.65782pt too wide) in paragraph at lines 1056--1060 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1056--1060 \T1/lmr/m/n/10 (-20) max- @@ -4857,7 +4857,7 @@ Overfull \hbox (7.44649pt too wide) in paragraph at lines 1064--1068 \T1/lmr/m/n/10 (-20) ing[][][] Overfull \hbox (22.65782pt too wide) in paragraph at lines 1064--1068 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1064--1068 \T1/lmr/m/n/10 (-20) max- @@ -6789,7 +6789,7 @@ Overfull \hbox (3.56793pt too wide) in paragraph at lines 1038--1043 \T1/lmr/m/n/10 (-20) tion. Overfull \hbox (22.65782pt too wide) in paragraph at lines 1038--1043 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1038--1043 \T1/lmr/m/n/10 (-20) max- @@ -6861,7 +6861,7 @@ Overfull \hbox (19.63986pt too wide) in paragraph at lines 1046--1050 \T1/lmr/m/n/10 (-20) method, Overfull \hbox (22.65782pt too wide) in paragraph at lines 1046--1050 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1046--1050 \T1/lmr/m/n/10 (-20) max- @@ -6921,7 +6921,7 @@ Overfull \hbox (2.53915pt too wide) in paragraph at lines 1056--1060 \T1/lmr/m/n/10 (-20) tios, Overfull \hbox (22.65782pt too wide) in paragraph at lines 1056--1060 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1056--1060 \T1/lmr/m/n/10 (-20) max- @@ -6969,7 +6969,7 @@ Overfull \hbox (7.44649pt too wide) in paragraph at lines 1064--1068 \T1/lmr/m/n/10 (-20) ing[][][] Overfull \hbox (22.65782pt too wide) in paragraph at lines 1064--1068 -\T1/lmtt/m/n/10 bscfmax +\T1/lmtt/m/n/10 bootstrap_current_fraction_max Overfull \hbox (5.74335pt too wide) in paragraph at lines 1064--1068 \T1/lmr/m/n/10 (-20) max- diff --git a/documentation/proc-pages/physics-models/plasma_current.md b/documentation/proc-pages/physics-models/plasma_current.md deleted file mode 100644 index 1a2a3f3029..0000000000 --- a/documentation/proc-pages/physics-models/plasma_current.md +++ /dev/null @@ -1,119 +0,0 @@ -# Plasma Current Scaling Laws - -A number of plasma current scaling laws are available in PROCESS [^1]. These are calculated in -routine `culcur`, which is called by `physics`. The safety factor $q_{95}$ required to prevent -disruptive MHD instabilities dictates the plasma current Ip: - -$$\begin{aligned} -I_p = \frac{2\pi}{\mu_0} B_t \frac{a^2 f_q}{Rq_{95}} -\end{aligned}$$ - -The factor $f_q$ makes allowance for toroidal effects and plasma shaping (elongation and -triangularity). Several formulae for this factor are available [2,3] depending on the value of -the switch `icurr`, as follows: - -| `icurr` | Description | -| :-: | - | -| 1 | Peng analytic fit | -| 2 | Peng double null divertor scaling (ST)[^4] | -| 3 | Simple ITER scaling | -| 4 | Revised ITER scaling[^5] $f_q = \frac{1.17-0.65\epsilon}{2(1-\epsilon^2)^2} (1 + \kappa_{95}^2 (1+2\delta_{95}^2 - 1.2\delta_{95}^3) )$| -| 5 | Todd empirical scaling, I | -| 6 | Todd empirical scaling, II | -| 7 | Connor-Hastie model | -| 8 | Sauter model, allows negative $\delta$ | -| 9 | Scaling for spherical tokamaks, based on a fit to a family of equilibria derived by Fiesta: $f_q = 0.538 (1 + 2.44\epsilon^{2.736}) \kappa^{2.154} \delta^{0.06}$| - -## Plasma Current Profile Consistency - -A limited degree of self-consistency between the plasma current profile and other parameters [^6] can be -enforced by setting switch `iprofile = 1`. This sets the current -profile peaking factor $\alpha_J$ (`alphaj`), the normalised internal inductance $l_i$ (`rli`) and beta limit $g$-factor (`dnbeta`) using the -safety factor on axis `q0` and the cylindrical safety factor $q*$ (`qstar`): - -$$\begin{aligned} -\alpha_J = \frac{q*}{q_0} - 1 -\end{aligned}$$ - -$$\begin{aligned} -l_i = \rm{ln}(1.65+0.89\alpha_J) -\end{aligned}$$ - -$$\begin{aligned} -g = 4 l_i -\end{aligned}$$ - -It is recommended that current scaling law `icurr = 4` is used if `iprofile = 1`. -This relation is only applicable to large aspect ratio tokamaks. - -For spherical tokamaks, the normalised internal inductance can be set from the elongation using `iprofile = 4` or `iprofile = 5`: - -$$\begin{aligned} -l_i = 3.4 - \kappa_x -\end{aligned}$$ - -Further desciption of `iprofile` is given in [Beta Limit](./plasma_beta.md). - -## Bootstrap, Diamagnetic and Pfirsch-Schlüter Current Scalings - -The fraction of the plasma current provided by the bootstrap effect -can be either input into the code directly, or calculated using one of four -methods, as summarised here. Note that methods `ibss = 1-3` do not take into account the -existence of pedestals, whereas the Sauter et al. scaling -(`ibss = 4`) allows general profiles to be used. - -| `ibss` | Description | -| :-: | - | -| 1 | ITER scaling -- To use the ITER scaling method for the bootstrap current fraction. Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). This method is valid at high aspect ratio only. -| 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` - to $(-1)$ times the required value (e.g. -0.73 sets the bootstrap faction to 0.73). - -The diamagnetic current fraction $f_{dia}$ is strongly related to $\beta$ and is typically small, -hence it is usually neglected. For high $\beta$ plasmas, such as those at tight -aspect ratio, it should be included and two scalings are offered. If the diamagnetic -current is expected to be above one per cent of the plasma current, a warning -is issued to calculate it. - -`idia = 0` Diamagnetic current fraction is zero. - -`idia = 1` Diamagnetic current fraction is calculated using a fit to spherical tokamak calculations by Tim Hender: - -$$f_{dia} = \frac{\beta}{2.8}$$ - -`idia = 2` Diamagnetic current fraction is calculated using a SCENE fit for all aspect ratios: - -$$f_{dia} = 0.414 \space \beta \space (\frac{0.1 q_{95}}{q_0} + 0.44)$$ - -A similar scaling is available for the Pfirsch-Schlüter current fraction $f_{PS}$. This is -typically smaller than the diamagnetic current, but is negative. - -`ips = 0` Pfirsch-Schlüter current fraction is set to zero. - -`ips = 1` Pfirsch-Schlüter current fraction is calculated using a SCENE fit for all aspect ratios: - -$$ f_{PS} = -0.09 \beta $$ - -There is no ability to input the diamagnetic and Pfirsch-Schlüter current -directly. In this case, it is recommended to turn off these two scalings -and to use the method of fixing the bootstrap current fraction. - -[^1]: D.J. Ward, 'PROCESS Fast Alpha Pressure', Work File Note F/PL/PJK/PROCESS/CODE/050 -[^2]: Albajar, Nuclear Fusion **41** (2001) 665 -[^3]: M. Kovari, R. Kemp, H. Lux, P. Knight, J. Morris, D.J. Ward, '“PROCESS”: A systems code for fusion power plants—Part 1: Physics' Fusion Engineering and Design 89 (2014) 3054–3069 -[^4]: J.D. Galambos, 'STAR Code : Spherical Tokamak Analysis and Reactor Code', -Unpublished internal Oak Ridge document. -[^5]: W.M. Nevins, 'Summary Report: ITER Specialists' Meeting on Heating and -Current Drive', ITER-TN-PH-8-4, 13--17 June 1988, Garching, FRG -[^6]: Y. Sakamoto, 'Recent progress in vertical stability analysis in JA', -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 -[^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. diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md new file mode 100644 index 0000000000..e3c4472200 --- /dev/null +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -0,0 +1,504 @@ +## Overview + +Bootstrap current in tokamaks originates from the pressure gradients within the plasma and the resulting collisions between particles. As the plasma pressure varies radially, it creates a differential in particle velocities, leading to a net drift of charged particles. This drift generates a toroidal current, known as the bootstrap current, which flows parallel to the magnetic field lines. The phenomenon is a consequence of the neoclassical transport theory, where the collisional processes in a magnetically confined plasma lead to a self-sustaining current. This current is particularly advantageous as it reduces the need for external current drive systems, thereby enhancing the efficiency and stability of the tokamak operation. The bootstrap current is proportional to the pressure gradient and the collisionality of the plasma, making it a critical factor in the design and operation of advanced tokamak reactors aiming for steady-state fusion. + +Some more info can be found [here](https://wiki.fusion.ciemat.es/wiki/Bootstrap_current) +## Selection + +The fraction of the plasma current provided by the bootstrap effect +can be either input into the code directly, or calculated using one of five +methods, as summarised here. Note that methods `i_bootstrap_current = 1-3 & 5` do not take into account the +existence of pedestals, whereas the Sauter et al. scaling +(`i_bootstrap_current = 4`) allows general profiles to be used. + +-------------- + +### ITER IPDG89 Scaling | `bootstrap_fraction_iter89()` + +Original empirical ITER bootstrap scaling from the 1989 Physics Design Guidelines[^0] +Is selected by setting `i_bootstrap_current = 1` + +Empirical fit for the bootstrap current fraction as: + +$$ +\frac{I_{\text{BS}}}{I} = C_{\text{BS}}\left(\epsilon^{0.5}\beta_{\text{pa}}\right)^{1.3} +$$ + +where: + +$$ +C_{\text{BS}} = 1.32 - 0.235\left(\frac{q_{95}}{q_0}\right) + 0.0185\left(\frac{q_{95}}{q_0}\right)^2 +$$ + +$$ +\beta_{\text{pa}} = \frac{\langle p \rangle}{\frac{B_{\text{pa}^2}}{2\mu_0}} = \beta_{\text{tot}}\left(\frac{B_0}{B_{\text{pa}}}\right)^2 +$$ + +$$ +B_{\text{pa}} = \frac{I}{5\langle a \rangle} +$$ + +$$ +\langle a \rangle = \left(\frac{V}{2\pi^2 R_0}\right)^{0.5} +$$ + +Here, $\beta_{\text{tot}}$ is the average total plasma (toroidal) beta. $I$ is given in $\text{MA}$ and $B_0$ is the on-axis toroidal field in Tesla. + +------------ + +### Nevins Scaling | `bootstrap_fraction_nevins()` + +The general Nevins scaling is normally cited from a ITER specialists meeting in 1989[^1] which is not publicly accessible. +However it can be found in the appendix [here](https://doi.org/10.1016/j.fusengdes.2014.07.009)[^2]. + +Is selected by setting `i_bootstrap_current = 2` + +$$ +f_{\text{BS}} = \frac{2.5\beta_{e0}R_pB_{T}q_{95}}{I_{\text{P}}}\int^1_0 B_{\text{int}}\ dy +$$ + +$$ +\beta_e = \frac{1.6022 \times 10^{-16}n_{\text{e}}T_{\text{e}}}{B_{T}^2 / 2\mu_0} +$$ + +$$ +\beta_{e0} = \frac{1.6022 \times 10^{-16}n_{\text{e0}}T_{\text{e0}}}{B_{T}^2 / 2\mu_0} +$$ + +$$ +\Delta = \epsilon y^{0.5} +$$ + +$$ +x = \frac{1.46 \sqrt{\Delta}+2.4\Delta}{\left(1-\Delta\right)^{1.5}} +$$ + +$$ +d = \sqrt{2}Z_{\text{eff}} + Z_{\text{eff}}^2 + x\left(0.75+2.657Z_{\text{eff}}+2Z_{\text{eff}}^2\right) \\ ++x^2\left(0.348+1.243Z_{\text{eff}}+Z_{\text{eff}}^2\right) +$$ + +$$ +A_1 = \left(\alpha_{\text{n}}+\alpha_{\text{T}}\right)\left(1-y\right)^{\alpha_{\text{n}}+\alpha_{\text{T}}-1} +$$ + +$$ +A_2 = \alpha_{\text{T}}\left(1-y\right)^{\alpha_{\text{n}}+\alpha_{\text{T}}-1} +$$ + +$$ +A_{l1} = \frac{x}{d}\left(0.754+2.21Z_{\text{eff}}+Z_{\text{eff}}^2+x\left(0.348+1.243Z_{\text{eff}}+Z_{\text{eff}}^2\right)\right) +$$ + +$$ +A_{l2}=-x\frac{0.884+2.0742Z_{\text{eff}}}{d} +$$ + +$$ +\alpha_i = -\frac{1.172}{1+0.462x} +$$ + +$$ +q = q_0+\left(q_{95}-q_0\right)\frac{y+y^2+y^3}{3} +$$ + +$$ +\beta_{\text{tot}} = \beta_{\text{T}}\frac{B_{\text{tot}}^2}{B_{\text{T}}^2} = \frac{\beta_{T}\left(B_{T}^2+\beta_{\text{P}}\right)}{B_{\text{T}}^2} +$$ + +$$ +P_{\text{ratio}} = \frac{\beta_{\text{T}}-\beta_e}{\beta_e} +$$ + +$$ +B_{\text{int}} = \frac{q}{q_{95}}\left[A_{l1}\left(A_1+P_{\text{ratio}}\left(A_1+\alpha_i A_2\right)\right)+A_{l2}A_2\right] +$$ + + + +---------------- + + +### Wilson Scaling | `bootstrap_fraction_wilson()` + +Wilson gives an empirical formula[^3] [^7] as a function of the pressure, +temperature and total current profiles as well as the poloidal beta and aspect ratio of a tokamak. This empirical formula was compared with an expression obtained by the ITER group; also compared with an analytical result (valid at large aspect ratio). It is found that the determined empirical result agreed well with the large aspect ratio result, but not so well with the empirical formula of the ITER group[^0] + +Is selected by setting `i_bootstrap_current = 3` + +Data is fitted to 3000 equilibria evenly distributed across the parameter range below: + +| Variable | $R$ | $A$ | $B_{\text{T}}$ | $\delta$ | $\kappa$ | $\alpha_{\text{P}}$ | $\alpha_{\text{T}}$ | $\alpha_{\text{J}}$ | $Z$ | +|----------|-----|-----|----------------|-----------|-----------|-------------------|-------------------|-------------------|-----| +| Value | 5.31 | 1.1-5 | 6.2 | 0.2 | 2.0 | 1-3 | 0.1-$\alpha_{\text{P}}$ | 0.5-2.0 | 1-3 | + + +Using the relation: + +$$ +\frac{I_{\text{b}}}{I_{\text{P}}} = \beta_{\text{P}}\epsilon_0^{\frac{1}{2}} \sum_{i=1}^{12}a_i(\alpha_{\text{J}}, Z)b_i +$$ + +In the paper definition $\epsilon_0$ is not the standard inverse aspect ration but is defined as: + +$$ +\epsilon_0 = \frac{R_2-R_1}{R_2+R_1} +$$ + +Where $R_2$ and $R_1$ are the maximum and minimum radii of the plasma. + +The poloidal beta term is defined as: + +$$ +\beta_{\text{P}} = \frac{2\mu_0\langle p \rangle}{\langle \langle B_{\text{P}} \rangle \rangle^2} +$$ + +Where $\langle p \rangle$ is the volume averaged pressure and $\langle \langle B_{\text{P}} \rangle \rangle$ is the plasma surface average of the polidal field. + +The Wilson method extrapolates on the analytically derived expression for the bootstrap current fraction in the limit of large aspect ratio with circular cross-section. This large aspect expression assumed coincidence of a constant $r$ surface with a constant flux surface. The allowed the treatment of the temperature, pressure and current density as functions of $r$ only with the simple parabolic profile form. +For the Wilson method this is expanded to the arbitrary aspect ratio case with D-shaped plasmas with a given triangularity and elongation. The profiles for pressure and temperature are now of the form: + +$$ +P=P_0\psi^{\alpha_{\text{P}}} +$$ + +Where $\psi$ is the flux function. For the current density we cannot model it in the form above as it varies across a flux surface. Instead it is convenient to consider a flux surface average of the total parallel current density which is taken to vary as: + +$$ +\frac{\langle j \cdot B \rangle}{\langle B^2 \rangle^{\frac{1}{2}}} = J_0 \psi^{\alpha_{\text{J}}} +$$ + +Describing this flux surface average by using a standard **parabolic profile** for the safety factor $(q)$. + +$$ +q(r) = q_0 +\left(q_{95}-q_0\right)\left(\frac{r}{a}\right)^2 +$$ + +To relate the temperature, pressure and current density profiles across a flux surface average of the flux function we re-arrange their standard parabolic forms into one that can be integrated into the $(q)$ profile function. Using the temperature profile as an example: + +$$ +T(r) = T_0 \left(1-\left(\frac{r}{a}\right)^2\right)^{\alpha_{\text{T}}} +$$ + +$$ +\frac{T(r)}{T_0} = \left(1-\left(\frac{r}{a}\right)^2\right)^{\alpha_{\text{T}}} +$$ + +$$ +\left(\frac{r}{a}\right)^2 = 1-\left(\frac{T(r)}{T_0}\right)^{\frac{1}{\alpha_{\text{T}}}} +$$ + +Substituting the above into the $q$ profile and taking a median profile value of half the core value we get: + +$$ +q_0+(q_{95}-q_0)\times (1-0.5^{\frac{1}{\alpha_{\text{T}}}}) +$$ + +To find the new flux surface average where the $q$ profile shares points where the temperature, pressure and current density are half of their core values, we relate the function above of the calculated $q$ value to the core $q_0$ value. Taking the natural log of the $q$ profile point at 50% of the core value for the temperature, pressure and current density we relate these two ratios compared to $q_{95}$. Dividing $\left(\ln 0.5\right)$ by this result provides new values of the profile exponents $\alpha_{\text{T}}, \alpha_{\text{P}} ,\alpha_{\text{J}}$ which are now averaged across a shared flux function. + + +$$ +\alpha = \frac{\ln{0.5}}{\ln{\frac{\ln \left(\frac{q_0+(q_{95}-q_0)\times (1-0.5^{\frac{1}{\alpha}})}{q95}\right)}{\ln{\frac{q_0}{q_{95}}}}}} +$$ + +These new calculated profile exponents $(\alpha)$ are what is used in the matrix values below and **not** the standard parabolic profile indexes. + +$$ +b_1 = 1 \ \ b_2 = \alpha_{\text{P}} \ \ b_3 = \alpha_{\text{T}} \ \ b_4 = \alpha_{\text{P}}\alpha_{\text{T}} \\ +b_5 = \epsilon_0^{\frac{1}{2}} \ \ b_6 = \alpha_{\text{P}}\epsilon_0^{\frac{1}{2}} \ \ b_7 = \alpha_{\text{T}}\epsilon_0^{\frac{1}{2}} \ \ b_8 = \alpha_{\text{P}}\alpha_{\text{T}}\epsilon_0^{\frac{1}{2}} \\ +b_9 = \epsilon_0 \ \ b_{10} = \alpha_{\text{P}}\epsilon_0 \ \ b_{11} = \alpha_{\text{T}}\epsilon_0 \ \ b_{12} = \alpha_{\text{P}}\alpha_{\text{T}}\epsilon_0 +$$ + +$$ +a_1 = 1.41\left(1.0 - 0.28 \alpha_{\text{J}}^{\frac{1}{2}}\right)\left(1.0 + \frac{0.12}{Z}\right) \\ +a_2 = 0.36 \left(1.0 - 0.59 \alpha_{\text{J}}^{\frac{1}{2}}\right) \left(1.0 + \frac{0.8}{Z}\right) \\ +a_3 = -0.27 \left(1.0 - 0.47 \alpha_{\text{J}}^{\frac{1}{2}}\right) \left(1.0 + \frac{3.0}{Z}\right) \\ +a_4 = 0.0053 \left(1.0 + \frac{5.0}{Z}\right) \\ +a_5 = -0.93 \left(1.0 - 0.34 \alpha_{\text{J}}^{\frac{1}{2}}\right) \left(1.0 + \frac{0.15}{Z}\right) \\ +a_6 = -0.26 \left(1.0 - 0.57 \alpha_{\text{J}}^{\frac{1}{2}}\right) \left(1.0 - 0.27 Z\right) \\ +a_7 = 0.064 \left(1.0 - 0.6 \alpha_{\text{J}} + 0.15 \alpha_{\text{J}}^2\right) \left(1.0 + \frac{7.6}{Z}\right) \\ +a_8 = -0.0011 \left(1.0 + \frac{9.0}{Z}\right) \\ +a_9 = -0.33 \left(1.0 - \alpha_{\text{J}} + 0.33 \alpha_{\text{J}}^2\right) \\ +a_{10} = -0.26 \left(1.0 - \frac{0.87}{\alpha_{\text{J}}^{\frac{1}{2}}} - 0.16 \alpha_{\text{J}}\right) \\ +a_{11} = -0.14 \left(1.0 - \frac{1.14}{\alpha_{\text{J}}^{\frac{1}{2}}} - 0.45 \alpha_{\text{J}}^{\frac{1}{2}}\right) \\ +a_{12} = -0.0069 \\ +$$ + +Coefficients are obtained by a least squares fit to the 3000 equilibria numerical solutions. Error distribution shows an average error of 3.6% and a maximum error of 20%. These larger errors appear to come from the cases where the temperature profile is approximately equal to the pressure profile. For these cases the density profile is very flat over much of the plasma radius and (as it is forced to be zero at the plasma edge) this +means that it falls off sharply at the plasma edge. + +!!! quote "Excerpt from Wilson[^3]" + + *"The ITER group gives an expression for the bootstrap current[^0] which differs significantly from that + presented here. Their relatively simple expression, + shows no variation with density and temperature profiles. It also does not reproduce the large aspect ratio + scaling with $\beta_{\text{P}}$, and $\epsilon$ in this limit."* + +!!! warning "Flux average profile indexes" + + The implementation of calculating the new profile indexes to be used in the $a$ & $b$ coefficients is not explicitly given in Wilson[^3]. It is an adaption to make the method more aligned with what `PROCESS` calculates. + + +-------------------- + +### Sauter Scaling | `bootstrap_fraction_sauter()` + +Sauter et al.[^4] [^5] provides a formula using the exact Fokker–Planck operator and +without any approximation on the plasma geometry or collisionality. In this way we have been able to accurately determine the neoclassical resistivity and the coefficients for the +bootstrap current which allows one to calculate the bootstrap fraction. + +Is selected by setting `i_bootstrap_current = 4` + +$$ +\left\langle j_{\|} B\right\rangle= \sigma_{\text {neo }}\left\langle E_{\|} B\right\rangle-I(\psi) p(\psi)\left[\mathcal{L}_{31} \frac{\partial \ln n_e}{\partial \psi} \\ ++R_{p e}\left(\mathcal{L}_{31}+\mathcal{L}_{32}\right) \frac{\partial \ln T_e}{\partial \psi}+\left(1-R_{p e}\right) \times\left(1+\frac{\mathcal{L}_{34}}{\mathcal{L}_{31}} \alpha\right) \mathcal{L}_{31} \frac{\partial \ln T_i}{\partial \psi}\right] +$$ + +Note that the above $\left\langle j_{\|} B\right\rangle$ given by Sauter et.al.[^4] gives the component of the current density in the direction of the field – not the toroidal component of the current. The error is second order in the pitch angle, but can be important, especially in the outboard region of a low aspect ratio tokamak, where the pitch angle is large. Moreover this error will always have the effect of overestimating the current. This is accounted for via poloidal correction function implemented as: [`beta_poloidal_sauter()`](#calculate-electron-only-poloidal-beta-correction-beta_poloidal_sauter) and[`beta_poloidal_total_sauter()`](#calculate-ion-and-electron-poloidal-beta-correction-beta_poloidal_total_sauter) + +The correction is of the form: + +$$ +I_{\phi}^{\text{bs}} = 2\pi \int d\psi \frac{q(\psi)}{\langle B^{2}\rangle} \left\langle j_{\|} B\right\rangle +$$ + +where $q(\psi)$ is the safety factor. + +The reconstructed implementation in PROCESS looks as such: + +$$ +\frac{I_{\text{b}}}{I_{\text{P}}} = \sum_2^{\rho_{\text{max}}} \left(2\pi\left[\rho\right]_{-1} \times \left(\left[\rho\right] -\left[\rho\right]_{-1}\right) \right) \times \\ +\left(0.5 \times \left[\mathcal{L}_{31} \frac{\partial \ln n_e}{\partial \psi} ++\left(\mathcal{L}_{31} + \mathcal{L}_{32}\right) \frac{\partial \ln T_e}{\partial \psi} ++ \left(1 + \frac{\mathcal{L}_{34}}{\mathcal{L}_{31}} \alpha\right) \mathcal{L}_{31} \frac{\partial \ln T_i}{\partial \psi}\right] \times \\ + 1 \times 10^6 \times \frac{-B_{0}\left[\rho\right]_{-1}\left[\frac{1}{q}\right]_{-1}}{0.2\pi R_0}\right) +$$ + +In this case square brackets denote array variables equal in length to $\rho_{\text{max}}$ representing the normalized radius elements across the profile. The $-1$ subscript denotes the previous array element in the summation. + +It is not known fully if the $\left(\sigma_{\text {neo }}\left\langle E_{\|} B\right\rangle-I(\psi) p(\psi)\right)$ term is properly implemented into the `PROCESS` version. The $R_{pe}$ value is assumingly taken to be 0.5 as it is stated to approximately to be in Sauter et.al[^4]. + +!!! warning "Validity of the Sauter Bootstrap Scaling" + + In its current state the several base functions called by the Sauter scaling have no reference and cannot be verified. The ad-hoc adaption of the Sauter scaling for use in `PROCESS` is done knowing that `PROCESS` does not calculate flux surfaces across the plasma. + +----------- + +#### Calculate the trapped particle fraction | `_trapped_particle_fraction_sauter()` + +This function calculates the trapped particle fraction $\left(f_t\right)$ used within other key internal Sauter scaling functions. + +$$ +f_t = \frac{1.0 - (1.0-\epsilon_{-1})\sqrt{1.0-\epsilon_{-1}}}{\left(1.0+1.46 \sqrt{\epsilon_{-1}}\right)} +$$ + +$\epsilon$ in this case is the local aspect ratio at that normalised radial point in the profile given by $\epsilon = \rho \left(\frac{a}{R}\right)$. The value of $\rho$ varies from 0 to 1 across the profile. + +The $-1$ subscript in this case refers to the value of the variable in the previous array index value. + +------------- + +#### Calculate electron density coefficient | `_calculate_l31_coefficient()` + +This function calculates and returns the $\mathcal{L}_{31}$ coefficient value for $\frac{\partial \ln n_e}{\partial \psi}$ + +$$ +\mathcal{L}_{31}= F_{31}\left(X=f_{\text {teff }}^{31}\right) \equiv\left(1+\frac{1.4}{Z+1}\right) X-\frac{1.9}{Z+1} X^2+\frac{0.3}{Z+1} X^3 +\frac{0.2}{Z+1} X^4 \\ +f_{\text {teff }}^{31}\left(\nu_{e *}\right)= \frac{f_t}{1+\left(1-0.1 f_t\right) \sqrt{\nu_{e *}}+0.5\left(1-f_t\right) \nu_{e *} / Z} +$$ + +The returned value is $\mathcal{L}_{31} \times$ [`_beta_poloidal_total_sauter()`](#calculate-ion-and-electron-poloidal-beta-correction-beta_poloidal_total_sauter) + +--------------- + +#### Calculate electron temperature coefficient | `_calculate_l31_32_coefficient()` + +This function calculates and returns the $\left(\mathcal{L}_{31}+\mathcal{L}_{32}\right)$ coefficient value for $\frac{\partial \ln T_{\text{e}}}{\partial \psi}$. + +$$ +\begin{align*} +\mathcal{L}_{32} &= F_{32 \_e e}\left(X=f_{\text {teff }}^{32 \_e e}\right)+F_{32 \_e i}\left(Y=f_{\text {teff }}^{32 \_e i}\right) \\ \\ +F_{32 \_e e}(X) &= \frac{0.05+0.62 Z}{Z(1+0.44 Z)}\left(X-X^4\right)+\frac{1}{1+0.22 Z}\left[X^2-X^4\right. \\ +& \quad \left.-1.2\left(X^3-X^4\right)\right]+\frac{1.2}{1+0.5 Z} X^4 \\ +f_{\text {teff }}^{32_{-e e}}\left(\nu_{e *}\right) &= \frac{f_t}{1+0.26\left(1-f_t\right) \sqrt{\nu_{e *}}+0.18\left(1-0.37 f_t\right) \frac{\nu_{e *}}{\sqrt{Z}}} +\end{align*} +$$ + +$$ +F_{32 \_e i}(Y) = -\frac{0.56+1.93 Z}{Z(1+0.44 Z)}\left(Y-Y^4\right)+\frac{4.95}{1+2.48 Z}\left[Y^2-Y^4\right] \\ +-0.55\left(Y^3-Y^4\right)-\frac{1.2}{1+0.5 Z} Y^4 \\ +f_{\text {teff }}^{32 \_e i}\left(\nu_{e *}\right) = \frac{f_t}{1+\left(1+0.6 f_t\right) \sqrt{\nu_{e *}}+0.85\left(1-0.37 f_t\right) \nu_{e *}(1+Z)} +$$ + +The above is added to a call of [`_calculate_l31_coefficient()`](#calculate-electron-density-coefficient-calculate_l31_coefficient). This is then multiplied by [`_beta_poloidal_sauter()`](#calculate-electron-only-poloidal-beta-correction-beta_poloidal_sauter). + +This product above is then multiplied by ([`_beta_poloidal_sauter()`](#calculate-electron-only-poloidal-beta-correction-beta_poloidal_sauter) divided by [`_beta_poloidal_total_sauter()`](#calculate-ion-and-electron-poloidal-beta-correction-beta_poloidal_total_sauter)) + +--------------- + +#### Calculate ion temperature coefficient | `_calculate_l34_alpha_31_coefficient()` + +This function calculates and returns the $\left(1+\frac{\mathcal{L}_{34}}{\mathcal{L}_{31}}\alpha\right)\mathcal{L}_{31}$ coefficient value for $\frac{\partial \ln T_{\text{i}}}{\partial \psi}$. + +$$ +\mathcal{L}_{34}=F_{31}\left(X=f_{\text {teff }}^{34}\right) +$$ + +$$ +f_{\text {teff }}^{34}\left(\nu_{e *}\right)=\frac{f_t}{1+\left(1-0.1 f_t\right) \sqrt{\nu_{e *}}+0.5\left(1-0.5 f_t\right) \nu_{e *} / Z} +$$ + +$$ +\alpha_0=-\frac{1.17\left(1-f_t\right)}{1-0.22 f_t-0.19 f_t^2} +$$ + +$$ +\alpha\left(\nu_{i *}\right)=\left[\frac{\alpha_0+0.25\left(1-f_t^2\right) \sqrt{\nu_{i *}}}{1+0.5 \sqrt{\nu_{i *}}} \\ ++0.315 \nu_{i *}^2 f_t^6\right] \frac{1}{1+0.15 \nu_{i *}^2 f_t^6} +$$ + +The definition of $\alpha\left(\nu_{i *}\right)$ is that found in the erratum paper which changes the value of $-0.315\nu_{i *}^2 f_t^6$ to positive.[^5] + +The return sequence is ([`_beta_poloidal_total_sauter()`](#calculate-ion-and-electron-poloidal-beta-correction-beta_poloidal_total_sauter) - [`_beta_poloidal_sauter()`](#calculate-electron-only-poloidal-beta-correction-beta_poloidal_sauter)) $\times (\mathcal{L}_{34} + \alpha)$ + [`_calculate_l31_coefficient()`](#calculate-electron-density-coefficient-calculate_l31_coefficient) $\times$ (1.0 - [`_beta_poloidal_sauter()`](#calculate-electron-only-poloidal-beta-correction-beta_poloidal_sauter) divided by [`_beta_poloidal_total_sauter()`](#calculate-ion-and-electron-poloidal-beta-correction-beta_poloidal_total_sauter)) + + + +------------- + +#### Calculate the Coulomb logarithm | `_coulomb_logarithm_sauter()` + +$$ +\ln \Lambda = 15.9 -0.5 \times \ln{n_{\text{e}}}+\ln{T_{\text{e}}} +$$ + +----------- + +#### Calculate frequency of electron collisions | `_electron_collisions_sauter()` + +Using the Coulomb logarithm ($\ln \Lambda$) calculated from [`_coulomb_logarithm_sauter()`](#calculate-the-coulomb-logarithm--_coulomb_logarithm_sauter) we get: + +$$ +\nu_{\text{e}} = 670 \times \frac{\ln \Lambda \times n_{\text{e}}}{T_{\text{e}}^{3/2}} +$$ + +------------ + +#### Calculate electron collisionality | `_electron_collisionality_sauter()` + +The origins of the coefficients values are not known, but thought to be derived from a condition of the [Bohm diffusion coefficient](https://en.wikipedia.org/wiki/Bohm_diffusion) + +Using the electron collision frequency ($\nu_{\text{e}}$) calculated from [`_electron_collisions_sauter()`](#calculate-frequency-of-electron-collisions--_electron_collisions_sauter) we get: +$$ +\nu_{\text{e*}} = \frac{1.4 \ R \ \nu_{\text{e}} \ Z_{\text{eff}}}{\left|\frac{1}{q}\epsilon^{3/2}\sqrt{T_{\text{e}}}\times 1.875\times10^7\right|} +$$ + +------------- + +#### Calculate frequency of ion collisions | `_ion_collisions_sauter()` + + +$$ +\nu_{\text{i}} = 320 \times \frac{Z_{\text{eff}}^4n_{\text{i}}}{T_{\text{i}}^{3/2}\sqrt{a_{\text{i}}}} +$$ + +----- + +#### Calculate ion collisionality | `_ion_collisionality_sauter()` + +The origins of the coefficients values are not known, but thought to be derived from a condition of the [Bohm diffusion coefficient](https://en.wikipedia.org/wiki/Bohm_diffusion) + +Using the ion collision frequency ($\nu_{\text{i}}$) calculated from [`_ion_collisions_sauter()`](#calculate-frequency-of-ion-collisions--_ion_collisions_sauter) we get: + +$$ +\nu_{\text{e*}} = \frac{3.2\times10^{-6} \nu_{\text{i}} R}{\left|\left(\frac{1}{q}+0.0001\right)\epsilon^{3/2} \sqrt{\frac{T_{\text{i}}}{a_{\text{i}}}} \right|} +$$ + +---------------- + +#### Calculate electron only poloidal beta correction | `_beta_poloidal_sauter()` + +This function returns an electron only local poloidal beta correction dependant on the array index of the profile. + +If the current index is not equal to the size value (or end of the array) then the following is returned: + +$$ +\frac{1.6\times 10^{-4}\pi R \left(n_{\text{e}}+n_{\text{e-1}}\right)\times \left(T_{\text{e}}+T_{\text{e-1}}\right)}{\left(B_{\text{T}}\rho_{-1}\left|\left(\frac{1}{q}\right)_{-1}+1\times 10^{-4}\right|\right)^2} +$$ + +Otherwise the following is returned: + +$$ +\frac{6.4\times 10^{-4}\pi R \left(n_{\text{e-1}}T_{\text{e-1}}\right)}{\left(B_{\text{T}}\rho_{-1}\left|\left(\frac{1}{q}\right)_{-1}+1\times 10^{-4}\right|\right)^2} +$$ + +The $-1$ subscript in this case refers to the value of the variable in the previous array index value + + + +--------------- + +#### Calculate ion and electron poloidal beta correction | `_beta_poloidal_total_sauter()` + +This function returns the local poloidal beta correction with both electron and ion pressure dependant on the array index of the profile. + +If the current index is not equal to the size value (or end of the array) then the following is returned: + +$$ +\frac{1.6\times 10^{-4}\pi R \left[\left(\left(n_{\text{e}}+n_{\text{e-1}}\right)\times \left(T_{\text{e}}+T_{\text{e-1}}\right)\right)+\left(\left(n_{\text{i}}+n_{\text{i-1}}\right)\times \left(T_{\text{i}}+T_{\text{i-1}}\right)\right)\right]}{\left(B_{\text{T}}\rho_{-1}\left|\left(\frac{1}{q}\right)_{-1}+1\times 10^{-4}\right|\right)^2} +$$ + +Otherwise the following is returned: + +$$ +\frac{6.4\times 10^{-4}\pi R \left[\left(n_{\text{e-1}}T_{\text{e-1}}\right)+\left(n_{\text{i-1}}T_{\text{i-1}}\right)\right]}{\left(B_{\text{T}}\rho_{-1}\left|\left(\frac{1}{q}\right)_{-1}+1\times 10^{-4}\right|\right)^2} +$$ + +The $-1$ subscript in this case refers to the value of the variable in the previous array index value + + +------------------ + +### Sakai Scaling | `bootstrap_fraction_sakai()` + +Is selected by setting `i_bootstrap_current = 5`[^5] + +$$ +f_{\text{BS}} = 10^{0.951 \epsilon - 0.948} \cdot \beta_p^{1.226 \epsilon + 1.584} \cdot l_i^{-0.184\epsilon - 0.282} \cdot \left(\frac{q_{95}}{q_0}\right)^{-0.042 \epsilon - 0.02} \\ +\cdot \alpha_n^{0.13 \epsilon + 0.05} \cdot \alpha_t^{0.502 \epsilon - 0.273} +$$ + +The model includes the toroidal diamagnetic current in the calculation due to the dataset, so `i_diamagnetic_current = 0` can only be used with it + +--------------------- + +## Setting of maximum desirable bootstrap current fraction + +The variable `bootstrap_current_fraction_max` can be set to the value of maximum desirable bootstrap current fraction for a specific design. When optimising if the value of the calculated `bootstrap_current_fraction` for the model selected with `i_bootstrap_current` exceeds this value, then `bootstrap_current_fraction` is set to the value of `bootstrap_current_fraction_max`. + +An error is also raised to the user in the terminal output at the end of the run saying "Bootstrap fraction upper limit enforced". + +## Fixing the bootstrap current fraction + +If the user wants to set the value of the bootstrap current fraction directly then the value can be set by assigning the negative of the desired value to `bootstrap_current_fraction_max`. + + +```txt +>>> IN.DAT + +# Setting a fixed bootstrap current fraction of 80% + +bootstrap_current_fraction_max = -0.8 +``` + + +[^0]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', +[^1]: Nevins, W. M. "Summary report: ITER specialists’ meeting on heating and current drive." ITER-TN-PH-8-4, June 1988. 1988. +[^2]: Keii Gi, Makoto Nakamura, Kenji Tobita, Yasushi Ono, Bootstrap current fraction scaling for a tokamak reactor design study, +Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2014.07.009. +[^3]: Wilson, H.R. (1992). Bootstrap current scaling in tokamaks. Nuclear Fusion, 32(2), pp.257–263. doi:https://doi.org/10.1088/0029-5515/32/2/i05. +[^4]: O. Sauter, C. Angioni, Y. R. Lin-Liu; Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. Phys. Plasmas 1 July 1999; 6 (7): 2834–2839. https://doi.org/10.1063/1.873240 +[^5]: O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: “Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime” [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 +[^6]: 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. +[^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + diff --git a/documentation/proc-pages/physics-models/plasma_current/diamagnetic_current.md b/documentation/proc-pages/physics-models/plasma_current/diamagnetic_current.md new file mode 100644 index 0000000000..a8d2c376da --- /dev/null +++ b/documentation/proc-pages/physics-models/plasma_current/diamagnetic_current.md @@ -0,0 +1,53 @@ + + +## Overview + +The diamagnetic current fraction $f_{\text{dia}}$ is strongly related to $\beta$ and is typically small, +hence it is usually neglected. For high $\beta$ plasmas, such as those at tight +aspect ratio, it should be included and two scalings are offered. If the diamagnetic +current is expected to be above one per cent of the plasma current, a warning +is issued to calculate it. + +$$ +J_{\text{dia}} = \frac{RB_{\theta}}{B^2}\frac{dP}{d\psi} +$$ + +$$ +\frac{I_{\text{dia}}}{I_{\text{p}}} = \frac{\beta_0}{2}\int_0^1 \frac{\rho^2q(a)}{q(\rho)}\frac{d\left(\frac{P}{P(0)}\right)}{dx} dx +$$ + +Where $\rho$ is the normalised radius = $\frac{r}{a}$ + +----------------------------------- + +### No diamagnetic current + +To have it so that the diamagnetic current is not calculated you can set `i_diamagnetic_current = 0` + +------------------------------------ + +### T.Hender fit for ST's: + +This model can be used by setting: `i_diamagnetic_current = 1` + + +$$f_{\text{dia}} = \frac{\beta}{2.8}$$ + + + +------------------------------------ + +### SCENE fit: + +This model can be used by setting: `i_diamagnetic_current = 2` + +This model is based off of 108 equilibria from SCENE. +Overall the equilibria cover: + +- $A$ = 1.6 to 3.2 +- $\beta$ = 0.5% to 26% +- $\frac{P(0)}{\langle P \rangle}$ = 1.8 to 7.2 +- $l_i$(2) = 0.21 to 1.0 +- $\frac{q_{95}}{q(0)}$ = 0.9 to 16 \ \ (i.e. deeply hollow current to very peaked) + +$$f_{\text{dia}} = 0.414 \space \beta \space \left(\frac{0.1 q_{95}}{q_0} + 0.44\right)$$ \ No newline at end of file diff --git a/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md b/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md new file mode 100644 index 0000000000..15c914a43d --- /dev/null +++ b/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md @@ -0,0 +1,7 @@ +Currently in `PROCESS` the inductive current fraction from the CS is not calculated directly but is just equal to ($1 - \mathtt{fvsbrnni}$). Where $\mathtt{fvsbrnni}$ is the sum of the fractions of current driven by non inductive means. + +This calculated fraction (`inductive_current_fraction`) is then used in the `vscalc()` and `burn()` functions to calculate the volt-second requirements and the burn time for a pulsed machine. + +!!! info "Inductive plasma current fraction refactor" + + It is hoped for the near future to have a more engineering based calculation of the fraction of the plasma current driven by the central solenoid. This would hopefully allow the setting of required ramp and flat-top times for which a given inductive current fraction can be given based on the operational performance and margin in the central solenoid. \ No newline at end of file diff --git "a/documentation/proc-pages/physics-models/plasma_current/pfirsch_schl\303\274ter_current_drive.md" "b/documentation/proc-pages/physics-models/plasma_current/pfirsch_schl\303\274ter_current_drive.md" new file mode 100644 index 0000000000..8ef50f3ccd --- /dev/null +++ "b/documentation/proc-pages/physics-models/plasma_current/pfirsch_schl\303\274ter_current_drive.md" @@ -0,0 +1,34 @@ +## Overview + +A similar SCENE scaling like that for the diamagnetic current is available for the Pfirsch-Schlüter current fraction $f_{\text{PS}}$. This is +typically smaller than the diamagnetic current, but is negative. + +There is no ability to input the diamagnetic and Pfirsch-Schlüter current +directly. In this case, it is recommended to turn off these two scalings +and to use the method of fixing the bootstrap current fraction. + +-------------- + +### No Pfirsch-Schlüter current + +To have it so that the Pfirsch-Schlüter current is not calculated you can set `i_pfirsch_schluter_current = 0` + +------------------ + +### SCENE Fit: + +This model can be used by setting: `i_pfirsch_schluter_current = 2` + +This model is based off of 108 equilibria from SCENE. +Overall the equilibria cover: + +- $A$ = 1.6 to 3.2 +- $\beta$ = 0.5% to 26% +- $\frac{P(0)}{\langle P \rangle}$ = 1.8 to 7.2 +- $l_i$(2) = 0.21 to 1.0 +- $\frac{q_{95}}{q(0)}$ = 0.9 to 16 \ \ (i.e. deeply hollow current to very peaked) + +This fit is derived from the same dataset used to derive the diamagnetic current fraction for [`i_diamagnetic_current = 2`](diamagnetic_current.md#scene-fit) + +$$ f_{\text{PS}} = -0.09 \beta $$ + diff --git a/documentation/proc-pages/physics-models/plasma_current/plasma_current.md b/documentation/proc-pages/physics-models/plasma_current/plasma_current.md new file mode 100644 index 0000000000..957e954ace --- /dev/null +++ b/documentation/proc-pages/physics-models/plasma_current/plasma_current.md @@ -0,0 +1,716 @@ +# Plasma Current + +## Overview + +In tokamaks, the plasma current ($I_{\text{p}}$) plays a crucial role in confining the plasma and maintaining the stability of the magnetic fields. The plasma current generates a poloidal magnetic field, which, when combined with the toroidal magnetic field, creates a helical magnetic field structure. + +The plasma current is typically driven by inductive means, where a central solenoid induces a current in the plasma, similar to the secondary winding of a transformer. However, non-inductive methods such as neutral beam injection and radio frequency (RF) heating are also employed, especially for steady-state operation. + +More information about the different heating and current drive models can be found [here](../../eng-models/heating_and_current_drive/heating-and-current-drive.md). + +In `PROCESS` there are several inductive and non-inductive contributions to the overall plasma current. The fractions of these contributions are checked to add up to 1.0. + +$$ +1.0 = [\underbrace{\overbrace{I_{\text{F,Boot}} + I_{\text{F,Dia}} + I_{\text{F,PS}}}^{\text{Generated by plasma}} + I_{\text{F,CD}}}_{\text{Non-inductive}}] + \underbrace{I_{\text{F,CS}}}_{\text{Inductive}} +$$ + +Here $I_{\text{F,Boot}}$ is the [bootstrap current fraction](bootstrap_current.md), $I_{\text{F,Dia}}$ is the [diamagnetic current fraction](diamagnetic_current.md), $I_{\text{F,PS}}$ is the [Pfirsch-Schlüter current fraction](pfirsch_schlüter_current_drive.md), $I_{\text{F,CD}}$ is the external [heating and current drive fraction](../../eng-models/heating_and_current_drive/heating-and-current-drive.md) and $I_{\text{F,CS}}$ is the [current fraction driven by the central solenoid](inductive_plasma_current.md). + + +-------------------------- + +### Derivation of plasma current relation + +Starting with the definition of $q$ in an axisymmetric equilibria: + +$$ +q = \frac{\Delta \phi}{2\pi} +$$ + +This is defined as the rate of returning to the same position in the poloidal plane after a change of toroidal angle $\Delta \phi$. + +In order to calculate the value of $q$ it is necessary to use the equation of the field line: + +$$ +\frac{R d\phi}{ds} = \frac{B_{\text{T}}}{B_{\text{p}}} +$$ + +where $ds$ is the distance moved in the poloidal direction while moving through a toroidal angle $d\phi$ and $B_{\text{p}}$ and $B_{\text{T}}$ are the poloidal and toroidal magnetic fields. + +Re-arranging the two equations above and substituting for $q$ we get: + +$$ +q = \frac{1}{2\pi} \oint \frac{1}{R}\frac{B_{\text{T}}}{B_{\text{p}}} ds +$$ + +where the integral is carried out over a single poloidal circuit around the flux surface. + +Assuming a large aspect ratio tokamak that has a circular plasma cross-section the integral simplifies to: + +$$ +q = \frac{r}{R_0}\frac{B_{\text{T}}}{B_{\text{p}}} +$$ + +where $r$ is the minor radius of the flux surface and $R_0$ is the major radius of the tokamak. Above is done assuming $ds = 2\pi r$. + + +From above, **assuming the toroidal field to be a constant radially**, the only variation of $q$ across the flux surfaces is caused solely by the poloidal magnetic field produced by the plasma current. Again assuming a large aspect ratio machine with circular plasma and flux surface cross sections, the toroidal current density profile $j(r)$, dictates $q$. + +Writing the plasma current via [Amperes's law](https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law): + +$$ +2\pi r B_{\text{p}} = \mu_0 I(r) +$$ + +where the current inside $r$ is: + +$$ +I(r) = 2\pi \int_0^r j(r^{\prime})r^{\prime} dr^{\prime} +$$ + +Substituting into the value of $q$ above: + +$$ +q(r) = \frac{2\pi r^2 B_{\text{T}}}{\mu_0 I(r) R} +$$ + +Taking the limit of $r$ to be the plasma minor radius, $a$: + +$$ +q_a = \frac{2\pi a^2 B_{\text{T}}}{\mu_0 I R} +$$ + +Re-arranging for $I$ we get a function for the total plasma current: + +$$ +I = \frac{2\pi a^2B_{\text{T}}}{\mu_0q_{95} R} +$$ + +Instead of $q_a$, $q_{95}$ is used as in plasma configurations with divertors the poloidal magnetic field has a null at the X-point. Thus as the separatrix is approached, $q$ tends to infinity. + +----------------- + +## Plasma Current Calculation | `calculate_plasma_current()` + +This function calculates the plasma current shaping factor ($f_q$), plasma current ($I_{\text{p}}$), qstar ($q^*$), normalized beta ($\beta_{\text{N}}$) and then poloidal field and the profile settings for $\mathtt{alphaj}$ ($\alpha_J$) and $\mathtt{rli}$ ($l_{\mathtt{i}}$). This is done in 5 separate steps which are shown in the following numbered sections. + + +$$\begin{aligned} +I_{\text{p}} = f_q \frac{2\pi}{\mu_0} \frac{a^2 B_{\text{T}}}{R \ q_{95}} +\end{aligned}$$ + + +The factor $f_q$ makes allowance for toroidal effects and plasma shaping (elongation and +triangularity). Several formulae for this factor are available depending on the value of +the switch `i_plasma_current`, as follows: + +--------------- + +### 1. Calculate plasma current shaping function $f_q$ + +------------ + +#### Peng analytic fit | `calculate_current_coefficient_peng()` + +Switch value: `i_plasma_current = 1` + +The formula for calculating `fq` is: + +$$f_q = \left(\frac{{1.22 - 0.68 \epsilon}}{{(1.0 - \epsilon^2)^2}}\right) \mathtt{{sf}}^2$$ + +Where $\epsilon$ is the inverse [aspect ratio](../plasma_geometry.md) ($\mathtt{eps}$) and $\mathtt{sf}$ is the shaping factor calculated in the [poloidal perimeter](../plasma_geometry.md#poloidal-perimeter) function in `plasma_geometry.py` + +----------- + +#### STAR, Peng double null divertor scaling | `calculate_plasma_current_peng()` + +Switch value: `i_plasma_current = 2` [^3] [^4] + +This is currently the only scaling in which the calculated plasma current does not follow the form of that show above. The bounds of applicability is for aspect ratios $\le 3$. + +The plasma current is given by: + +$$ +I_{\text{p}} = \frac{5a B_{\text{T}}\kappa}{2\pi^2 \bar{q}}(F_1+F_2)\left(\frac{\arcsin{E_1}}{E_1}+\frac{\arcsin{E_2}}{E_2}\right) +$$ + +The values of $F_1$, $F_2$, $d_1$ & $d_2$ are first calculated from the [`_plascar_bpol()`](#_plasc_bpol) function. + +The values of $E_1$ & $E_2$ are then calculated such as + +$$ +E_1 = \frac{2\kappa}{d_1(1+\delta)} +$$ + +$$ +E_2 = \frac{2\kappa}{d_2(1.0 - \delta)} +$$ + +$I_{\text{p}}$ from above is then calculated from these values +$\bar{q}$ is defined as the average safety factor[^4]. +In the original STAR code model $\bar{q}$ is defined as: + +$$ +\bar{q} = \bar{q}_0(1+2.6\epsilon^{2.8}) +$$ + +The STAR code document states that the $\bar{q}_0$ is normally equal to 3.0 +In PROCESS this is implemented differently as the function: + +$$ +\bar{q} = q_{95} \times 1.3(1-\epsilon)^{0.6} +$$ + +This is to allow the use of the available $q_{95}$ parameter as **the origin and definition of $\bar{q}_0$ is not fully known.** + +The coefficient values $q_{95}$ and $\bar{q}_0$ for the PROCESS and STAR code implementations can be experimented with in the graph below: + + + + + + + + +My Bokeh Plot + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + +---------------- + +#### Simple ITER scaling + +Switch value: `i_plasma_current = 3` [^5] + +The simple cylindrical case is assumed so: + +$$ +f_q = 1 +$$ + +----------------- +#### ITER IPDG89 scaling | `calculate_current_coefficient_ipdg89()` + +Switch value: `i_plasma_current = 4`[^5] [^7] + +The formula for calculating `fq` is: + +$$ +f_q = \left(\frac{{0.5 \cdot (1.17 - 0.65 \cdot \epsilon)}}{{(1.0 - \epsilon^2)^2}} \cdot \left(1.0 + \kappa_{95}^2 \cdot \left(1.0 + 2.0 \cdot \delta_{95}^2 - 1.2 \cdot \delta_{95}^3\right)\right)\right) +$$ + + + + + +-------------- + +#### Todd empirical scaling, I | `calculate_current_coefficient_todd()` + +Switch value: `i_plasma_current = 5`[^6] [^7] + +The formula for calculating `fq` is: + +$$ +F_{T1} | f_q = \left( + (1.0 + 2.0\epsilon^2) + \cdot 0.5 + \cdot (1.0 + \kappa_{95}^2) + \cdot ( + 1.24 + - 0.54 \cdot \kappa_{95} + + 0.3 \cdot (\kappa_{95}^2 + \delta_{95}^2) + + 0.125 \cdot \delta_{95} + ) + \right) +$$ + + +------------------- + +#### Todd empirical scaling, II + +Switch value: `i_plasma_current = 6` [^6] [^7] + +This function is similar to the previous [Todd scaling](#todd-empirical-scaling-i) except it is multiplied by a new elongation dependant term + + +$$ +f_q = F_{T1} \times \left(1+\left(\kappa-1.2\right)^3\right) +$$ + +------------------ + + +#### Connor-Hastie model | `calculate_current_coefficient_hastie()` + +Switch value: `i_plasma_current = 7` [^7] [^8] + +Asymptotically correct in the range of: + +- $\epsilon \ll 1$ +- $\delta \ll 1$ +- $(\kappa - 1) \ll 1$ + +$$ +f_q = \left(\frac{(\kappa+1)^2}{2}\right)\left(1+\left(\frac{\kappa+1}{2}\right)^2\epsilon^2+\frac{1}{2}\Delta^{\prime 2}+2\frac{\Delta}{R_0} \\ ++ \frac{1}{2}\left(E^{\prime 2}+\frac{E^2}{r^2}\right)+\frac{1}{2}\left(T^{\prime 2}+\frac{4T^2}{r^2}\right)\right) +$$ + +where: + +$$ +\frac{T}{r} = \frac{\kappa \delta}{(\kappa+1)^2} \ \ \ ; \ \ \ \frac{E}{r} = \frac{\kappa-1}{\kappa+1} +$$ + +$$ +T^{\prime} = 2 \left(\frac{1+\lambda}{1+\frac{\lambda}{2}}\right)\frac{T}{r} \ \ \ ; \ \ \ E^{\prime} = 2 \left(\frac{1+\lambda}{1+\frac{\lambda}{3}}\right)\frac{E}{r} +$$ + +$$ +\Delta^{\prime} = \frac{\kappa+1}{2}\epsilon\frac{l_i}{2}+\frac{\beta_0 (1+\lambda)^2}{(\frac{\kappa+1}{2}\epsilon)(1+\nu)} +$$ + +with + +$$ +l_i = \frac{1+\lambda}{\lambda}\left(\frac{1+\lambda}{\lambda}\text{ln}(1-\lambda)-1\right) +$$ + +$$ +\frac{\Delta}{R_0} = \frac{\beta_0}{6}\left[1+\frac{5}{6}\lambda+\frac{1}{4}\lambda^2\right]+\left(\frac{\kappa+1}{2}\epsilon\right)^2 \frac{1}{8}\left(1-\frac{\lambda^2}{3}\right) +$$ + +The parameters $\lambda$ and $\nu$ characterise the current profile $J_{\phi} = \frac{j_0}{(1+\lambda r^2)^2}$ and pressure profile $p = p_0(1-r^2)^{\nu}$. The pressure profile used by Connor-Hastie is still of the normal parabolic type. + + +!!! warning Assumed current profiles + + $\lambda$ in this case is treated as the current profile index within the Connor-Hastie model. This profile is not the same as the standard [parabolic profile](../profiles/plasma_profiles.md`) used within PROCESS for the current which uses the $\alpha_{\text{J}}$ index. The difference between these two assumed current profiles can be experimented with below. + + + + + + + + +My Bokeh Plot + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + +--------------- + +#### Sauter model, allows negative triangularity | `calculate_current_coefficient_sauter()` + +Switch value: `i_plasma_current = 8`[^9] + +Assumes zero squareness with the $w_{07}$ parameter ($w_{07}$ = 1). + +The values for $\kappa$ and $\delta$ is taken at the separatrix and is not the 95% values. + +!!! note "$w_{07}$ setting" + + There is currently no parameter to set in the input file in order to change $w_{07}$ directly. This can be done directly through the code in `physics.py` + +$$ +f_q = \frac{4.1 \times 10^6}{\frac{2\pi}{\mu_0}}[1+1.2(\kappa-1)+0.56(\kappa-1)^2] \\ +\times (1+0.09\delta +0.16\delta^2)\frac{1+0.45\delta \epsilon}{1-0.74\epsilon}[1+0.55(w_{07}-1)] +$$ + +---------------------- + +#### FIESTA for ST's | `calculate_current_coefficient_fiesta()` + +Switch value: `i_plasma_current = 9`[^10] + +Assumptions: + + - D-shaped plasmas with $A < 3$: + - X-pont values of $\kappa$ & $\delta$ used + +A set of equilibria from FIESTA were created and compared to the calculations from the [Peng double null divertor scaling](plasma_current.md#peng-double-null-divertor-scaling-st) For the low elongation equilibria, the calculated values for +the plasma current were close to those from FIESTA, however moving to +higher elongations causes an underestimate by up to 20%. + +Given the parameter dependencies a new plasma current relation based on fits to FIESTA +equilibria was created. It showed no bias with any parameter fitted and that the fit is accurate to 10%. +The linear relation between these and the 95% values expressed in does not hold at high values of elongation and triangularity as per the [`ishape = 0`](../plasma_geometry.md#plasma-geometry-parameters-geomty) relation. + + +$$ +f_q = 0.538 (1.0 + 2.440 \epsilon^{2.736}) \kappa^{2.154}\delta^{0.060} +$$ + +----------------------------- + +### 2. Calculate the cylindrical safety factor + +In reactor design the total plasma current $I_{\text{p}}$, is limited by considerations of tearing mode stability, and major disruptions. These considerations place a limit on the safety factor at the plasma boundary, $q_{\psi}$. In a cylindrical plasma model there is a simple relation between this value, $q_{\text{cyl}}$, and the plasma current, given by + +$$ +q_{\psi} = q_{\text{cyl}} \equiv \frac{5a^2B_0}{R_0I_p} +$$ + +where the minor and major radii are expressed in metres, the toroidal magnetic field $B_0$, in Tesla, and the plasma current $I_{\text{p}}$, in megaamps. Thus a stability requirement of the form $q_{\psi} > q_{\text{crit}}$ (usually considered to be ~ 2.0) places a limit on $I_{\text{p}}$. +The effect of toroidicity and of shaping of the plasma cross section is to modify the relation above between $q_{\psi}$ and the plasma current, $I_{\text{p}}$. In general this relation may be written in the form: + + +$$ +q_{\psi} = \frac{5a^2B_0}{R_0I_p}F\left(\epsilon,\kappa,\delta; \text{profiles}\right) +$$ + +It is not possible to derive a general analytic expression for $F$, so it has been customary to employ an empirical, or semi-empirical expression involving $\epsilon$,$\kappa$ and $\delta$, chosen to fit data taken from numerical solutions of the Grad Shafranov equation. + +The cylindrical equivalent $q$ definition used currently is the one given below from IPDG89[^5] + +$$ +q_* = \frac{5 \times 10^6a^2B_T}{RI_{\text{p}}}\frac{(1+\kappa^2(1+2\delta^2-1.2\delta^3))}{2} +$$ + + +!!! note "$q_{\text{cyl}}$ definition and use in PROCESS" + + Though this seems to not be the definition of the cylindrical safety factor due to the elongation and triangularity scaling. + Zohm[^12], Wesson[^13], and AEA FUS 172[^7] all say that cylindrical safety factor is just the first fraction. + And if we assume a cylindrical plasma where $\kappa = 1$ and $\delta = 0$ then the second term reduces to 1. + +-------------- + +### 3. Calculate the normalized beta + +The total normalized beta is calculated as per: + +$$ +\beta_N = \beta\frac{1\times10^8 a B_{\text{T}}}{I_{\text{P}}} +$$ + + +----------------- + +### 4. Plasma Current Poloidal Field + +For calculating the poloidal magnetic field created due to the presence of the plasma current, [Ampere's law](https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law) can be used. In this case the poloidal field is simply returned as: + +$$ +B_{\text{p}} = \frac{\mu_0 I_{\text{p}}}{\mathtt{pperim}} +$$ + +Where `pperim` is the plasma poloidal perimeter calculated [here](../plasma_geometry.md#poloidal-perimeter). + +In the case of using the Peng double null scaling ([`i_plasma_current = 2`](plasma_current.md#star-peng-double-null-divertor-scaling-st)), the values $F_1$ and $F_2$ are calculated from [_plasc_bpol](plasma_current.md#_plasc_bpol) and used to calculated the poloidal field from the toroidal field as per: + +$$ +B_{\text{p}} = B_{\text{T}}\frac{F_1 + F_2}{2\pi \overline{q}} +$$ + +------------ + +### 5. Plasma Current Profile Consistency + +A limited degree of self-consistency between the plasma current profile and other parameters can be +enforced by setting switch `iprofile = 1`. This sets the current +profile peaking factor $\alpha_J$ (`alphaj`), the normalised internal inductance $l_i$ (`rli`) and beta limit $g$-factor (`dnbeta`) using the +safety factor on axis `q0` and the cylindrical safety factor $q_*$ (`qstar`): + +$$\begin{aligned} +\alpha_J = \frac{q*}{q_0} - 1 +\end{aligned}$$ + +$$\begin{aligned} +l_i = \rm{ln}(1.65+0.89\alpha_J) +\end{aligned}$$ + +$$\begin{aligned} +g = 4 l_i +\end{aligned}$$ + +It is recommended that current scaling law `i_plasma_current = 4` is used if `iprofile = 1`. +This relation is only applicable to large aspect ratio tokamaks. + +For spherical tokamaks, the normalised internal inductance can be set from the elongation using `iprofile = 4` or `iprofile = 5` or `iprofile = 6`[^11]: + +$$\begin{aligned} +l_i = 3.4 - \kappa_x +\end{aligned}$$ + +Further description of `iprofile` is given in [Beta Limit](../plasma_beta.md). + +----------------------- + +## _plasc_bpol + +This internal function is used to calculate the plasma shape and poloidal coefficients for calculating the plasma current in the [Peng double null scaling from the STAR code](plasma_current.md#star-peng-double-null-divertor-scaling-st). If this scaling is selected the coefficients are also used to calculate the [poloidal field from the plasma current](plasma_current.md#plasma-current-poloidal-field). + + +where if $A < \frac{\kappa^2}{(1+\delta)}+\delta$: + +$$ +F_1 = f_1\left(g-h_1\ln \left(\frac{1+y_1}{1-y_1}\right)\right) +$$ + +$$ +f_1 = \frac{d_1(1+\delta)\epsilon}{(1+\epsilon)(c_1\epsilon -1)} +$$ + +$$ +h_1 = \frac{1+(1-c_1)\frac{\epsilon}{2}}{(1+\epsilon)(c_1 \epsilon -1)} +$$ + +$$ +y_1 = \frac{\sqrt{c_1\epsilon-1}}{1+\epsilon}\frac{1+\delta}{\kappa} +$$ + +and if $A > \frac{\kappa^2}{(1+\delta)}+\delta$: + +$$ +F_1 = f_1\left(-g +2h_1 \arctan(y_1)\right) +$$ + +$$ +f_1 = -\frac{d_1(1+\delta)\epsilon}{(1+\epsilon)(c_1\epsilon -1)} +$$ + +$$ +h_1 = \frac{1+(1-c_1)\frac{\epsilon}{2}}{(1+\epsilon)(1-c_1 \epsilon)} +$$ + +$$ +y_1 = \frac{\sqrt{1-c_1\epsilon}}{1+\epsilon}\frac{1+\delta}{\kappa} +$$ + +where both conditions share: + +$$ +g = \frac{\epsilon \kappa}{(1-\epsilon \delta)} +$$ + +$$ +F_2 = f_2\left(g +2h_2 \arctan(y_2)\right) +$$ + +$$ +f_2 = -\frac{d_2(1-\delta)\epsilon}{(1-\epsilon)(c_2\epsilon -1)} +$$ + +$$ +h_2 = \frac{1+(c_2-1)\frac{\epsilon}{2}}{(1-\epsilon)(c_2 \epsilon+1)} +$$ + +$$ +E_1 = \frac{2\kappa}{d_1(1+\delta)} ; \ \ E_2 = \frac{2\kappa}{d_2(1-\delta)} +$$ + +$$ +c_1 = \frac{\kappa^2}{(1+\delta)}+\delta \ ; \ \ c_2 = \frac{\kappa^2}{(1-\delta)}+\delta +$$ + +$$ +d_1 = \left(\frac{\kappa}{1+\delta}\right)^2+1 \ ; \ \ d_2 = \left(\frac{\kappa}{1-\delta}\right)^2+1 +$$ + +$$ +y_2 = \frac{\sqrt{c_2\epsilon+1}}{1-\epsilon}\frac{1-\delta}{\kappa} +$$ + +----------------------- + +## Key Constraints + +----------------------- + +### Plasma current ramp-up time lower limit + +This constraint can be activated by stating `icc = 41` in the input file. + +The value of `tohsm` can be set to the required minimum plasma current ramp up time at the start of a pulse. The scaling value `ftohs` can be varied also + +The calculated plasma current ramp up time `tohs` is dictated by the [pulsed plant operation configuration](../pulsed-plant.md). + +This constraint will ensure that the value of `tohs` is always greater than or equal to `tohsm` + +-------------------- + +### Plasma and Rod current fraction upper limit + +This constraint can be activated by stating `icc = 46` in the input file. + +The constraint condition is[^14]: + +$$ +\frac{I_{\text{p}}}{I_{\text{cp}}} < 1.0 + 4.91\left(\epsilon - 0.62\right) +$$ + +In this case $I_{\text{cp}}$ is the total current going up the centrepost in a spherical tokamak. +This constraint was initially though to prevent instabilities and act as a guideline to limit power dissipation when generating new designs. The scaling value for the constraint, `fipir` can be varied also. + + The origins of the relation should be seen in early spherical tokamak papers not yet referenced here. + + + +[^3]: Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), 1729–1738. https://doi.org/10.13182/FST92-A29971 +[^4]: J.D. Galambos, 'STAR Code : Spherical Tokamak Analysis and Reactor Code', +Unpublished internal Oak Ridge document. +[^5]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', +[^6]: D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 +[^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 +[^8]: J.W.Connor and R.J.Hastie, Culham Lab Report CLM-M106 (1985). https://scientific-publications.ukaea.uk/wp-content/uploads/CLM-M106-1.pdf +[^9]: O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, ISSN 0920-3796, +https://doi.org/10.1016/j.fusengdes.2016.04.033. +[^10]: Stuart I. Muldrew, Hanni Lux, Geof Cunningham, Tim C. Hender, Sebastien Kahn, Peter J. Knight, Bhavin Patel, Garry M. Voss, Howard R. Wilson, “PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. +[^11]: J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, doi: https://doi.org/10.1088/0029-5515/56/10/106023. +[^12]: Zohm Hartmut, 2019, "On the size of tokamak fusion power plants", Phil. Trans. R. Soc. A.37720170437 +http://doi.org/10.1098/rsta.2017.0437 +[^13]: Wesson, J. and Campbell, D. J. (2004) Tokamaks. Clarendon Press (International series of monographs on physics). Available at: https://books.google.co.uk/books?id=iPlAwZI6HIYC. +[^14]: Stuart I. Muldrew, Hanni Lux, Geof Cunningham, Tim C. Hender, Sebastien Kahn, Peter J. Knight, Bhavin Patel, Garry M. Voss, Howard R. Wilson,“PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, Volume 154, 2020, +111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. + + + + diff --git a/documentation/proc-pages/physics-models/plasma_geometry.md b/documentation/proc-pages/physics-models/plasma_geometry.md index a290f8b5ca..a0229dea46 100644 --- a/documentation/proc-pages/physics-models/plasma_geometry.md +++ b/documentation/proc-pages/physics-models/plasma_geometry.md @@ -377,7 +377,7 @@ $$ \mathtt{pperim} = 2.0 \times (\mathtt{xo} \times \mathtt{thetao} + \mathtt{xi} \times \mathtt{thetai}) $$ -The shaping factor for `icurr=1` is also claculated here: +The shaping factor for `i_plasma_current = 1` is also calculated here: $$ \mathtt{sf} = \frac{\mathtt{pperim}}{ 2.0\pi a} diff --git a/documentation/proc-pages/stylesheets/vardes.css b/documentation/proc-pages/stylesheets/vardes.css index 39787d98ed..9d4bb49acc 100644 --- a/documentation/proc-pages/stylesheets/vardes.css +++ b/documentation/proc-pages/stylesheets/vardes.css @@ -19,6 +19,11 @@ td { padding: 10px; border: 1px solid; } +hr { + border: 0; + height: 0.1px; /* Adjust the height to change the thickness */ + background: rgba(128, 128, 128, 0.25); /* Adjust the color and transparency as needed */ +} html { font-size: 1.3rem; diff --git a/examples/data/csv_output_large_tokamak_MFILE.DAT b/examples/data/csv_output_large_tokamak_MFILE.DAT index c7cb3dbc66..ac73645d0c 100644 --- a/examples/data/csv_output_large_tokamak_MFILE.DAT +++ b/examples/data/csv_output_large_tokamak_MFILE.DAT @@ -331,8 +331,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 3.8398E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.1738E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 1.8882E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.6699E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.6699E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.8387E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.1898E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -7.3601E-01 @@ -520,7 +520,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -3.2449E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 4.1436E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.9627E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 4.0562E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 3.0010E+03 OP @@ -535,7 +535,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 7.5000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 2.1293E-01 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.5000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.5000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 2.0990E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 9.8586E+03 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.6299E-02 OP @@ -544,10 +544,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1436E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 5.9038E-04 - Inductive_fraction______________________________________________________ (facoh)_______________________ 5.8505E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 5.9038E-04 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.8505E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.1495E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 7.5213E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 2.0000E+02 @@ -1568,13 +1568,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -1649,7 +1649,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/examples/data/large_tokamak_1_MFILE.DAT b/examples/data/large_tokamak_1_MFILE.DAT index e34aae237f..b27363afbd 100644 --- a/examples/data/large_tokamak_1_MFILE.DAT +++ b/examples/data/large_tokamak_1_MFILE.DAT @@ -332,8 +332,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 3.8398E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.1738E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 1.8882E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.6631E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.6631E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9805E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2275E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -7.4373E-01 @@ -516,7 +516,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -3.0283E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 4.1920E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.7450E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 3.9755E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 3.0619E+03 OP @@ -533,7 +533,7 @@ Ratio_of_power_for_flat-top_to_start-up_(MW)____________________________ (startupratio)________________ 1.0000E+00 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 7.5000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 5.1429E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.5000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.5000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 2.0060E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 2.3908E+05 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.6488E-02 OP @@ -542,10 +542,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (facoh)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 2.0000E+02 @@ -1562,13 +1562,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -1643,7 +1643,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/examples/data/large_tokamak_2_MFILE.DAT b/examples/data/large_tokamak_2_MFILE.DAT index 0b53497232..422e223a73 100644 --- a/examples/data/large_tokamak_2_MFILE.DAT +++ b/examples/data/large_tokamak_2_MFILE.DAT @@ -332,8 +332,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 3.8398E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.1738E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 1.8882E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.6631E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.6631E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9805E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2275E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -7.4373E-01 @@ -516,7 +516,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -3.0283E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 4.1920E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.7450E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 3.9755E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 3.0619E+03 OP @@ -533,7 +533,7 @@ Ratio_of_power_for_flat-top_to_start-up_(MW)____________________________ (startupratio)________________ 1.0000E+00 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 7.5000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 5.1429E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.5000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.5000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 2.0060E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 2.3908E+05 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.6488E-02 OP @@ -542,10 +542,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (facoh)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 2.0000E+02 @@ -1562,13 +1562,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -1643,7 +1643,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/examples/data/large_tokamak_3_MFILE.DAT b/examples/data/large_tokamak_3_MFILE.DAT index db22bbd540..53207efed5 100644 --- a/examples/data/large_tokamak_3_MFILE.DAT +++ b/examples/data/large_tokamak_3_MFILE.DAT @@ -332,8 +332,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 3.8398E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.1738E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 1.8882E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.6631E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.6631E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9805E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2275E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -7.4373E-01 @@ -516,7 +516,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -3.0283E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 4.1920E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.7450E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 3.9755E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 3.0619E+03 OP @@ -533,7 +533,7 @@ Ratio_of_power_for_flat-top_to_start-up_(MW)____________________________ (startupratio)________________ 1.0000E+00 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 7.5000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 5.1429E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.5000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.5000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 2.0060E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 2.3908E+05 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.6488E-02 OP @@ -542,10 +542,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (facoh)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 2.0000E+02 @@ -1562,13 +1562,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -1643,7 +1643,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/examples/data/large_tokamak_4_MFILE.DAT b/examples/data/large_tokamak_4_MFILE.DAT index 723ffcc459..c2cb109cf9 100644 --- a/examples/data/large_tokamak_4_MFILE.DAT +++ b/examples/data/large_tokamak_4_MFILE.DAT @@ -332,8 +332,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 3.8398E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.1738E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 1.8882E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.6631E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.6631E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9805E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2275E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -7.4373E-01 @@ -516,7 +516,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -3.0283E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 4.1920E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.7450E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 3.9755E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 3.0619E+03 OP @@ -533,7 +533,7 @@ Ratio_of_power_for_flat-top_to_start-up_(MW)____________________________ (startupratio)________________ 1.0000E+00 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 7.5000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 5.1429E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.5000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.5000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 2.0060E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 2.3908E+05 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.6488E-02 OP @@ -542,10 +542,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (facoh)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 2.0000E+02 @@ -1562,13 +1562,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -1643,7 +1643,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/examples/data/large_tokamak_IN.DAT b/examples/data/large_tokamak_IN.DAT index a1298b175f..d686f0f0c0 100644 --- a/examples/data/large_tokamak_IN.DAT +++ b/examples/data/large_tokamak_IN.DAT @@ -372,13 +372,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -453,7 +453,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/examples/data/scan_MFILE.DAT b/examples/data/scan_MFILE.DAT index 12406c6309..85850a9f25 100644 --- a/examples/data/scan_MFILE.DAT +++ b/examples/data/scan_MFILE.DAT @@ -187,8 +187,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -373,7 +373,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -388,7 +388,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -397,10 +397,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -1182,8 +1182,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -1368,7 +1368,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -1383,7 +1383,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -1392,10 +1392,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -2177,8 +2177,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -2363,7 +2363,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -2378,7 +2378,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -2387,10 +2387,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -3172,8 +3172,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -3358,7 +3358,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -3373,7 +3373,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -3382,10 +3382,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -4167,8 +4167,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -4353,7 +4353,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -4368,7 +4368,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -4377,10 +4377,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -5162,8 +5162,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -5348,7 +5348,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -5363,7 +5363,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -5372,10 +5372,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -6157,8 +6157,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -6343,7 +6343,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -6358,7 +6358,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -6367,10 +6367,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -7152,8 +7152,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -7338,7 +7338,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -7353,7 +7353,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -7362,10 +7362,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -8147,8 +8147,8 @@ Plasma_cross-sectional_area_(m2)________________________________________ (xarea)_______________________ 4.7287E+01 Plasma_surface_area_(m2)________________________________________________ (sarea)_______________________ 1.4956E+03 OP Plasma_volume_(m3)______________________________________________________ (vol)_________________________ 2.6696E+03 OP - Plasma_current_scaling_law_used_________________________________________ (icurr)_______________________ 4 - Plasma_current_(MA)_____________________________________________________ (plascur/1D6)_________________ 1.8078E+01 + Plasma_current_scaling_law_used_________________________________________ (i_plasma_current)_______________________ 4 + Plasma_current_(MA)_____________________________________________________ (plasma_current_MA)_________________ 1.8078E+01 Current_density_profile_factor__________________________________________ (alphaj)______________________ 1.9256E+00 Plasma_internal_inductance,_li__________________________________________ (rli)_________________________ 1.2131E+00 Vertical_field_at_plasma_(T)____________________________________________ (bvert)_______________________ -6.9441E-01 @@ -8333,7 +8333,7 @@ Pfirsch-Schlueter_fraction_(SCENE)______________________________________ (pscf_scene)__________________ -2.7628E-03 Bootstrap_fraction_(enforced)___________________________________________ (bootipf.)____________________ 3.8470E-01 Diamagnetic_fraction_(enforced)_________________________________________ (diaipf.)_____________________ 0.0000E+00 - Pfirsch-Schlueter_fraction_(enforced)___________________________________ (psipf.)______________________ 0.0000E+00 + Pfirsch-Schlueter_fraction_(enforced)___________________________________ (ps_current_fraction.)______________________ 0.0000E+00 Loop_voltage_during_burn_(V)____________________________________________ (vburn)_______________________ 3.1974E-02 OP Plasma_resistance_(ohm)_________________________________________________ (rplas)_______________________ 2.8865E-09 OP Resistive_diffusion_time_(s)____________________________________________ (res_time)____________________ 4.8432E+03 OP @@ -8348,7 +8348,7 @@ Current_drive_efficiency_model__________________________________________ (iefrf)_______________________ 10 Auxiliary_power_used_for_plasma_heating_only_(MW)_______________________ (pheat)_______________________ 5.0000E+01 Power_injected_for_current_drive_(MW)___________________________________ (pcurrentdrivemw)_____________ 1.0000E+00 - Maximum_Allowed_Bootstrap_current_fraction______________________________ (bscfmax)_____________________ 9.9000E-01 + Maximum_Allowed_Bootstrap_current_fraction______________________________ (bootstrap_current_fraction_max)_____________________ 9.9000E-01 Fusion_gain_factor_Q____________________________________________________ (bigq)________________________ 3.9003E+01 OP Auxiliary_current_drive_(A)_____________________________________________ (auxiliary_cd)________________ 4.5950E+04 OP Current_drive_efficiency_(A/W)__________________________________________ (effcd)_______________________ 4.5950E-02 OP @@ -8357,10 +8357,10 @@ ECRH_plasma_heating_efficiency__________________________________________ (gamma_ecrh)__________________ 3.0000E-01 Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 - Pfirsch-Schlueter_fraction______________________________________________ (psipf)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (faccd)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (facoh)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+faccd+facoh)_________ 1.0000E+00 + Pfirsch-Schlueter_fraction______________________________________________ (ps_current_fraction)_______________________ 0.0000E+00 + Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV Electron_cyclotron_injected_power_(MW)__________________________________ (echpwr)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (pinjalw)_____________________ 5.1000E+01 @@ -9144,7 +9144,7 @@ ucme = 3.0d8 * Unit cost of maintenance equipment ($/w**0;3) *-------------Current Drive Variables--------------* -bscfmax = 0.99 * Maximum fraction of plasma current from bootstrap; +bootstrap_current_fraction_max = 0.99 * Maximum fraction of plasma current from bootstrap; iefrf = 10 * Switch for current drive efficiency model; gamma_ecrh = 0.30 * ECRH gamma_CD (user input) etaech = 0.4 * ECRH wall-plug efficiency @@ -9238,9 +9238,9 @@ fkzohm = 1.0245 * Zohm elongation scaling adjustment factor (ishape=2; 3) fvsbrnni = 0.4434 * Fraction of the plasma current produced by gamma = 0.3 * Ejima coefficient for resistive startup v-s formula hfact = 1.1 * H factor on energy confinement times (iteration variable 10) -ibss = 4 * Switch for bootstrap current scaling; +i_bootstrap_current = 4 * Switch for bootstrap current scaling; iculbl = 1 * Switch for beta limit scaling (constraint equation 24); -icurr = 4 * Switch for plasma current scaling to use; +i_plasma_current = 4 * Switch for plasma current scaling to use; idensl = 7 * Switch for density limit to enforce (constraint equation 5); ifalphap = 1 * Switch for fast alpha pressure calculation; ifispact = 0 * Switch for neutronics calculations; diff --git a/examples/data/scan_example_file_IN.DAT b/examples/data/scan_example_file_IN.DAT index 0cae5c4a77..9d42b64771 100644 --- a/examples/data/scan_example_file_IN.DAT +++ b/examples/data/scan_example_file_IN.DAT @@ -372,13 +372,13 @@ fkzohm = 1.02 gamma = 0.3 * Switch for bootstrap current scaling -ibss = 4 +i_bootstrap_current = 4 * Switch for beta limit scaling iculbl = 1 * Switch for plasma current scaling -icurr = 4 +i_plasma_current = 4 * Switch for density limit to enforce idensl = 7 @@ -453,7 +453,7 @@ tramp = 500.0 *---------------* * Maximum fraction of plasma current from bootstrap -bscfmax = 0.95 +bootstrap_current_fraction_max = 0.95 * Switch for current drive efficiency model iefrf = 10 diff --git a/mkdocs.yml b/mkdocs.yml index 9a9290374e..5658112ef1 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -52,9 +52,13 @@ nav: - Beta Limit: physics-models/plasma_beta.md - Fast Alpha: physics-models/plasma_alpha.md - Density Limit: physics-models/plasma_density.md - - Composition & Impurities: physics-models/plasma_radiation_impurities.md - - Radiation: physics-models/plasma_radiation_impurities.md - - Plasma Current: physics-models/plasma_current.md + - Composition, Impurities & Radiation: physics-models/plasma_radiation_impurities.md + - Plasma Current: + - Overview: physics-models/plasma_current/plasma_current.md + - Bootstrap Current: physics-models/plasma_current/bootstrap_current.md + - Diamagnetic Current: physics-models/plasma_current/diamagnetic_current.md + - Pfirsch-Schlüter Current: physics-models/plasma_current/pfirsch_schlüter_current_drive.md + - Inductive Current: physics-models/plasma_current/inductive_plasma_current.md - Confinement time: physics-models/plasma_confinement.md - Plasma Core Power Balance: physics-models/plasma_power_balance.md - Pulsed Plant Operation: physics-models/pulsed-plant.md diff --git a/process/current_drive.py b/process/current_drive.py index 64e5393463..87d6e12b67 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -41,7 +41,7 @@ def cudriv(self, output: bool): pinjemwfix = 0.0 pinjimwfix = 0.0 auxiliary_cdfix = 0.0 - faccdfix = 0.0 + aux_current_fraction_fix = 0.0 gamcdfix = 0.0e0 # To stop issues with input file we force @@ -284,7 +284,9 @@ def cudriv(self, output: bool): ) * 1.0e6 ) - faccdfix = auxiliary_cdfix / physics_variables.plascur + aux_current_fraction_fix = ( + auxiliary_cdfix / physics_variables.plasma_current + ) elif current_drive_variables.iefrffix in [3, 7, 10, 12, 13]: # Injected power pinjemwfix = current_drive_variables.pinjfixmw @@ -306,7 +308,9 @@ def cudriv(self, output: bool): ) * 1.0e6 ) - faccdfix = auxiliary_cdfix / physics_variables.plascur + aux_current_fraction_fix = ( + auxiliary_cdfix / physics_variables.plasma_current + ) elif current_drive_variables.iefrffix in [5, 8]: # Account for first orbit losses # (power due to particles that are ionised but not thermalised) [MW]: @@ -353,7 +357,9 @@ def cudriv(self, output: bool): ) * 1.0e6 ) - faccdfix = auxiliary_cdfix / physics_variables.plascur + aux_current_fraction_fix = ( + auxiliary_cdfix / physics_variables.plasma_current + ) # Fenstermacher Lower Hybrid model if current_drive_variables.iefrf == 1: @@ -553,15 +559,21 @@ def cudriv(self, output: bool): ) # Compute current drive wall plug and injected powers (MW) and efficiencies - auxiliary_cd = physics_variables.faccd * physics_variables.plascur + auxiliary_cd = ( + physics_variables.aux_current_fraction + * physics_variables.plasma_current + ) # LHCD or ICCD if current_drive_variables.iefrf in [1, 2, 4, 6]: # Injected power current_drive_variables.plhybd = ( 1.0e-6 - * (physics_variables.faccd - faccdfix) - * physics_variables.plascur + * ( + physics_variables.aux_current_fraction + - aux_current_fraction_fix + ) + * physics_variables.plasma_current / effrfss + current_drive_variables.pheat ) @@ -585,8 +597,11 @@ def cudriv(self, output: bool): # Injected power (set to close to close the Steady-state current equilibrium) current_drive_variables.echpwr = ( 1.0e-6 - * (physics_variables.faccd - faccdfix) - * physics_variables.plascur + * ( + physics_variables.aux_current_fraction + - aux_current_fraction_fix + ) + * physics_variables.plasma_current / effrfss + current_drive_variables.pheat ) @@ -604,8 +619,11 @@ def cudriv(self, output: bool): # MDK. See Gitlab issue #248, and scanned note. power1 = ( 1.0e-6 - * (physics_variables.faccd - faccdfix) - * physics_variables.plascur + * ( + physics_variables.aux_current_fraction + - aux_current_fraction_fix + ) + * physics_variables.plasma_current / effnbss + current_drive_variables.pheat ) @@ -766,7 +784,7 @@ def cudriv(self, output: bool): po.oblnkl(self.outfile) - if abs(physics_variables.facoh) > 1.0e-8: + if abs(physics_variables.inductive_current_fraction) > 1.0e-8: po.ocmmnt(self.outfile, "Current is driven by both inductive") po.ocmmnt(self.outfile, "and non-inductive means.") @@ -793,8 +811,8 @@ def cudriv(self, output: bool): po.ovarre( self.outfile, "Maximum Allowed Bootstrap current fraction", - "(bscfmax)", - current_drive_variables.bscfmax, + "(bootstrap_current_fraction_max)", + current_drive_variables.bootstrap_current_fraction_max, ) if current_drive_variables.iefrffix != 0: po.ovarre( @@ -909,52 +927,52 @@ def cudriv(self, output: bool): po.ovarrf( self.outfile, "Bootstrap fraction", - "(bootipf)", - current_drive_variables.bootipf, + "(bootstrap_current_fraction)", + current_drive_variables.bootstrap_current_fraction, "OP ", ) po.ovarrf( self.outfile, "Diamagnetic fraction", - "(diaipf)", - current_drive_variables.diaipf, + "(diamagnetic_current_fraction)", + current_drive_variables.diamagnetic_current_fraction, "OP ", ) po.ovarrf( self.outfile, "Pfirsch-Schlueter fraction", - "(psipf)", - current_drive_variables.psipf, + "(ps_current_fraction)", + current_drive_variables.ps_current_fraction, "OP ", ) po.ovarrf( self.outfile, "Auxiliary current drive fraction", - "(faccd)", - physics_variables.faccd, + "(aux_current_fraction)", + physics_variables.aux_current_fraction, "OP ", ) po.ovarrf( self.outfile, "Inductive fraction", - "(facoh)", - physics_variables.facoh, + "(inductive_current_fraction)", + physics_variables.inductive_current_fraction, "OP ", ) # Add total error check. po.ovarrf( self.outfile, "Total", - "(plasipf+faccd+facoh)", - current_drive_variables.plasipf - + physics_variables.faccd - + physics_variables.facoh, + "(plasma_current_internal_fraction+aux_current_fraction+inductive_current_fraction)", + current_drive_variables.plasma_current_internal_fraction + + physics_variables.aux_current_fraction + + physics_variables.inductive_current_fraction, ) if ( abs( - current_drive_variables.plasipf - + physics_variables.faccd - + physics_variables.facoh + current_drive_variables.plasma_current_internal_fraction + + physics_variables.aux_current_fraction + + physics_variables.inductive_current_fraction - 1.0e0 ) > 1.0e-8 @@ -970,7 +988,10 @@ def cudriv(self, output: bool): ) if ( - abs(current_drive_variables.bootipf - current_drive_variables.bscfmax) + abs( + current_drive_variables.bootstrap_current_fraction + - current_drive_variables.bootstrap_current_fraction_max + ) < 1.0e-8 ): po.ocmmnt(self.outfile, "Warning : bootstrap current fraction is at") diff --git a/process/io/mfile_comparison.py b/process/io/mfile_comparison.py index fc77b479fc..91f53915ad 100755 --- a/process/io/mfile_comparison.py +++ b/process/io/mfile_comparison.py @@ -57,7 +57,7 @@ "blnkith", "blnkoth", "powfmw", - "plascur/1d6", + "plasma_current_MA", "bt", "q95", "betap", @@ -84,9 +84,9 @@ "pnucshld", "pdivt", "pheat", - "bootipf", - "faccd", - "facoh", + "bootstrap_current_fraction", + "aux_current_fraction", + "inductive_current_fraction", "gamnb", "enbeam", "powerht", @@ -110,7 +110,7 @@ "vol", "n_tf", "powfmw", - "plascur/1d6", + "plasma_current_MA", "bt", "q95", "beta", @@ -173,9 +173,9 @@ "pnetelmw", "pinjmw", "pheat", - "bootipf", - "faccd", - "facoh", + "bootstrap_current_fraction", + "aux_current_fraction", + "inductive_current_fraction", "gamnb", "enbeam", "powerht", @@ -216,7 +216,7 @@ "triang", "triang95", "powfmw", - "plascur/1d6", + "plasma_current_MA", "bt", "q95", "beta", @@ -229,7 +229,7 @@ "ralpne", "pinnerzoneradmw", "pradmw", - "bootipf", + "bootstrap_current_fraction", "pdivmax/rmajor", "fimp(14", "etath", diff --git a/process/io/obsolete_vars.py b/process/io/obsolete_vars.py index 80df8b359e..3a5a1f5a2b 100644 --- a/process/io/obsolete_vars.py +++ b/process/io/obsolete_vars.py @@ -113,6 +113,12 @@ "theat": "t_fusion_ramp", "ieped": None, "eped_sf": None, + "icurr": "i_plasma_current", + "idia": "i_diamagnetic_current", + "ibss": "i_bootstrap_current", + "ips": "i_pfirsch_schluter_current", + "bootipf": "bootstrap_current_fraction", + "bscfmax": "bootstrap_current_fraction_max", "vgap2": "vgap_vv_thermalshield", "vgap": "vgap_xpoint_divertor", } diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 836a7fba95..40dad0a559 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2407,7 +2407,7 @@ def plot_physics_info(axis, mfile_data, scan): data = [ ("powfmw", "Fusion power", "MW"), ("bigq", "$Q_{p}$", ""), - ("plascur/1d6", "$I_p$", "MA"), + ("plasma_current_MA", "$I_p$", "MA"), ("bt", "Vacuum $B_T$ at $R_0$", "T"), ("q95", r"$q_{\mathrm{95}}$", ""), ("normalised_thermal_beta", r"$\beta_N$, thermal", "% m T MA$^{-1}$"), @@ -2749,9 +2749,9 @@ def plot_current_drive_info(axis, mfile_data, scan): data = [ (pinjie, "Steady state auxiliary power", "MW"), ("pheat", "Power for heating only", "MW"), - ("bootipf", "Bootstrap fraction", ""), - ("faccd", "Auxiliary fraction", ""), - ("facoh", "Inductive fraction", ""), + ("bootstrap_current_fraction", "Bootstrap fraction", ""), + ("aux_current_fraction", "Auxiliary fraction", ""), + ("inductive_current_fraction", "Inductive fraction", ""), ("powerht", "Plasma heating used for H factor", "MW"), ( "effcd", @@ -2779,9 +2779,9 @@ def plot_current_drive_info(axis, mfile_data, scan): data = [ (pinjie, "Steady state auxiliary power", "MW"), ("pheat", "Power for heating only", "MW"), - ("bootipf", "Bootstrap fraction", ""), - ("faccd", "Auxiliary fraction", ""), - ("facoh", "Inductive fraction", ""), + ("bootstrap_current_fraction", "Bootstrap fraction", ""), + ("aux_current_fraction", "Auxiliary fraction", ""), + ("inductive_current_fraction", "Inductive fraction", ""), ("gamnb", "NB gamma", "$10^{20}$ A W$^{-1}$ m$^{-2}$"), ("enbeam", "NB energy", "keV"), ("powerht", "Plasma heating used for H factor", "MW"), @@ -2805,9 +2805,9 @@ def plot_current_drive_info(axis, mfile_data, scan): data = [ (pinjie, "Steady state auxiliary power", "MW"), ("pheat", "Power for heating only", "MW"), - ("bootipf", "Bootstrap fraction", ""), - ("faccd", "Auxiliary fraction", ""), - ("facoh", "Inductive fraction", ""), + ("bootstrap_current_fraction", "Bootstrap fraction", ""), + ("aux_current_fraction", "Auxiliary fraction", ""), + ("inductive_current_fraction", "Inductive fraction", ""), ("powerht", "Plasma heating used for H factor", "MW"), ( "gamcd", @@ -2834,9 +2834,9 @@ def plot_current_drive_info(axis, mfile_data, scan): data = [ (pinjie, "Steady state auxiliary power", "MW"), ("pheat", "Power for heating only", "MW"), - ("bootipf", "Bootstrap fraction", ""), - ("faccd", "Auxiliary fraction", ""), - ("facoh", "Inductive fraction", ""), + ("bootstrap_current_fraction", "Bootstrap fraction", ""), + ("aux_current_fraction", "Auxiliary fraction", ""), + ("inductive_current_fraction", "Inductive fraction", ""), ("powerht", "Plasma heating used for H factor", "MW"), ( "gamcd", @@ -2863,9 +2863,9 @@ def plot_current_drive_info(axis, mfile_data, scan): data = [ (pinjie, "Steady state auxiliary power", "MW"), ("pheat", "Power for heating only", "MW"), - ("bootipf", "Bootstrap fraction", ""), - ("faccd", "Auxiliary fraction", ""), - ("facoh", "Inductive fraction", ""), + ("bootstrap_current_fraction", "Bootstrap fraction", ""), + ("aux_current_fraction", "Auxiliary fraction", ""), + ("inductive_current_fraction", "Inductive fraction", ""), ("powerht", "Plasma heating used for H factor", "MW"), ( "gamcd", diff --git a/process/io/plot_radial_build.py b/process/io/plot_radial_build.py index 69ff228741..1b75e0e628 100644 --- a/process/io/plot_radial_build.py +++ b/process/io/plot_radial_build.py @@ -160,7 +160,7 @@ def main(args=None): "te", "boundu(15)", "dnbeta", - "bscfmax", + "bootstrap_current_fraction_max", "boundu(10)", "fiooic", "fjprot", diff --git a/process/io/plot_scans.py b/process/io/plot_scans.py index 2176755244..67e57d1f2e 100644 --- a/process/io/plot_scans.py +++ b/process/io/plot_scans.py @@ -218,7 +218,7 @@ def main(args=None): nsweep_dict[9] = "te" nsweep_dict[10] = "boundu(15)" nsweep_dict[11] = "dnbeta" - nsweep_dict[12] = "bscfmax" + nsweep_dict[12] = "bootstrap_current_fraction_max" nsweep_dict[13] = "boundu(10)" nsweep_dict[14] = "fiooic" nsweep_dict[15] = "fjprot" @@ -632,12 +632,12 @@ def main(args=None): plt.tight_layout() # Output file naming - if output_name == "plascur/1d6": + if output_name == "plasma_current_MA": plt.savefig( - f"{args.outputdir}/scan_{scan_var_name}_vs_plascur" + f"{args.outputdir}/scan_{scan_var_name}_vs_plasma_current" + f"_vs_{output_name2}" if output_names2 != [] - else f"{args.outputdir}/scan_{scan_var_name}_vs_plascur" + else f"{args.outputdir}/scan_{scan_var_name}_vs_plasma_current" + f".{save_format}" ) elif output_name == "pdivtbt/qar": diff --git a/process/io/variable_metadata.py b/process/io/variable_metadata.py index 744fba437f..116ffbf93a 100644 --- a/process/io/variable_metadata.py +++ b/process/io/variable_metadata.py @@ -161,7 +161,7 @@ class VariableMetadata: "f_tf_steel": VariableMetadata( latex=r"f_\mathrm{steel}^\mathrm{TF}", description="TF steel fraction", units="" ), - "plascur/1d6": VariableMetadata( + "plasma_current_MA": VariableMetadata( latex=r"$I_{\mathrm{p}}$[$MA$]", description="Plasma current", units="MA" ), "n_cycle": VariableMetadata( @@ -230,13 +230,13 @@ class VariableMetadata: latex=r"$\eta_{\mathrm{CD}}$[$A/W$]", description="CD efficiency", units="A/W" ), "bigq": VariableMetadata(latex=r"$Q$", description="Plasma Q value", units=""), - "faccd": VariableMetadata( + "aux_current_fraction": VariableMetadata( latex=r"$f_{\mathrm{CD}}$", description="CD factor", units="" ), - "facoh": VariableMetadata( + "inductive_current_fraction": VariableMetadata( latex=r"$f_{\mathrm{CD,ind}}$", description="Inductive CD factor", units="" ), - "bootipf": VariableMetadata( + "bootstrap_current_fraction": VariableMetadata( latex=r"$f_{\mathrm{BS}}$", description="Bootstrap current fraction", units="" ), "pdivt": VariableMetadata( diff --git a/process/pfcoil.py b/process/pfcoil.py index 4cafb95b39..c52299e1c7 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -296,11 +296,11 @@ def pfcoil(self): elif pfv.ipfloc[i] == 2: # PF coil is on top of the TF coil - pf.ccls[i] = 0.3e0 * pv.aspect**1.6e0 * pv.plascur + pf.ccls[i] = 0.3e0 * pv.aspect**1.6e0 * pv.plasma_current elif pfv.ipfloc[i] == 3: # PF coil is radially outside the TF coil - pf.ccls[i] = -0.4e0 * pv.plascur + pf.ccls[i] = -0.4e0 * pv.plasma_current else: eh.idiags[0] = i @@ -310,7 +310,7 @@ def pfcoil(self): # Vertical field (T) pv.bvert = ( -1.0e-7 - * pv.plascur + * pv.plasma_current / pv.rmajor * ( math.log(8.0e0 * pv.aspect) @@ -343,7 +343,7 @@ def pfcoil(self): # This is a fixed current for this calculation -- RK 07/12 pf.ccls[i] = ( - pv.plascur + pv.plasma_current * 2.0e0 * (1.0e0 - (pv.kappa * pv.rminor) / abs(pf.zcls[i, 0])) ) @@ -391,7 +391,7 @@ def pfcoil(self): bzin[0] = ( -1.0e-7 - * pv.plascur + * pv.plasma_current / pv.rmajor * ( math.log(8.0e0 * pv.aspect) @@ -767,9 +767,9 @@ def pfcoil(self): # Plasma wave form pfv.cpt[pfv.ncirt - 1, 0] = 0.0e0 pfv.cpt[pfv.ncirt - 1, 1] = 0.0e0 - pfv.cpt[pfv.ncirt - 1, 2] = pv.plascur - pfv.cpt[pfv.ncirt - 1, 3] = pv.plascur - pfv.cpt[pfv.ncirt - 1, 4] = pv.plascur + pfv.cpt[pfv.ncirt - 1, 2] = pv.plasma_current + pfv.cpt[pfv.ncirt - 1, 3] = pv.plasma_current + pfv.cpt[pfv.ncirt - 1, 4] = pv.plasma_current pfv.cpt[pfv.ncirt - 1, 5] = 0.0e0 def efc( @@ -1112,7 +1112,7 @@ def ohcalc(self): # Calculation of CS fatigue # this is only valid for pulsed reactor design - if pv.facoh > 0.0e-4: + if pv.inductive_current_fraction > 0.0e-4: csfv.n_cycle, csfv.t_crack_radial = self.cs_fatigue.ncycle( pf.sig_hoop, csfv.residual_sig_hoop, @@ -1348,7 +1348,7 @@ def peakb(self, i, ii, it): kk = kk + 1 pf.rfxf[kk - 1] = pv.rmajor pf.zfxf[kk - 1] = 0.0e0 - pf.cfxf[kk - 1] = pv.plascur + pf.cfxf[kk - 1] = pv.plasma_current # Calculate the field at the inner and outer edges # of the coil of interest @@ -2126,7 +2126,7 @@ def outpf(self): tfv.tmargmin_cs, ) # only output CS fatigue model for pulsed reactor design - if pv.facoh > 0.0e-4: + if pv.inductive_current_fraction > 0.0e-4: op.ovarre( self.outfile, "Residual hoop stress in CS Steel (Pa)", diff --git a/process/physics.py b/process/physics.py index 0aa4345617..cff586ad89 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1,10 +1,13 @@ +import math +from typing import Tuple import numpy as np import numba as nb -import scipy.integrate as integrate from scipy.optimize import root_scalar -import math -import process.physics_functions as physics_funcs from process.utilities.f2py_string_patch import f2py_compatible_to_string + +import scipy.integrate as integrate + +import process.physics_functions as physics_funcs from process.fortran import ( constraint_variables, reinke_variables, @@ -27,159 +30,6 @@ ) -@nb.jit(nopython=True, cache=True) -def coulg(j, tempe, ne): - """Coulomb logarithm - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the Coulomb logarithm, valid - for e-e collisions (T_e > 0.01 keV), and for - e-i collisions (T_e > 0.01*Zeff^2) (Alexander, 9/5/1994). -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - C. A. Ordonez and M. I. Molina, Phys. Plasmas 1 (1994) 2515 - Rev. Mod. Phys., V.48, Part 1 (1976) 275 - """ - return 15.9 - 0.5 * np.log(ne[j - 1]) + np.log(tempe[j - 1]) - - -@nb.jit(nopython=True, cache=True) -def nuee(j, tempe, ne): - """Frequency of electron-electron collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the frequency of electron-electron - collisions (Hz): NUEE = 4*SQRT(pi)/3*Ne*e**4*lambd/ - SQRT(Me)/Te**1.5 -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - Yushmanov, 25th April 1987 (?), - updated by Pereverzev, 9th November 1994 (?) - """ - return ( - 670.0 * coulg(j, tempe, ne) * ne[j - 1] / (tempe[j - 1] * np.sqrt(tempe[j - 1])) - ) - - -@nb.jit(nopython=True, cache=True) -def nui(j, zmain, ni, tempi, amain): - """Full frequency of ion collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the full frequency of ion - collisions (Hz). -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - None - """ - # Coulomb logarithm = 15 is used - return ( - zmain[j - 1] ** 4 - * ni[j - 1] - * 322.0 - / (tempi[j - 1] * np.sqrt(tempi[j - 1] * amain[j - 1])) - ) - - -@nb.jit(nopython=True, cache=True) -def _plascar_bpol(aspect, eps, kappa, delta): - # Original coding, only suitable for TARTs [STAR Code] - - c1 = kappa**2 / (1.0 + delta) + delta - c2 = kappa**2 / (1.0 - delta) - delta - - d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 - d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 - - c1_aspect = (c1 * eps - 1.0) if aspect < c1 else (1.0 - c1 * eps) - - y1 = np.sqrt(c1_aspect / (1.0 + eps)) * (1.0 + delta) / kappa - y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * (1.0 - delta) / kappa - - h2 = (1.0 + (c2 - 1.0) * eps / 2.0) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) - f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * (c2 * eps + 1.0)) - g = eps * kappa / (1.0 - eps * delta) - ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) - - h1 = (1.0 + (1.0 - c1) * eps / 2.0) / np.sqrt((1.0 + eps) * c1_aspect) - f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - - if aspect < c1: - ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) - else: - ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) - - return ff1, ff2, d1, d2 - - -@nb.jit(nopython=True, cache=True) -def bpol(icurr, ip, qbar, aspect, eps, bt, kappa, delta, perim, rmu0): - """Function to calculate poloidal field - author: J Galambos, FEDC/ORNL - author: P J Knight, CCFE, Culham Science Centre - icurr : input integer : current scaling model to use - ip : input real : plasma current (A) - qbar : input real : edge q-bar - aspect : input real : plasma aspect ratio - bt : input real : toroidal field on axis (T) - kappa : input real : plasma elongation - delta : input real : plasma triangularity - perim : input real : plasma perimeter (m) - This function calculates the poloidal field in Tesla, - using a simple calculation using Stoke's Law for conventional - tokamaks, or for TARTs, a scaling from Peng, Galambos and - Shipe (1992). - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, - Fusion Technology, 21, 1729 - """ - if icurr != 2: - return rmu0 * ip / perim - - ff1, ff2, _, _ = _plascar_bpol(aspect, eps, kappa, delta) - - return bt * (ff1 + ff2) / (2.0 * np.pi * qbar) - - -@nb.jit(nopython=True, cache=True) -def plasc(qbar, aspect, eps, rminor, bt, kappa, delta): - """Function to calculate plasma current (Peng scaling) - author: J Galambos, FEDC/ORNL - author: P J Knight, CCFE, Culham Science Centre - aspect : input real : plasma aspect ratio - bt : input real : toroidal field on axis (T) - delta : input real : plasma triangularity - kappa : input real : plasma elongation - qbar : input real : edge q-bar - rminor : input real : plasma minor radius (m) - This function calculates the plasma current in MA, - using a scaling from Peng, Galambos and Shipe (1992). - It is primarily used for Tight Aspect Ratio Tokamaks and is - selected via icurr=2. - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, - Fusion Technology, 21, 1729 - """ - - ff1, ff2, d1, d2 = _plascar_bpol(aspect, eps, kappa, delta) - - e1 = 2.0 * kappa / (d1 * (1.0 + delta)) - e2 = 2.0 * kappa / (d2 * (1.0 - delta)) - - return ( - rminor - * bt - / qbar - * 5.0 - * kappa - / (2.0 * np.pi**2) - * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) - * (ff1 + ff2) - ) - - @nb.jit(nopython=True, cache=True) def rether(alphan, alphat, dene, dlamie, te, ti, zeffai): """Routine to find the equilibration power between the @@ -209,12 +59,12 @@ def rether(alphan, alphat, dene, dlamie, te, ti, zeffai): def vscalc( csawth, eps, - facoh, + inductive_current_fraction, gamma, kappa, rmajor, rplas, - plascur, + plasma_current, t_fusion_ramp, tburn, rli, @@ -224,10 +74,10 @@ def vscalc( author: P J Knight, CCFE, Culham Science Centre csawth : input real : coefficient for sawteeth effects eps : input real : inverse aspect ratio - facoh : input real : fraction of plasma current produced inductively + inductive_current_fraction : input real : fraction of plasma current produced inductively gamma : input real : Ejima coeff for resistive start-up V-s component kappa : input real : plasma elongation - plascur: input real : plasma current (A) + plasma_current: input real : plasma current (A) rli : input real : plasma normalised inductivity rmajor : input real : plasma major radius (m) rplas : input real : plasma resistance (ohm) @@ -245,12 +95,12 @@ def vscalc( # Internal inductance rlpint = rmu0 * rmajor * rli / 2.0 - phiint = rlpint * plascur + phiint = rlpint * plasma_current # Start-up resistive component # Uses ITER formula without the 10 V-s add-on - vsres = gamma * rmu0 * plascur * rmajor + vsres = gamma * rmu0 * plasma_current * rmajor # Hirshman, Neilson: Physics of Fluids, 29 (1986) p790 # fit for external inductance @@ -267,14 +117,14 @@ def vscalc( # Inductive V-s component - vsind = rlp * plascur + vsind = rlp * plasma_current vsstt = vsres + vsind # Loop voltage during flat-top # Include enhancement factor in flattop V-s requirement # to account for MHD sawtooth effects. - vburn = plascur * rplas * facoh * csawth + vburn = plasma_current * rplas * inductive_current_fraction * csawth # N.B. tburn on first iteration will not be correct # if the pulsed reactor option is used, but the value @@ -287,12 +137,12 @@ def vscalc( @nb.jit(nopython=True, cache=True) -def culblm(bt, dnbeta, plascur, rminor): +def culblm(bt, dnbeta, plasma_current, rminor): """Beta scaling limit author: P J Knight, CCFE, Culham Science Centre bt : input real : toroidal B-field on plasma axis (T) dnbeta : input real : Troyon-like g coefficient - plascur : input real : plasma current (A) + plasma_current : input real : plasma current (A) rminor : input real : plasma minor axis (m) betalim : output real : beta limit as defined below This subroutine calculates the beta limit, using @@ -309,30 +159,300 @@ def culblm(bt, dnbeta, plascur, rminor): AEA FUS 172: Physics Assessment for the European Reactor Study """ - return 0.01 * dnbeta * (plascur / 1.0e6) / (rminor * bt) + return 0.01 * dnbeta * (plasma_current / 1.0e6) / (rminor * bt) + + +# ----------------------------------------------------- +# Plasma Current & Poloidal Field Calculations +# ----------------------------------------------------- @nb.jit(nopython=True, cache=True) -def conhas(alphaj, alphap, bt, delta95, eps, kappa95, p0, rmu0): - """Routine to calculate the F coefficient used for scaling the - plasma current - author: P J Knight, CCFE, Culham Science Centre - alphaj : input real : current profile index - alphap : input real : pressure profile index - bt : input real : toroidal field on axis (T) - delta95 : input real : plasma triangularity 95% - eps : input real : inverse aspect ratio - kappa95 : input real : plasma elongation 95% - p0 : input real : central plasma pressure (Pa) - fq : output real : scaling for edge q from circular - cross-section cylindrical case - This routine calculates the F coefficient used for scaling the - plasma current, using the Connor-Hastie scaling given in - AEA FUS 172: Physics Assessment for the European Reactor Study +def _plascar_bpol( + aspect: float, eps: float, kappa: float, delta: float +) -> Tuple[float, float, float, float]: + """ + Calculate the poloidal field coefficients for determining the plasma current + and poloidal field. + + Parameters: + - aspect: float, plasma aspect ratio + - eps: float, inverse aspect ratio + - kappa: float, plasma elongation + - delta: float, plasma triangularity + + Returns: + - Tuple[float, float, float, float], coefficients ff1, ff2, d1, d2 + + This internal function calculates the poloidal field coefficients, + which is used to calculate the poloidal field and the plasma current. + + References: + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729–1738. https://doi.org/10.13182/FST92-A29971 + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + """ + # Original coding, only suitable for TARTs [STAR Code] + + c1 = (kappa**2 / (1.0 + delta)) + delta + c2 = (kappa**2 / (1.0 - delta)) - delta + + d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 + d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 + + c1_aspect = ((c1 * eps) - 1.0) if aspect < c1 else (1.0 - (c1 * eps)) + + y1 = np.sqrt(c1_aspect / (1.0 + eps)) * ((1.0 + delta) / kappa) + y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * ((1.0 - delta) / kappa) + + h2 = (1.0 + (c2 - 1.0) * (eps / 2.0)) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) + f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * ((c2 * eps) + 1.0)) + g = (eps * kappa) / (1.0 - (eps * delta)) + ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) + + h1 = (1.0 + (1.0 - c1) * (eps / 2.0)) / np.sqrt((1.0 + eps) * c1_aspect) + f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) + + if aspect < c1: + ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) + else: + ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) + + return ff1, ff2, d1, d2 + + +def calculate_poloidal_field( + i_plasma_current: int, + ip: float, + q95: float, + aspect: float, + eps: float, + bt: float, + kappa: float, + delta: float, + perim: float, + rmu0: float, +) -> float: + """ + Function to calculate poloidal field from the plasma current + + Parameters: + - i_plasma_current: int, current scaling model to use + - ip: float, plasma current (A) + - q95: float, 95% flux surface safety factor + - aspect: float, plasma aspect ratio + - eps: float, inverse aspect ratio + - bt: float, toroidal field on axis (T) + - kappa: float, plasma elongation + - delta: float, plasma triangularity + - perim: float, plasma perimeter (m) + - rmu0: float, vacuum permeability (H/m) + + Returns: + - float, poloidal field in Tesla + + This function calculates the poloidal field from the plasma current in Tesla, + using a simple calculation using Ampere's law for conventional + tokamaks, or for TARTs, a scaling from Peng, Galambos and + Shipe (1992). + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729–1738. https://doi.org/10.13182/FST92-A29971 + """ + # Use Ampere's law using the plasma poloidal cross-section + if i_plasma_current != 2: + return rmu0 * ip / perim + else: + # Use the relation from Peng, Galambos and Shipe (1992) [STAR code] otherwise + ff1, ff2, _, _ = _plascar_bpol(aspect, eps, kappa, delta) + + # Transform q95 to qbar + qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 + + return bt * (ff1 + ff2) / (2.0 * np.pi * qbar) + + +def calculate_current_coefficient_peng(eps: float, sf: float) -> float: + """ + Calculate the plasma current scaling coefficient for the Peng scaling from the STAR code. + + Parameters: + - eps: float, plasma inverse aspect ratio + - sf: float, shaping factor calculated in the poloidal perimeter function + + References: + - None + """ + return (1.22 - 0.68 * eps) / ((1.0 - eps * eps) ** 2) * sf**2 + + +def calculate_plasma_current_peng( + q95: float, + aspect: float, + eps: float, + rminor: float, + bt: float, + kappa: float, + delta: float, +) -> float: + """ + Function to calculate plasma current (Peng scaling from the STAR code) + + Parameters: + - q95: float, 95% flux surface safety factor + - aspect: float, plasma aspect ratio + - eps: float, inverse aspect ratio + - rminor: float, plasma minor radius (m) + - bt: float, toroidal field on axis (T) + - kappa: float, plasma elongation + - delta: float, plasma triangularity + + Returns: + - float, plasma current in MA + + This function calculates the plasma current in MA, + using a scaling from Peng, Galambos and Shipe (1992). + It is primarily used for Tight Aspect Ratio Tokamaks and is + selected via i_plasma_current=2. + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729–1738. https://doi.org/10.13182/FST92-A29971 + """ + + # Transform q95 to qbar + qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 + + ff1, ff2, d1, d2 = _plascar_bpol(aspect, eps, kappa, delta) + + e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) + e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) + + return ( + rminor + * bt + / qbar + * 5.0 + * kappa + / (2.0 * np.pi**2) + * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) + * (ff1 + ff2) + ) + + +@nb.jit(nopython=True, cache=True) +def calculate_current_coefficient_ipdg89( + eps: float, kappa95: float, triang95: float +) -> float: + """ + Calculate the fq coefficient from the IPDG89 guidlines used in the plasma current scaling. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa95: float, plasma elongation 95% + - triang95: float, plasma triangularity 95% + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient used in the IPDG89 plasma current scaling, + based on the given plasma parameters. + + References: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 0.5 + * (1.17 - 0.65 * eps) + / ((1.0 - eps * eps) ** 2) + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + ) + + +@nb.jit(nopython=True, cache=True) +def calculate_current_coefficient_todd( + eps: float, kappa95: float, triang95: float, model: int +) -> float: + """ + Calculate the fq coefficient used in the two Todd plasma current scalings. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa95: float, plasma elongation 95% + - triang95: float, plasma triangularity 95% + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient based on the given plasma parameters for the two Todd scalings. + + References: + - D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 """ + # Calculate the Todd scaling based on the model + base_scaling = ( + (1.0 + 2.0 * eps**2) + * ((1.0 + kappa95**2) / 0.5) + * ( + 1.24 + - 0.54 * kappa95 + + 0.3 * (kappa95**2 + triang95**2) + + 0.125 * triang95 + ) + ) + if model == 1: + return base_scaling + elif model == 2: + return base_scaling * (1.0 + (abs(kappa95 - 1.2)) ** 3) - # Exponent in Connor-Hastie current profile - matching total - # current gives the following trivial relation + +@nb.jit(nopython=True, cache=True) +def calculate_current_coefficient_hastie( + alphaj: float, + alphap: float, + bt: float, + delta95: float, + eps: float, + kappa95: float, + p0: float, + rmu0: float, +) -> float: + """ + Routine to calculate the f_q coefficient for the Connor-Hastie model used for scaling the plasma current. + + Parameters: + - alphaj: float, the current profile index + - alphap: float, the pressure profile index + - bt: float, the toroidal field on axis (T) + - delta95: float, the plasma triangularity 95% + - eps: float, the inverse aspect ratio + - kappa95: float, the plasma elongation 95% + - p0: float, the central plasma pressure (Pa) + - rmu0: float, the vacuum permeability (H/m) + + Returns: + - float, the F coefficient + + This routine calculates the f_q coefficient used for scaling the plasma current, + using the Connor-Hastie scaling + + Reference: + - J.W.Connor and R.J.Hastie, Culham Lab Report CLM-M106 (1985). + https://scientific-publications.ukaea.uk/wp-content/uploads/CLM-M106-1.pdf + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + # Exponent in Connor-Hastie current profile lamda = alphaj # Exponent in Connor-Hastie pressure profile @@ -379,389 +499,966 @@ def conhas(alphaj, alphap, bt, delta95, eps, kappa95, p0, rmu0): ) -def bsinteg( - y, - dene, - ten, - bt, - rminor, - rmajor, - zeff, - alphat, - alphan, - q0, - q95, - betat, - echarge, - rmu0, -): - """Integrand function for Nevins et al bootstrap current scaling - author: P J Knight, CCFE, Culham Science Centre - y : input real : abscissa of integration, = normalised minor radius - This function calculates the integrand function for the - Nevins et al bootstrap current scaling, 4/11/90. +@nb.jit(nopython=True, cache=True) +def calculate_current_coefficient_sauter( + eps: float, + kappa: float, + triang: float, +) -> float: + """ + Routine to calculate the f_q coefficient for the Sauter model used for scaling the plasma current. + + Parameters: + - eps: float, inverse aspect ratio + - kappa: float, plasma elongation at the separatrix + - triang: float, plasma triangularity at the separatrix + + Returns: + - float, the fq coefficient + + Reference: + - O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, + Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, + ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2016.04.033. + """ + w07 = 1.0 # zero squareness - can be modified later if required + + fq = ( + (4.1e6 / 5.0e6) + * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) + * (1.0 + 0.09 * triang + 0.16 * triang**2) + * (1.0 + 0.45 * triang * eps) + / (1.0 - 0.74 * eps) + * (1.0 + 0.55 * (w07 - 1.0)) + ) + + return fq + + +@nb.jit(nopython=True, cache=True) +def calculate_current_coefficient_fiesta( + eps: float, kappa: float, triang: float +) -> float: + """ + Calculate the fq coefficient used in the FIESTA plasma current scaling. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa: float, plasma elongation at the separatrix + - triang: float, plasma triangularity at the separatrix + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient based on the given plasma parameters for the FIESTA scaling. + + References: + - S.Muldrew et.al,“PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, + Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. + """ + return 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060 + + +# -------------------------------- +# Bootstrap Current Calculations +# -------------------------------- + + +def _nevins_integral( + y: float, + dene: float, + te: float, + bt: float, + rminor: float, + rmajor: float, + zeff: float, + alphat: float, + alphan: float, + q0: float, + q95: float, + betat: float, +) -> float: + """ + Integrand function for Nevins et al bootstrap current scaling. + + Parameters: + - y: float, abscissa of integration, normalized minor radius + - dene: float, volume averaged electron density (/m^3) + - te: float, volume averaged electron temperature (keV) + - bt: float, toroidal field on axis (T) + - rminor: float, plasma minor radius (m) + - rmajor: float, plasma major radius (m) + - zeff: float, plasma effective charge + - alphat: float, temperature profile index + - alphan: float, density profile index + - q0: float, normalized safety factor at the magnetic axis + - q95: float, normalized safety factor at 95% of the plasma radius + - betat: float, Toroidal plasma beta + + Returns: + - float, the integrand value + + This function calculates the integrand function for the Nevins et al bootstrap current scaling. + + Reference: See appendix of: + Keii Gi, Makoto Nakamura, Kenji Tobita, Yasushi Ono, + Bootstrap current fraction scaling for a tokamak reactor design study, + Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, + ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2014.07.009. + + Nevins, W. M. "Summary report: ITER specialists’ meeting on heating and current drive." + ITER-TN-PH-8-4, June 1988. 1988. + """ - # Constants for fit to q-profile - c1 = 1.0 - c2 = 1.0 - c3 = 1.0 # Compute average electron beta - betae = dene * ten * 1.0e3 * echarge / (bt**2 / (2.0 * rmu0)) + betae = dene * te * 1.0e3 * constants.echarge / (bt**2 / (2.0 * constants.rmu0)) nabla = rminor * np.sqrt(y) / rmajor x = (1.46 * np.sqrt(nabla) + 2.4 * nabla) / (1.0 - nabla) ** 1.5 z = zeff d = ( 1.414 * z - + z * z - + x * (0.754 + 2.657 * z + 2.0 * z * z) - + x * x * (0.348 + 1.243 * z + z * z) + + z**2 + + x * (0.754 + 2.657 * z + (2.0 * z**2)) + + (x**2 * (0.348 + 1.243 * z + z**2)) ) - al2 = -x * (0.884 + 2.074 * z) / d + a1 = (alphan + alphat) * (1.0 - y) ** (alphan + alphat - 1.0) a2 = alphat * (1.0 - y) ** (alphan + alphat - 1.0) + al1 = (x / d) * (0.754 + 2.21 * z + z**2 + x * (0.348 + 1.243 * z + z**2)) + al2 = -x * ((0.884 + 2.074 * z) / d) alphai = -1.172 / (1.0 + 0.462 * x) - a1 = (alphan + alphat) * (1.0 - y) ** (alphan + alphat - 1.0) - al1 = x * (0.754 + 2.21 * z + z * z + x * (0.348 + 1.243 * z + z * z)) / d # q-profile - q = q0 + (q95 - q0) * (c1 * y + c2 * y * y + c3 * y**3) / (c1 + c2 + c3) + q = q0 + (q95 - q0) * ((y + y**2 + y**3) / (3.0)) pratio = (betat - betae) / betae - return (q / q95) * (al1 * (a1 + pratio * (a1 + alphai * a2)) + al2 * a2) + return (q / q95) * (al1 * (a1 + (pratio * (a1 + alphai * a2))) + al2 * a2) -def diamagnetic_fraction_hender(beta): - """author: S.I. Muldrew, CCFE, Culham Science Centre - Diamagnetic contribution at tight aspect ratio. - Tim Hender fit +# ----------------------------------------------------- +# Diamagnetic and Pfirsch-Schlüter Current Calculations +# ----------------------------------------------------- + + +@nb.jit(nopython=True, cache=True) +def diamagnetic_fraction_hender(beta: float) -> float: + """ + Calculate the diamagnetic fraction based on the Hender fit. + + Parameters: + - beta: float, the plasma beta value + + Returns: + - float, the diamagnetic fraction """ return beta / 2.8 -def diamagnetic_fraction_scene(beta, q95, q0): - """author: S.I. Muldrew, CCFE, Culham Science Centre - Diamagnetic fraction based on SCENE fit by Tim Hender - See Issue #992 +@nb.jit(nopython=True, cache=True) +def diamagnetic_fraction_scene(beta: float, q95: float, q0: float) -> float: + """ + Calculate the diamagnetic fraction based on the SCENE fit by Tim Hender. + + Parameters: + - beta: float, the plasma beta value + - q95: float, the normalized safety factor at 95% of the plasma radius + - q0: float, the normalized safety factor at the magnetic axis + + Returns: + - float, the diamagnetic fraction + """ + return beta * (0.1 * q95 / q0 + 0.44) * 0.414 + + +@nb.jit(nopython=True, cache=True) +def ps_fraction_scene(beta: float) -> float: """ - return beta * (0.1 * q95 / q0 + 0.44) * 4.14e-1 + Calculate the Pfirsch-Schlüter fraction based on the SCENE fit by Tim Hender 2019. + Parameters: + - beta: float, the plasma beta value -def ps_fraction_scene(beta): - """author: S.I. Muldrew, CCFE, Culham Science Centre - Pfirsch-Schlüter fraction based on SCENE fit by Tim Hender - See Issue #992 + Returns: + - float, the Pfirsch-Schlüter current fraction """ return -9e-2 * beta +# -------------------------------------------------- +# Sauter Bootstrap Current Scaling Functions +# -------------------------------------------------- + + @nb.jit(nopython=True, cache=True) -def nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni): - """Relative frequency of ion collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the relative frequency of ion - collisions: NU* = Nui*q*Rt/eps**1.5/Vti - The full ion collision frequency NUI is used. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - Yushmanov, 30th April 1987 (?) +def _coulomb_logarithm_sauter( + radial_elements: int, tempe: np.ndarray, ne: np.ndarray +) -> np.ndarray: + """ + Calculate the Coulomb logarithm used in the arrays for the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, the radial element indexes in the range 1 to nr + - tempe: np.ndarray, the electron temperature array + - ne: np.ndarray, the electron density array + + Returns: + - np.ndarray, the Coulomb logarithm at each array point + + This function calculates the Coulomb logarithm, which is valid for e-e collisions (T_e > 0.01 keV) + and for e-i collisions (T_e > 0.01*Zeff^2) (Alexander, 9/5/1994). + + Reference: + - C. A. Ordonez, M. I. Molina; + Evaluation of the Coulomb logarithm using cutoff and screened Coulomb interaction potentials. + Phys. Plasmas 1 August 1994; 1 (8): 2515–2518. https://doi.org/10.1063/1.870578 + - Y. R. Shen, “Recent advances in nonlinear optics,” Reviews of Modern Physics, vol. 48, no. 1, + pp. 1–32, Jan. 1976, doi: https://doi.org/10.1103/revmodphys.48.1. + """ + return ( + 15.9 + - 0.5 * np.log(ne[radial_elements - 1]) + + np.log(tempe[radial_elements - 1]) + ) + + +@nb.jit(nopython=True, cache=True) +def _electron_collisions_sauter( + radial_elements: np.ndarray, tempe: np.ndarray, ne: np.ndarray +) -> np.ndarray: + """ + Calculate the frequency of electron-electron collisions used in the arrays for the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, the radial element index in the range 1 to nr + - tempe: np.ndarray, the electron temperature array + - ne: np.ndarray, the electron density array + + Returns: + - float, the frequency of electron-electron collisions (Hz) + + This function calculates the frequency of electron-electron collisions + + Reference: + - Yushmanov, 25th April 1987 (?), updated by Pereverzev, 9th November 1994 (?) + """ + return ( + 670.0 + * _coulomb_logarithm_sauter(radial_elements, tempe, ne) + * ne[radial_elements - 1] + / (tempe[radial_elements - 1] * np.sqrt(tempe[radial_elements - 1])) + ) + + +@nb.jit(nopython=True, cache=True) +def _electron_collisionality_sauter( + radial_elements: np.ndarray, + rmajor: float, + zeff: np.ndarray, + inverse_q: np.ndarray, + sqeps: np.ndarray, + tempe: np.ndarray, + ne: np.ndarray, +) -> np.ndarray: + """ + Calculate the electron collisionality used in the arrays for the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, the radial element index in the range 1 to nr + - rmajor: float, the plasma major radius (m) + - zeff: np.ndarray, the effective charge array + - inverse_q: np.ndarray, inverse safety factor profile + - sqeps: np.ndarray, the square root of the inverse aspect ratio array + - tempe: np.ndarray, the electron temperature array + - ne: np.ndarray, the electron density array + + Returns: + - np.ndarray, the relative frequency of electron collisions + + Reference: + - Yushmanov, 30th April 1987 (?) + """ + return ( + _electron_collisions_sauter(radial_elements, tempe, ne) + * 1.4 + * zeff[radial_elements - 1] + * rmajor + / np.abs( + inverse_q[radial_elements - 1] + * (sqeps[radial_elements - 1] ** 3) + * np.sqrt(tempe[radial_elements - 1]) + * 1.875e7 + ) + ) + + +@nb.jit(nopython=True, cache=True) +def _ion_collisions_sauter( + radial_elements: np.ndarray, + zeff: np.ndarray, + ni: np.ndarray, + tempi: np.ndarray, + amain: np.ndarray, +) -> np.ndarray: + """ + Calculate the full frequency of ion collisions used in the arrays for the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, the radial element indexes in the range 1 to nr + - zeff: np.ndarray, the effective charge array + - ni: np.ndarray, the ion density array + - tempi: np.ndarray, the ion temperature array + - amain: np.ndarray, the atomic mass of the main ion species array + + Returns: + - np.ndarray, the full frequency of ion collisions (Hz) + + This function calculates the full frequency of ion collisions using the Coulomb logarithm of 15. + + Reference: + - None + """ + return ( + zeff[radial_elements - 1] ** 4 + * ni[radial_elements - 1] + * 322.0 + / ( + tempi[radial_elements - 1] + * np.sqrt(tempi[radial_elements - 1] * amain[radial_elements - 1]) + ) + ) + + +@nb.jit(nopython=True, cache=True) +def _ion_collisionality_sauter( + radial_elements: np.ndarray, + rmajor: float, + inverse_q: np.ndarray, + sqeps: np.ndarray, + tempi: np.ndarray, + amain: np.ndarray, + zeff: np.ndarray, + ni: np.ndarray, +) -> float: + """ + Calculate the ion collisionality to be used in the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, the radial element indexes in the range 1 to nr + - rmajor: float, the plasma major radius (m) + - inverse_q: np.ndarray, inverse safety factor profile + - sqeps: np.ndarray, the square root of the inverse aspect ratio profile + - tempi: np.ndarray, the ion temperature profile (keV) + - amain: np.ndarray, the atomic mass of the main ion species profile + - zeff: np.ndarray, the effective charge of the main ion species + - ni: np.ndarray, the ion density profile (/m^3) + + Returns: + - float, the ion collisionality + + Reference: + - None """ return ( 3.2e-6 - * nui(j, zmain, ni, tempi, amain) + * _ion_collisions_sauter(radial_elements, zeff, ni, tempi, amain) * rmajor / ( - np.abs(mu[j - 1] + 1.0e-4) - * sqeps[j - 1] ** 3 - * np.sqrt(tempi[j - 1] / amain[j - 1]) + np.abs(inverse_q[radial_elements - 1] + 1.0e-4) + * sqeps[radial_elements - 1] ** 3 + * np.sqrt(tempi[radial_elements - 1] / amain[radial_elements - 1]) ) ) @nb.jit(nopython=True, cache=True) -def dcsa(j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps): - """Grad(ln(ne)) coefficient in the Sauter bootstrap scaling - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 2 to nr - nr : input integer : maximum value of j - This function calculates the coefficient scaling grad(ln(ne)) - in the Sauter bootstrap current scaling. - Code by Angioni, 29th May 2002. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 - - DCSA $\\equiv \\mathcal{L}_{31}$, Eq.14a, Sauter et al, 1999 +def _calculate_l31_coefficient( + radial_elements: np.ndarray, + number_of_elements: int, + rmajor: float, + bt: float, + triang: float, + ne: np.ndarray, + ni: np.ndarray, + tempe: np.ndarray, + tempi: np.ndarray, + inverse_q: np.ndarray, + rho: np.ndarray, + zeff: np.ndarray, + sqeps: np.ndarray, +) -> float: + """ + L31 coefficient before Grad(ln(ne)) in the Sauter bootstrap scaling. + + Parameters: + - radial_elements: np.ndarray, radial element indexes in range 2 to nr + - number_of_elements: int, maximum value of radial_elements + - rmajor: float, plasma major radius (m) + - bt: float, toroidal field on axis (T) + - triang: float, plasma triangularity + - ne: np.ndarray, electron density profile (/m^3) + - ni: np.ndarray, ion density profile (/m^3) + - tempe: np.ndarray, electron temperature profile (keV) + - tempi: np.ndarray, ion temperature profile (keV) + - inverse_q: np.ndarray, inverse safety factor profile + - rho: np.ndarray, normalized minor radius profile + - zeff: np.ndarray, effective charge profile + - sqeps: np.ndarray, square root of inverse aspect ratio profile + + Returns: + - float, the coefficient scaling grad(ln(ne)) in the Sauter bootstrap current scaling. + + This function calculates the coefficient scaling grad(ln(ne)) in the Sauter bootstrap current scaling. + + Reference: + - O. Sauter, C. Angioni, Y. R. Lin-Liu; + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. + Phys. Plasmas 1 July 1999; 6 (7): 2834–2839. https://doi.org/10.1063/1.873240 + - O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime + [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 + """ + # Prevents first element being 0 + charge_profile = zeff[radial_elements - 1] + + # Calculate trapped particle fraction + f_trapped = _trapped_particle_fraction_sauter(radial_elements, triang, sqeps) + + # Calculated electron collisionality; nu_e* + electron_collisionality = _electron_collisionality_sauter( + radial_elements, rmajor, zeff, inverse_q, sqeps, tempe, ne + ) + + # $f^{31}_{teff}(\nu_{e*})$, Eq.14b + f31_teff = f_trapped / ( + (1.0 + (1.0 - 0.1 * f_trapped) * np.sqrt(electron_collisionality)) + + (0.5 * (1.0 - f_trapped) * electron_collisionality) / charge_profile + ) + + l31_coefficient = ( + ((1.0 + 1.4 / (charge_profile + 1.0)) * f31_teff) + - (1.9 / (charge_profile + 1.0) * f31_teff**2) + + ((0.3 * f31_teff**3 + 0.2 * f31_teff**4) / (charge_profile + 1.0)) + ) + + # Corrections suggested by Fable, 15/05/2015 + return l31_coefficient * _beta_poloidal_total_sauter( + radial_elements, + number_of_elements, + rmajor, + bt, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + ) + + +@nb.jit(nopython=True, cache=True) +def _calculate_l31_32_coefficient( + radial_elements: np.ndarray, + number_of_elements: int, + rmajor: float, + bt: float, + triang: float, + ne: np.ndarray, + ni: np.ndarray, + tempe: np.ndarray, + tempi: np.ndarray, + inverse_q: np.ndarray, + rho: np.ndarray, + zeff: np.ndarray, + sqeps: np.ndarray, +) -> float: + """ + L31 & L32 coefficient before Grad(ln(Te)) in the Sauter bootstrap scaling. + + Parameters: + - radial_elements: np.ndarray, radial element indexes in range 2 to nr + - number_of_elements: int, maximum value of radial_elements + - rmajor: float, plasma major radius (m) + - bt: float, toroidal field on axis (T) + - triang: float, plasma triangularity + - ne: np.ndarray, electron density profile (/m^3) + - ni: np.ndarray, ion density profile (/m^3) + - tempe: np.ndarray, electron temperature profile (keV) + - tempi: np.ndarray, ion temperature profile (keV) + - inverse_q: np.ndarray, inverse safety factor profile + - rho: np.ndarray, normalized minor radius profile + - zeff: np.ndarray, effective charge profile + - sqeps: np.ndarray, square root of inverse aspect ratio profile + + Returns: + - float, the L31 & L32 coefficient scaling grad(ln(Te)) in the Sauter bootstrap current scaling. + + This function calculates the coefficient scaling grad(ln(Te)) in the Sauter bootstrap current scaling. + + Reference: + - O. Sauter, C. Angioni, Y. R. Lin-Liu; + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. + Phys. Plasmas 1 July 1999; 6 (7): 2834–2839. https://doi.org/10.1063/1.873240 + - O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime + [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 """ - zz = zef[j - 1] - zft = tpf(j, triang, sqeps) - _nues = nues(j, rmajor, zef, mu, sqeps, tempe, ne) - zdf = (1.0 + (1.0 - 0.1 * zft) * np.sqrt(_nues)) + (0.5 * (1.0 - zft) * _nues) / zz - zft /= zdf # $f^{31}_{teff}(\nu_{e*})$, Eq.14b - zft2 = zft**2 - zft3 = zft2 * zft - dcsa = ( - ((1.0 + 1.4 / (zz + 1.0)) * zft) - - (1.9 / (zz + 1.0) * zft2) - + ((0.3 * zft3 + 0.2 * zft3 * zft) / (zz + 1.0)) + + # Prevents first element being 0 + charge_profile = zeff[radial_elements - 1] + + # Calculate trapped particle fraction + f_trapped = _trapped_particle_fraction_sauter(radial_elements, triang, sqeps) + + # Calculated electron collisionality; nu_e* + electron_collisionality = _electron_collisionality_sauter( + radial_elements, rmajor, zeff, inverse_q, sqeps, tempe, ne + ) + + # $f^{32\_ee}_{teff}(\nu_{e*})$, Eq.15d + f32ee_teff = f_trapped / ( + ( + 1.0 + + 0.26 * (1.0 - f_trapped) * np.sqrt(electron_collisionality) + + ( + 0.18 + * (1.0 - 0.37 * f_trapped) + * electron_collisionality + / np.sqrt(charge_profile) + ) + ) + ) + + # $f^{32\_ei}_{teff}(\nu_{e*})$, Eq.15e + f32ei_teff = f_trapped / ( + (1.0 + (1.0 + 0.6 * f_trapped) * np.sqrt(electron_collisionality)) + + ( + 0.85 + * (1.0 - 0.37 * f_trapped) + * electron_collisionality + * (1.0 + charge_profile) + ) + ) + + # $F_{32\_ee}(X)$, Eq.15b + big_f32ee_teff = ( + ( + (0.05 + 0.62 * charge_profile) + / charge_profile + / (1.0 + 0.44 * charge_profile) + * (f32ee_teff - f32ee_teff**4) + ) + + ( + ( + f32ee_teff**2 + - f32ee_teff**4 + - 1.2 * (f32ee_teff**3 - f32ee_teff**4) + ) + / (1.0 + 0.22 * charge_profile) + ) + + (1.2 / (1.0 + 0.5 * charge_profile) * f32ee_teff**4) + ) + + # $F_{32\_ei}(Y)$, Eq.15c + big_f32ei_teff = ( + ( + -(0.56 + 1.93 * charge_profile) + / charge_profile + / (1.0 + 0.44 * charge_profile) + * (f32ei_teff - f32ei_teff**4) + ) + + ( + 4.95 + / (1.0 + 2.48 * charge_profile) + * ( + f32ei_teff**2 + - f32ei_teff**4 + - 0.55 * (f32ei_teff**3 - f32ei_teff**4) + ) + ) + - (1.2 / (1.0 + 0.5 * charge_profile) * f32ei_teff**4) ) + # big_f32ee_teff + big_f32ei_teff = L32 coefficient + # Corrections suggested by Fable, 15/05/2015 - return dcsa * beta_poloidal_local_total( - j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + return _beta_poloidal_sauter( + radial_elements, number_of_elements, rmajor, bt, ne, tempe, inverse_q, rho + ) * (big_f32ee_teff + big_f32ei_teff) + _calculate_l31_coefficient( + radial_elements, + number_of_elements, + rmajor, + bt, + triang, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + zeff, + sqeps, + ) * _beta_poloidal_sauter( + radial_elements, number_of_elements, rmajor, bt, ne, tempe, inverse_q, rho + ) / _beta_poloidal_total_sauter( + radial_elements, + number_of_elements, + rmajor, + bt, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, ) @nb.jit(nopython=True, cache=True) -def hcsa(j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps): - """Grad(ln(Te)) coefficient in the Sauter bootstrap scaling - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 2 to nr - nr : input integer : maximum value of j - This function calculates the coefficient scaling grad(ln(Te)) - in the Sauter bootstrap current scaling. - Code by Angioni, 29th May 2002. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 +def _calculate_l34_alpha_31_coefficient( + radial_elements: np.ndarray, + number_of_elements: int, + rmajor: float, + bt: float, + triang: float, + inverse_q: np.ndarray, + sqeps: np.ndarray, + tempi: np.ndarray, + tempe: np.ndarray, + amain: float, + zmain: float, + ni: np.ndarray, + ne: np.ndarray, + rho: np.ndarray, + zeff: np.ndarray, +) -> float: """ - zz = zef[j - 1] - zft = tpf(j, triang, sqeps) - _nues = nues(j, rmajor, zef, mu, sqeps, tempe, ne) - _nues_sqrt = np.sqrt(_nues) - zdf = (1.0 + 0.26 * (1.0 - zft) * _nues_sqrt) + ( - 0.18 * (1.0 - 0.37 * zft) * _nues / np.sqrt(zz) - ) - zfte = zft / zdf # $f^{32\_ee}_{teff}(\nu_{e*})$, Eq.15d - zfte2 = zfte * zfte - zfte3 = zfte * zfte2 - zfte4 = zfte2 * zfte2 - - zdf = (1.0 + (1.0 + 0.6 * zft) * _nues_sqrt) + ( - 0.85 * (1.0 - 0.37 * zft) * _nues * (1.0 + zz) - ) - zfti = zft / zdf # $f^{32\_ei}_{teff}(\nu_{e*})$, Eq.15e - zfti2 = zfti * zfti - zfti3 = zfti * zfti2 - zfti4 = zfti2 * zfti2 + L34, alpha and L31 coefficient before Grad(ln(Ti)) in the Sauter bootstrap scaling. + + Parameters: + - radial_elements: np.ndarray, radial element indexes in range 2 to nr + - number_of_elements: int, maximum value of radial_elements + - rmajor: float, plasma major radius (m) + - bt: float, toroidal field on axis (T) + - triang: float, plasma triangularity + - inverse_q: np.ndarray, inverse safety factor profile + - sqeps: np.ndarray, square root of inverse aspect ratio profile + - tempi: np.ndarray, ion temperature profile (keV) + - tempe: np.ndarray, electron temperature profile (keV) + - amain: float, atomic mass of the main ion + - zmain: float, charge of the main ion + - ni: np.ndarray, ion density profile (/m^3) + - ne: np.ndarray, electron density profile (/m^3) + - rho: np.ndarray, normalized minor radius profile + - zeff: np.ndarray, effective charge profile + + Returns: + - float, the the L34, alpha and L31 coefficient scaling grad(ln(Ti)) in the Sauter bootstrap current scaling. + + This function calculates the coefficient scaling grad(ln(Ti)) in the Sauter bootstrap current scaling. + + Reference: + - O. Sauter, C. Angioni, Y. R. Lin-Liu; + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. + Phys. Plasmas 1 July 1999; 6 (7): 2834–2839. https://doi.org/10.1063/1.873240 + - O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime + [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 + """ + # Prevents first element being 0 + charge_profile = zeff[radial_elements - 1] - # $F_{32\_ee}(X)$, Eq.15b - hcee = ( - ((0.05 + 0.62 * zz) / zz / (1.0 + 0.44 * zz) * (zfte - zfte4)) - + ((zfte2 - zfte4 - 1.2 * (zfte3 - zfte4)) / (1.0 + 0.22 * zz)) - + (1.2 / (1.0 + 0.5 * zz) * zfte4) - ) + # Calculate trapped particle fraction + f_trapped = _trapped_particle_fraction_sauter(radial_elements, triang, sqeps) - # $F_{32\_ei}(Y)$, Eq.15c - hcei = ( - (-(0.56 + 1.93 * zz) / zz / (1.0 + 0.44 * zz) * (zfti - zfti4)) - + (4.95 / (1.0 + 2.48 * zz) * (zfti2 - zfti4 - 0.55 * (zfti3 - zfti4))) - - (1.2 / (1.0 + 0.5 * zz) * zfti4) + # Calculated electron collisionality; nu_e* + electron_collisionality = _electron_collisionality_sauter( + radial_elements, rmajor, zeff, inverse_q, sqeps, tempe, ne ) - # Corrections suggested by Fable, 15/05/2015 - return beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) * ( - hcee + hcei - ) + dcsa( - j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps - ) * beta_poloidal_local( - j, nr, rmajor, bt, ne, tempe, mu, rho - ) / beta_poloidal_local_total( - j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + # $f^{34}_{teff}(\nu_{e*})$, Eq.16b + f34_teff = f_trapped / ( + (1.0 + (1.0 - 0.1 * f_trapped) * np.sqrt(electron_collisionality)) + + 0.5 * (1.0 - 0.5 * f_trapped) * electron_collisionality / charge_profile ) - -@nb.jit(nopython=True, cache=True) -def xcsa( - j, - nr, - rmajor, - bt, - triang, - mu, - sqeps, - tempi, - tempe, - amain, - zmain, - ni, - ne, - rho, - zef, -): - """Grad(ln(Ti)) coefficient in the Sauter bootstrap scaling - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 2 to nr - nr : input integer : maximum value of j - This function calculates the coefficient scaling grad(ln(Ti)) - in the Sauter bootstrap current scaling. - Code by Angioni, 29th May 2002. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 - """ - zz = zef[j - 1] - zft = tpf(j, triang, sqeps) - _nues = nues(j, rmajor, zef, mu, sqeps, tempe, ne) - zdf = (1.0 + (1.0 - 0.1 * zft) * np.sqrt(_nues)) + 0.5 * ( - 1.0 - 0.5 * zft - ) * _nues / zz - zfte = zft / zdf # $f^{34}_{teff}(\nu_{e*})$, Eq.16b - # Eq.16a - xcsa = ((1.0 + 1.4 / (zz + 1.0)) * zfte - 1.9 / (zz + 1.0) * zfte * zfte) + ( - 0.3 * zfte * zfte + 0.2 * zfte * zfte * zfte - ) * zfte / (zz + 1.0) + l34_coefficient = ( + ((1.0 + (1.4 / (charge_profile + 1.0))) * f34_teff) + - ((1.9 / (charge_profile + 1.0)) * f34_teff**2) + + ((0.3 / (charge_profile + 1.0)) * f34_teff**3) + + ((0.2 / (charge_profile + 1.0)) * f34_teff**4) + ) # $\alpha_0$, Eq.17a - a0 = (-1.17 * (1.0 - zft)) / (1.0 - 0.22 * zft - 0.19 * zft * zft) + alpha_0 = (-1.17 * (1.0 - f_trapped)) / ( + 1.0 - (0.22 * f_trapped) - 0.19 * f_trapped**2 + ) - _nuis = nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni) - _nuis_sqrt = np.sqrt(_nuis) - a1 = _nuis**2 * zft**6 + # Calculate the ion collisionality + ion_collisionality = _ion_collisionality_sauter( + radial_elements, rmajor, inverse_q, sqeps, tempi, amain, zmain, ni + ) # $\alpha(\nu_{i*})$, Eq.17b - alp = ( - (a0 + 0.25 * (1.0 - zft * zft) * _nuis_sqrt) / (1.0 + 0.5 * _nuis_sqrt) - + 0.315 * a1 - ) / (1.0 + 0.15 * a1) + alpha = ( + ( + (alpha_0 + (0.25 * (1.0 - f_trapped**2)) * np.sqrt(ion_collisionality)) + / (1.0 + (0.5 * np.sqrt(ion_collisionality))) + + (0.315 * ion_collisionality**2 * f_trapped**6) + ) + ) / (1.0 + (0.15 * ion_collisionality**2 * f_trapped**6)) # Corrections suggested by Fable, 15/05/2015 + # Below calculates the L34 * alpha + L31 coefficient return ( - beta_poloidal_local_total(j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho) - - beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) - ) * (xcsa * alp) + dcsa( - j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps + _beta_poloidal_total_sauter( + radial_elements, + number_of_elements, + rmajor, + bt, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + ) + - _beta_poloidal_sauter( + radial_elements, + number_of_elements, + rmajor, + bt, + ne, + tempe, + inverse_q, + rho, + ) + ) * (l34_coefficient * alpha) + _calculate_l31_coefficient( + radial_elements, + number_of_elements, + rmajor, + bt, + triang, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + zeff, + sqeps, ) * ( 1.0 - - beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) - / beta_poloidal_local_total(j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho) + - _beta_poloidal_sauter( + radial_elements, + number_of_elements, + rmajor, + bt, + ne, + tempe, + inverse_q, + rho, + ) + / _beta_poloidal_total_sauter( + radial_elements, + number_of_elements, + rmajor, + bt, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + ) ) @nb.jit(nopython=True, cache=True) -def nues(j, rmajor, zef, mu, sqeps, tempe, ne): - """Relative frequency of electron collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the relative frequency of electron - collisions: NU* = Nuei*q*Rt/eps**1.5/Vte - The electron-ion collision frequency NUEI=NUEE*1.4*ZEF is - used. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - Yushmanov, 30th April 1987 (?) +def _beta_poloidal_sauter( + radial_elements: np.ndarray, + nr: int, + rmajor: float, + bt: float, + ne: np.ndarray, + tempe: np.ndarray, + inverse_q: np.ndarray, + rho: np.ndarray, +) -> np.ndarray: """ - return ( - nuee(j, tempe, ne) - * 1.4 - * zef[j - 1] - * rmajor - / np.abs(mu[j - 1] * (sqeps[j - 1] ** 3) * np.sqrt(tempe[j - 1]) * 1.875e7) - ) - - -@nb.jit(nopython=True, cache=True) -def beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho): - """Local beta poloidal calculation - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the local beta poloidal. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). -

beta poloidal = 4*pi*ne*Te/Bpo**2 - Pereverzev, 25th April 1989 (?) + Calculate the local beta poloidal using only electron profiles for the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, radial element indexes in range 1 to nr + - nr: int, maximum value of radial_elements + - rmajor: float, plasma major radius (m) + - bt: float, toroidal field on axis (T) + - ne: np.ndarray, electron density profile (/m^3) + - tempe: np.ndarray, electron temperature profile (keV) + - inverse_q: np.ndarray, inverse safety factor profile + - rho: np.ndarray, normalized minor radius profile + + Returns: + - np.ndarray, the local beta poloidal + + Reference: + - None """ return ( np.where( - j != nr, - 1.6e-4 * np.pi * (ne[j] + ne[j - 1]) * (tempe[j] + tempe[j - 1]), - 6.4e-4 * np.pi * ne[j - 1] * tempe[j - 1], + radial_elements != nr, + 1.6e-4 + * np.pi + * (ne[radial_elements] + ne[radial_elements - 1]) + * (tempe[radial_elements] + tempe[radial_elements - 1]), + 6.4e-4 * np.pi * ne[radial_elements - 1] * tempe[radial_elements - 1], ) - * (rmajor / (bt * rho[j - 1] * np.abs(mu[j - 1] + 1.0e-4))) ** 2 + * ( + rmajor + / ( + bt + * rho[radial_elements - 1] + * np.abs(inverse_q[radial_elements - 1] + 1.0e-4) + ) + ) + ** 2 ) @nb.jit(nopython=True, cache=True) -def beta_poloidal_local_total(j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho): - """Local beta poloidal calculation, including ion pressure - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the local total beta poloidal. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). -

beta poloidal = 4*pi*(ne*Te+ni*Ti)/Bpo**2 - where ni is the sum of all ion densities (thermal) - Pereverzev, 25th April 1989 (?) - E Fable, private communication, 15th May 2014 +def _beta_poloidal_total_sauter( + radial_elements: np.ndarray, + nr: int, + rmajor: float, + bt: float, + ne: np.ndarray, + ni: np.ndarray, + tempe: np.ndarray, + tempi: np.ndarray, + inverse_q: np.ndarray, + rho: np.ndarray, +) -> np.ndarray: + """ + Calculate the local beta poloidal including ion pressure for the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: np.ndarray, radial element indexes in range 1 to nr + - nr: int, maximum value of radial_elements + - rmajor: float, plasma major radius (m) + - bt: float, toroidal field on axis (T) + - ne: np.ndarray, electron density profile (/m^3) + - ni: np.ndarray, ion density profile (/m^3) + - tempe: np.ndarray, electron temperature profile (keV) + - tempi: np.ndarray, ion temperature profile (keV) + - inverse_q: np.ndarray, inverse safety factor profile + - rho: np.ndarray, normalized minor radius profile + + Returns: + - np.ndarray, the local total beta poloidal + + Reference: + - None """ - return ( np.where( - j != nr, + radial_elements != nr, 1.6e-4 * np.pi * ( - ((ne[j] + ne[j - 1]) * (tempe[j] + tempe[j - 1])) - + ((ni[j] + ni[j - 1]) * (tempi[j] + tempi[j - 1])) + ( + (ne[radial_elements] + ne[radial_elements - 1]) + * (tempe[radial_elements] + tempe[radial_elements - 1]) + ) + + ( + (ni[radial_elements] + ni[radial_elements - 1]) + * (tempi[radial_elements] + tempi[radial_elements - 1]) + ) + ), + 6.4e-4 + * np.pi + * ( + ne[radial_elements - 1] * tempe[radial_elements - 1] + + ni[radial_elements - 1] * tempi[radial_elements - 1] ), - 6.4e-4 * np.pi * (ne[j - 1] * tempe[j - 1] + ni[j - 1] * tempi[j - 1]), ) - * (rmajor / (bt * rho[j - 1] * np.abs(mu[j - 1] + 1.0e-4))) ** 2 + * ( + rmajor + / ( + bt + * rho[radial_elements - 1] + * np.abs(inverse_q[radial_elements - 1] + 1.0e-4) + ) + ) + ** 2 ) @nb.jit(nopython=True, cache=True) -def tpf(j, triang, sqeps, fit=1): - """Trapped particle fraction - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - fit : input integer : (1)=ASTRA method, 2=Equation from Sauter2002, 3=Equation from Sauter2013 - This function calculates the trapped particle fraction at - a given radius. -

A number of different fits are provided, but the one - to be used is hardwired prior to run-time. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter et al, Plasma Phys. Contr. Fusion 44 (2002) 1999 - O. Sauter, 2013: - http://infoscience.epfl.ch/record/187521/files/lrp_012013.pdf +def _trapped_particle_fraction_sauter( + radial_elements: np.ndarray, triang: float, sqeps: np.ndarray, fit: int = 0 +) -> np.ndarray: + """ + Calculates the trapped particle fraction to be used in the Sauter bootstrap current scaling. + + Parameters: + - radial_elements: list, radial element index in range 1 to nr + - triang: float, plasma triangularity + - sqeps: list, square root of local aspect ratio + - fit: int, fit method (1 = ASTRA method, 2 = Equation from Sauter 2002, 3 = Equation from Sauter 2016) + + Returns: + - list, trapped particle fraction + + This function calculates the trapped particle fraction at a given radius. + + References: + Used in this paper: + - O. Sauter, C. Angioni, Y. R. Lin-Liu; + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. + Phys. Plasmas 1 July 1999; 6 (7): 2834–2839. https://doi.org/10.1063/1.873240 + + - O. Sauter, R. J. Buttery, R. Felton, T. C. Hender, D. F. Howell, and contributors to the E.-J. Workprogramme, + “Marginal  -limit for neoclassical tearing modes in JET H-mode discharges,” + Plasma Physics and Controlled Fusion, vol. 44, no. 9, pp. 1999–2019, Aug. 2002, + doi: https://doi.org/10.1088/0741-3335/44/9/315. + + - O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, + Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, ISSN 0920-3796, + https://doi.org/10.1016/j.fusengdes.2016.04.033. + """ - s = sqeps[j - 1] - eps = s * s + # Prevent first element from being zero + sqeps_reduced = sqeps[radial_elements - 1] + eps = sqeps_reduced**2 - if fit == 1: + if fit == 0: # ASTRA method, from Emiliano Fable, private communication # (Excluding h term which dominates for inverse aspect ratios < 0.5, # and tends to take the trapped particle fraction to 1) zz = 1.0 - eps - return 1.0 - zz * np.sqrt(zz) / (1.0 + 1.46 * s) - elif fit == 2: - # Equation 4 of Sauter 2002 - # Similar to, but not quite identical to g above + return 1.0 - zz * np.sqrt(zz) / (1.0 + 1.46 * sqeps_reduced) + + elif fit == 1: + # Equation 4 of Sauter 2002; https://doi.org/10.1088/0741-3335/44/9/315. + # Similar to, but not quite identical to above + + return 1.0 - ( + ((1.0 - eps) ** 2) + / ((1.0 + 1.46 * sqeps_reduced) * np.sqrt(1.0 - eps**2)) + ) - return 1.0 - (1.0 - eps) ** 2 / (1.0 + 1.46 * s) / np.sqrt(1.0 - eps * eps) - elif fit == 3: + elif fit == 2: + # Sauter 2016; https://doi.org/10.1016/j.fusengdes.2016.04.033. # Includes correction for triangularity epseff = 0.67 * (1.0 - 1.4 * triang * np.abs(triang)) * eps - return 1.0 - np.sqrt((1.0 - eps) / (1.0 + eps)) * (1.0 - epseff) / ( - 1.0 + 2.0 * np.sqrt(epseff) + return 1.0 - ( + ((1.0 - epseff) / (1.0 + 2.0 * np.sqrt(epseff))) + * (np.sqrt((1.0 - eps) / (1.0 + eps))) ) raise RuntimeError(f"fit={fit} is not valid. Must be 1, 2, or 3") @@ -778,32 +1475,18 @@ def physics(self): Routine to calculate tokamak plasma physics information author: P J Knight, CCFE, Culham Science Centre None - This routine calculates all the primary plasma physics - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - - Part 1: Physics https://www.sciencedirect.com/science/article/pii/S0920379614005961 - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO - https://iopscience.iop.org/article/10.1088/0029-5515/53/7/073019 - T. Hartmann, 2013, Development of a modular systems code to analyse the - implications of physics assumptions on the design of a demonstration fusion power plant - https://inis.iaea.org/search/search.aspx?orig_q=RN:45031642 + This routine calculates all the primary plasma physics parameters for a tokamak fusion reactor. + + References: + - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - Part 1: Physics + https://www.sciencedirect.com/science/article/pii/S0920379614005961 + - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO + https://iopscience.iop.org/article/10.1088/0029-5515/53/7/073019 + - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant + https://inis.iaea.org/search/search.aspx?orig_q=RN:45031642 """ - # kappaa_IPB = physics_variables.vol / ( - # 2.0e0 - # * numpy.pi - # * numpy.pi - # * physics_variables.rminor - # * physics_variables.rminor - # * physics_variables.rmajor - # ) - - if physics_variables.icurr == 2: - physics_variables.q95 = ( - physics_variables.q * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 - ) - else: - physics_variables.q95 = ( - physics_variables.q - ) # i.e. input (or iteration variable) value + + physics_variables.q95 = physics_variables.q # Calculate plasma composition # Issue #261 Remove old radiation model (imprad_model=0) @@ -815,13 +1498,13 @@ def physics(self): physics_variables.rli, physics_variables.bp, physics_variables.qstar, - physics_variables.plascur, - ) = self.culcur( + physics_variables.plasma_current, + ) = self.calculate_plasma_current( physics_variables.alphaj, physics_variables.alphap, physics_variables.bt, physics_variables.eps, - physics_variables.icurr, + physics_variables.i_plasma_current, physics_variables.iprofile, physics_variables.kappa, physics_variables.kappa95, @@ -845,7 +1528,7 @@ def physics(self): physics_variables.neped = ( physics_variables.fgwped * 1.0e14 - * physics_variables.plascur + * physics_variables.plasma_current / (np.pi * physics_variables.rminor * physics_variables.rminor) ) @@ -853,7 +1536,7 @@ def physics(self): physics_variables.nesep = ( physics_variables.fgwsep * 1.0e14 - * physics_variables.plascur + * physics_variables.plasma_current / (np.pi * physics_variables.rminor * physics_variables.rminor) ) @@ -868,7 +1551,7 @@ def physics(self): # Set PF coil ramp times if pulse_variables.lpulse != 1: if times_variables.tohsin == 0.0e0: - times_variables.tohs = physics_variables.plascur / 5.0e5 + times_variables.tohs = physics_variables.plasma_current / 5.0e5 times_variables.tramp = times_variables.tohs times_variables.tqnch = times_variables.tohs else: @@ -877,7 +1560,7 @@ def physics(self): else: if times_variables.pulsetimings == 0.0e0: # times_variables.tramp is input - times_variables.tohs = physics_variables.plascur / 1.0e5 + times_variables.tohs = physics_variables.plasma_current / 1.0e5 times_variables.tqnch = times_variables.tohs else: @@ -918,24 +1601,66 @@ def physics(self): + times_variables.tdwell ) - # Calculate bootstrap current fraction using various models - current_drive_variables.bscf_iter89 = self.bootstrap_fraction_iter89( - physics_variables.aspect, - physics_variables.beta, - physics_variables.btot, - current_drive_variables.cboot, - physics_variables.plascur, - physics_variables.q95, - physics_variables.q0, - physics_variables.rmajor, - physics_variables.vol, + # ***************************** # + # DIAMAGNETIC CURRENT # + # ***************************** # + + # Hender scaling for diamagnetic current at tight physics_variables.aspect ratio + current_drive_variables.diacf_hender = diamagnetic_fraction_hender( + physics_variables.beta ) + # SCENE scaling for diamagnetic current + current_drive_variables.diacf_scene = diamagnetic_fraction_scene( + physics_variables.beta, physics_variables.q95, physics_variables.q0 + ) + + if physics_variables.i_diamagnetic_current == 1: + current_drive_variables.diamagnetic_current_fraction = ( + current_drive_variables.diacf_hender + ) + elif physics_variables.i_diamagnetic_current == 2: + current_drive_variables.diamagnetic_current_fraction = ( + current_drive_variables.diacf_scene + ) + + # ***************************** # + # PFIRSCH-SCHLÜTER CURRENT # + # ***************************** # + + # Pfirsch-Schlüter scaling for diamagnetic current + current_drive_variables.pscf_scene = ps_fraction_scene(physics_variables.beta) + + if physics_variables.i_pfirsch_schluter_current == 1: + current_drive_variables.ps_current_fraction = ( + current_drive_variables.pscf_scene + ) + + # ***************************** # + # BOOTSTRAP CURRENT # + # ***************************** # + + # Calculate bootstrap current fraction using various models + current_drive_variables.bscf_iter89 = ( + current_drive_variables.cboot + * self.bootstrap_fraction_iter89( + physics_variables.aspect, + physics_variables.beta, + physics_variables.btot, + physics_variables.plasma_current, + physics_variables.q95, + physics_variables.q0, + physics_variables.rmajor, + physics_variables.vol, + ) + ) + # Calculate the toroidal beta for the Nevins scaling betat = ( physics_variables.beta * physics_variables.btot**2 / physics_variables.bt**2 ) + current_drive_variables.bscf_nevins = ( current_drive_variables.cboot * self.bootstrap_fraction_nevins( @@ -944,12 +1669,12 @@ def physics(self): betat, physics_variables.bt, physics_variables.dene, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.q95, physics_variables.q0, physics_variables.rmajor, physics_variables.rminor, - physics_variables.ten, + physics_variables.te, physics_variables.zeff, ) ) @@ -972,19 +1697,6 @@ def physics(self): ) ) - # Hender scaling for diamagnetic current at tight physics_variables.aspect ratio - current_drive_variables.diacf_hender = diamagnetic_fraction_hender( - physics_variables.beta - ) - - # SCENE scaling for diamagnetic current - current_drive_variables.diacf_scene = diamagnetic_fraction_scene( - physics_variables.beta, physics_variables.q95, physics_variables.q0 - ) - - # Pfirsch-Schlüter scaling for diamagnetic current - current_drive_variables.pscf_scene = ps_fraction_scene(physics_variables.beta) - current_drive_variables.bscf_sauter = ( current_drive_variables.cboot * self.bootstrap_fraction_sauter(self.plasma_profile) @@ -1003,45 +1715,55 @@ def physics(self): ) ) - if current_drive_variables.bscfmax < 0.0e0: - current_drive_variables.bootipf = abs(current_drive_variables.bscfmax) - current_drive_variables.plasipf = current_drive_variables.bootipf + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: + current_drive_variables.bootstrap_current_fraction = abs( + current_drive_variables.bootstrap_current_fraction_max + ) + current_drive_variables.plasma_current_internal_fraction = ( + current_drive_variables.bootstrap_current_fraction + ) else: - if physics_variables.ibss == 1: - current_drive_variables.bootipf = current_drive_variables.bscf_iter89 - elif physics_variables.ibss == 2: - current_drive_variables.bootipf = current_drive_variables.bscf_nevins - elif physics_variables.ibss == 3: - 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: + if physics_variables.i_bootstrap_current == 1: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_iter89 + ) + elif physics_variables.i_bootstrap_current == 2: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_nevins + ) + elif physics_variables.i_bootstrap_current == 3: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_wilson + ) + elif physics_variables.i_bootstrap_current == 4: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_sauter + ) + elif physics_variables.i_bootstrap_current == 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 + # So the diamagnetic current calculation should be turned off when using, (i_diamagnetic_current = 0). + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_sakai + ) else: - error_handling.idiags[0] = physics_variables.ibss + error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) physics_module.err242 = 0 - if current_drive_variables.bootipf > current_drive_variables.bscfmax: - current_drive_variables.bootipf = min( - current_drive_variables.bootipf, current_drive_variables.bscfmax + if ( + current_drive_variables.bootstrap_current_fraction + > current_drive_variables.bootstrap_current_fraction_max + ): + current_drive_variables.bootstrap_current_fraction = min( + current_drive_variables.bootstrap_current_fraction, + current_drive_variables.bootstrap_current_fraction_max, ) physics_module.err242 = 1 - if physics_variables.idia == 1: - current_drive_variables.diaipf = current_drive_variables.diacf_hender - elif physics_variables.idia == 2: - current_drive_variables.diaipf = current_drive_variables.diacf_scene - - if physics_variables.ips == 1: - current_drive_variables.psipf = current_drive_variables.pscf_scene - - current_drive_variables.plasipf = ( - current_drive_variables.bootipf - + current_drive_variables.diaipf - + current_drive_variables.psipf + current_drive_variables.plasma_current_internal_fraction = ( + current_drive_variables.bootstrap_current_fraction + + current_drive_variables.diamagnetic_current_fraction + + current_drive_variables.ps_current_fraction ) # Plasma driven current fraction (Bootstrap + Diamagnetic @@ -1050,17 +1772,24 @@ def physics(self): # produced by non-inductive means (which also includes # the current drive proportion) physics_module.err243 = 0 - if current_drive_variables.plasipf > physics_variables.fvsbrnni: - current_drive_variables.plasipf = min( - current_drive_variables.plasipf, physics_variables.fvsbrnni + if ( + current_drive_variables.plasma_current_internal_fraction + > physics_variables.fvsbrnni + ): + current_drive_variables.plasma_current_internal_fraction = min( + current_drive_variables.plasma_current_internal_fraction, + physics_variables.fvsbrnni, ) physics_module.err243 = 1 # Fraction of plasma current produced by inductive means - physics_variables.facoh = max(1.0e-10, (1.0e0 - physics_variables.fvsbrnni)) + physics_variables.inductive_current_fraction = max( + 1.0e-10, (1.0e0 - physics_variables.fvsbrnni) + ) # Fraction of plasma current produced by auxiliary current drive - physics_variables.faccd = ( - physics_variables.fvsbrnni - current_drive_variables.plasipf + physics_variables.aux_current_fraction = ( + physics_variables.fvsbrnni + - current_drive_variables.plasma_current_internal_fraction ) # Auxiliary current drive power calculations @@ -1214,9 +1943,9 @@ def physics(self): physics_variables.rpfac, physics_variables.rplas, ) = self.pohm( - physics_variables.facoh, + physics_variables.inductive_current_fraction, physics_variables.kappa95, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.rmajor, physics_variables.rminor, physics_variables.ten, @@ -1288,7 +2017,7 @@ def physics(self): physics_variables.bt, physics_variables.idensl, physics_variables.pdivt, - physics_variables.plascur, + physics_variables.plasma_current, divertor_variables.prn1, physics_variables.qstar, physics_variables.q95, @@ -1325,7 +2054,7 @@ def physics(self): physics_variables.kappa95, physics_variables.pchargemw, current_drive_variables.pinjmw, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.pcoreradpv, physics_variables.rmajor, physics_variables.rminor, @@ -1355,12 +2084,12 @@ def physics(self): ) = vscalc( physics_variables.csawth, physics_variables.eps, - physics_variables.facoh, + physics_variables.inductive_current_fraction, physics_variables.gamma, physics_variables.kappa, physics_variables.rmajor, physics_variables.rplas, - physics_variables.plascur, + physics_variables.plasma_current, times_variables.t_fusion_ramp, times_variables.tburn, physics_variables.rli, @@ -1383,7 +2112,7 @@ def physics(self): physics_variables.deni, physics_variables.fusionrate, physics_variables.alpharate, - physics_variables.plascur, + physics_variables.plasma_current, sbar, physics_variables.dnalp, physics_variables.taueff, @@ -1431,7 +2160,7 @@ def physics(self): physics_variables.betalim = culblm( physics_variables.bt, physics_variables.dnbeta, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.rminor, ) @@ -1591,7 +2320,7 @@ def physics(self): @staticmethod def culdlm( - bt, idensl, pdivt, plascur, prn1, qcyl, q95, rmajor, rminor, sarea, zeff + bt, idensl, pdivt, plasma_current, prn1, qcyl, q95, rmajor, rminor, sarea, zeff ): """Density limit calculation author: P J Knight, CCFE, Culham Science Centre @@ -1599,7 +2328,7 @@ def culdlm( idensl : input/output integer : switch denoting which formula to enforce pdivt : input real : power flowing to the edge plasma via charged particles (MW) - plascur : input real : plasma current (A) + plasma_current : input real : plasma current (A) prn1 : input real : edge density / average plasma density qcyl : input real : equivalent cylindrical safety factor (qstar) q95 : input real : safety factor at 95% surface @@ -1682,7 +2411,7 @@ def culdlm( # Greenwald limit - dlimit[6] = 1.0e14 * plascur / (np.pi * rminor * rminor) + dlimit[6] = 1.0e14 * plasma_current / (np.pi * rminor * rminor) # Enforce the chosen density limit @@ -1904,7 +2633,7 @@ def phyaux( deni, fusionrate, alpharate, - plascur, + plasma_current, sbar, dnalp, taueff, @@ -1918,7 +2647,7 @@ def phyaux( dnalp : input real : alpha ash density (/m3) fusionrate : input real : fusion reaction rate (/m3/s) alpharate : input real : alpha particle production rate (/m3/s) - plascur: input real : plasma current (A) + plasma_current: input real : plasma current (A) sbar : input real : exponent for aspect ratio (normally 1) taueff : input real : global energy confinement time (s) vol : input real : plasma volume (m3) @@ -1933,7 +2662,7 @@ def phyaux( needed by other parts of the code """ - figmer = 1e-6 * plascur * aspect**sbar + figmer = 1e-6 * plasma_current * aspect**sbar dntau = taueff * dene @@ -1976,7 +2705,16 @@ def phyaux( return burnup, dntau, figmer, fusrat, qfuel, rndfuel, taup @staticmethod - def pohm(facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff): + def pohm( + inductive_current_fraction, + kappa95, + plasma_current, + rmajor, + rminor, + ten, + vol, + zeff, + ): # Density weighted electron temperature in 10 keV units t10 = ten / 10.0 @@ -2007,9 +2745,9 @@ def pohm(facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff): error_handling.report_error(83) # Ohmic heating power per unit volume - # Corrected from: pohmpv = (facoh*plascur)**2 * ... + # Corrected from: pohmpv = (inductive_current_fraction*plasma_current)**2 * ... - pohmpv = facoh * plascur**2 * rplas * 1.0e-6 / vol + pohmpv = inductive_current_fraction * plasma_current**2 * rplas * 1.0e-6 / vol # Total ohmic heating power @@ -2018,171 +2756,183 @@ def pohm(facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff): return pohmpv, pohmmw, rpfac, rplas @staticmethod - def culcur( - alphaj, - alphap, - bt, - eps, - icurr, - iprofile, - kappa, - kappa95, - p0, - pperim, - q0, - q95, - rli, - rmajor, - rminor, - sf, - triang, - triang95, - ): - """Routine to calculate the plasma current - author: P J Knight, CCFE, Culham Science Centre - alphaj : input/output real : current profile index - alphap : input real : pressure profile index - bt : input real : toroidal field on axis (T) - eps : input real : inverse aspect ratio - icurr : input integer : current scaling model to use - 1 = Peng analytic fit - 2 = Peng divertor scaling (TART) - 3 = simple ITER scaling - 4 = revised ITER scaling - 5 = Todd empirical scaling I - 6 = Todd empirical scaling II - 7 = Connor-Hastie model - 8 = Sauter scaling (allowing negative triangularity) Issue #392 - 'Geometric formulas for system codes including the effect of negative triangularity' - iprofile : input integer : switch for current profile consistency - 0 use input values for alphaj, rli, dnbeta - 1 make these consistent with input q, q_0 values (recommend `icurr=4` with this option) - 2 use input values for alphaj, rli. Scale dnbeta with aspect ratio (original scaling) - 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 - kappa : input real : plasma elongation - kappa95 : input real : plasma elongation at 95% surface - p0 : input real : central plasma pressure (Pa) - pperim : input real : plasma perimeter length (m) - q0 : input real : plasma safety factor on axis - q95 : input real : plasma safety factor at 95% flux (= q-bar for icurr=2) - rli : input/output real : plasma normalised internal inductance - rmajor : input real : major radius (m) - rminor : input real : minor radius (m) - sf : input real : shape factor for icurr=1 (=A/pi in documentation) - triang : input real : plasma triangularity - triang95 : input real : plasma triangularity at 95% surface - bp : output real : poloidal field (T) - qstar : output real : equivalent cylindrical safety factor (shaped) - plascur : output real : plasma current (A) - This routine calculates the plasma current based on the edge - safety factor q95. It will also make the current profile parameters - consistent with the q-profile if required. - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, - Fusion Technology, 21, 1729 - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - - Part 1: Physics https://www.sciencedirect.com/science/article/pii/S0920379614005961 - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO - https://iopscience.iop.org/article/10.1088/0029-5515/53/7/073019 - T. Hartmann, 2013, Development of a modular systems code to analyse the - implications of physics assumptions on the design of a demonstration fusion power plant - https://inis.iaea.org/search/search.aspx?orig_q=RN:45031642 - Sauter, Geometric formulas for systems codes..., FED 2016 + def calculate_plasma_current( + alphaj: float, + alphap: float, + bt: float, + eps: float, + i_plasma_current: int, + iprofile: int, + kappa: float, + kappa95: float, + p0: float, + pperim: float, + q0: float, + q95: float, + rli: float, + rmajor: float, + rminor: float, + sf: float, + triang: float, + triang95: float, + ) -> Tuple[float, float, float, float, float]: + """Calculate the plasma current. + + Args: + alphaj (float): Current profile index. + alphap (float): Pressure profile index. + bt (float): Toroidal field on axis (T). + eps (float): Inverse aspect ratio. + i_plasma_current (int): Current scaling model to use. + 1 = Peng analytic fit + 2 = Peng divertor scaling (TART,STAR) + 3 = Simple ITER scaling + 4 = IPDG89 scaling + 5 = Todd empirical scaling I + 6 = Todd empirical scaling II + 7 = Connor-Hastie model + 8 = Sauter scaling (allowing negative triangularity) + 9 = FIESTA ST scaling + iprofile (int): Switch for current profile consistency. + 0: Use input values for alphaj, rli, dnbeta. + 1: Make these consistent with input q, q_0 values. + 2: Use input values for alphaj, rli. Scale dnbeta with aspect ratio (original scaling). + 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. + kappa (float): Plasma elongation. + kappa95 (float): Plasma elongation at 95% surface. + p0 (float): Central plasma pressure (Pa). + pperim (float): Plasma perimeter length (m). + q0 (float): Plasma safety factor on axis. + q95 (float): Plasma safety factor at 95% flux (= q-bar for i_plasma_current=2). + rli (float): Plasma normalised internal inductance. + rmajor (float): Major radius (m). + rminor (float): Minor radius (m). + sf (float): Shape factor for i_plasma_current=1 (=A/pi in documentation). + triang (float): Plasma triangularity. + triang95 (float): Plasma triangularity at 95% surface. + + Returns: + Tuple[float, float, float, float, float]: Tuple containing bp, qstar, plasma_current, alphaj, rli. + + Raises: + ValueError: If invalid value for i_plasma_current is provided. + + Notes: + This routine calculates the plasma current based on the edge safety factor q95. + It will also make the current profile parameters consistent with the q-profile if required. + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729–1738. https://doi.org/10.13182/FST92-A29971 + - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - Part 1: Physics + - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO + - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant + - Sauter, Geometric formulas for systems codes..., FED 2016 """ # Aspect ratio + aspect_ratio = 1.0 / eps - asp = 1.0 / eps + # Only the Sauter scaling (i_plasma_current=8) is suitable for negative triangularity: + if i_plasma_current != 8 and triang < 0.0: + raise ValueError( + f"Triangularity is negative without i_plasma_current = 8 selected: {triang=}, {i_plasma_current=}" + ) # Calculate the function Fq that scales the edge q from the # circular cross-section cylindrical case - # Only the Sauter scaling (icurr=8) is suitable for negative triangularity: + # Peng analytical fit + if i_plasma_current == 1: + fq = calculate_current_coefficient_peng(eps, sf) - if icurr != 8 and triang < 0.0: - raise ValueError( - f"Triangularity is negative without icurr = 8: {triang=}, {icurr=}" + # Peng scaling for double null divertor; TARTs [STAR Code] + elif i_plasma_current == 2: + plasma_current = 1.0e6 * calculate_plasma_current_peng( + q95, aspect_ratio, eps, rminor, bt, kappa, triang ) - if icurr == 1: # Peng analytical fit - fq = (1.22 - 0.68 * eps) / ((1.0 - eps * eps) ** 2) * sf**2 - elif icurr == 2: # Peng scaling for double null divertor; TARTs [STAR Code] - plascur = 1.0e6 * plasc(q95, asp, eps, rminor, bt, kappa, triang) - elif icurr == 3: # Simple ITER scaling (simply the cylindrical case) + # Simple ITER scaling (simply the cylindrical case) + elif i_plasma_current == 3: fq = 1.0 - elif icurr == 4: # ITER formula (IPDG89) - fq = ( - 0.5 - * (1.17 - 0.65 * eps) - / ((1.0 - eps * eps) ** 2) - * ( - 1.0 - + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3) - ) - ) - elif icurr in [5, 6]: # Todd empirical scalings - fq = ( - (1.0 + 2.0 * eps * eps) - * 0.5 - * (1.0 + kappa95**2) - * ( - 1.24 - - 0.54 * kappa95 - + 0.3 * (kappa95**2 + triang95**2) - + 0.125 * triang95 - ) - ) - fq *= 1 if icurr == 7 else (1.0 + (abs(kappa95 - 1.2)) ** 3) - elif icurr == 7: # Connor-Hastie asymptotically-correct expression - # N.B. If iprofile=1, alphaj will be wrong during the first call (only) - fq = conhas(alphaj, alphap, bt, triang95, eps, kappa95, p0, constants.rmu0) - elif ( - icurr == 8 - ): # Sauter scaling allowing negative triangularity [FED May 2016] - # Assumes zero squareness, note takes kappa, delta at separatrix not _95 + # ITER formula (IPDG89) + elif i_plasma_current == 4: + fq = calculate_current_coefficient_ipdg89(eps, kappa95, triang95) - w07 = 1.0 # zero squareness - can be modified later if required + # Todd empirical scalings + # D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + elif i_plasma_current in [5, 6]: + fq = calculate_current_coefficient_todd(eps, kappa95, triang95, model=1) - fq = ( - (4.1e6 / 5.0e6) - * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) - * (1.0 + 0.09 * triang + 0.16 * triang**2) - * (1.0 + 0.45 * triang * eps) - / (1.0 - 0.74 * eps) - * (1.0 + 0.55 * (w07 - 1.0)) + if i_plasma_current == 6: + fq = calculate_current_coefficient_todd(eps, kappa95, triang95, model=2) + + # Connor-Hastie asymptotically-correct expression + elif i_plasma_current == 7: + # N.B. If iprofile=1, alphaj will be wrong during the first call (only) + fq = calculate_current_coefficient_hastie( + alphaj, alphap, bt, triang95, eps, kappa95, p0, constants.rmu0 ) - elif icurr == 9: - fq = 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060 + + # Sauter scaling allowing negative triangularity [FED May 2016] + # https://doi.org/10.1016/j.fusengdes.2016.04.033. + elif i_plasma_current == 8: + # Assumes zero squareness, note takes kappa, delta at separatrix not _95 + fq = calculate_current_coefficient_sauter(eps, kappa, triang) + + # FIESTA ST scaling + # https://doi.org/10.1016/j.fusengdes.2020.111530. + elif i_plasma_current == 9: + fq = calculate_current_coefficient_fiesta(eps, kappa, triang) else: - raise ValueError(f"Invalid value {icurr=}") + raise ValueError(f"Invalid value {i_plasma_current=}") - if icurr != 2: - plascur = 5.0e6 * rminor**2 / (rmajor * q95) * fq * bt - # == 2 case covered above + # Main plasma current calculation using the fq value from the different settings + if i_plasma_current != 2: + plasma_current = ( + (constants.twopi / constants.rmu0) + * rminor**2 + / (rmajor * q95) + * fq + * bt + ) + # i_plasma_current == 2 case covered above + # Calculate cyclindrical safety factor from IPDG89 qstar = ( 5.0e6 * rminor**2 - / (rmajor * plascur / bt) + / (rmajor * plasma_current / bt) * 0.5 - * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + * (1.0 + kappa**2 * (1.0 + 2.0 * triang**2 - 1.2 * triang**3)) ) physics_variables.normalised_total_beta = ( - 1.0e8 * physics_variables.beta * rminor * bt / plascur - ) - bp = bpol( - icurr, plascur, q95, asp, eps, bt, kappa, triang, pperim, constants.rmu0 + 1.0e8 * physics_variables.beta * rminor * bt / plasma_current + ) + + # Calculate the poloidal field generated by the plasma current + bp = calculate_poloidal_field( + i_plasma_current, + plasma_current, + q95, + aspect_ratio, + eps, + bt, + kappa, + triang, + pperim, + constants.rmu0, ) if iprofile == 1: # Ensure current profile consistency, if required - # This is as described in Hartmann and Zohm only if icurr = 4 as well... + # This is as described in Hartmann and Zohm only if i_plasma_current = 4 as well... # Tokamaks 4th Edition, Wesson, page 116 alphaj = qstar / q0 - 1.0 @@ -2193,7 +2943,7 @@ def culcur( # Menard et al. (2016), Nuclear Fusion, 56, 106023 rli = 3.4 - kappa - return alphaj, rli, bp, qstar, plascur + return alphaj, rli, bp, qstar, plasma_current def outtim(self): po.oheadr(self.outfile, "Times") @@ -2260,7 +3010,7 @@ def outplas(self): * physics_variables.kappa / ( physics_module.total_plasma_internal_energy**2 - * physics_variables.plascur + * physics_variables.plasma_current ) ) @@ -2513,15 +3263,15 @@ def outplas(self): po.ovarin( self.outfile, "Plasma current scaling law used", - "(icurr)", - physics_variables.icurr, + "(i_plasma_current)", + physics_variables.i_plasma_current, ) po.ovarrf( self.outfile, "Plasma current (MA)", - "(plascur/1D6)", - physics_variables.plascur / 1.0e6, + "(plasma_current_MA)", + physics_variables.plasma_current / 1.0e6, "OP ", ) @@ -2583,7 +3333,7 @@ def outplas(self): self.outfile, "Safety factor on axis", "(q0)", physics_variables.q0 ) - if physics_variables.icurr == 2: + if physics_variables.i_plasma_current == 2: po.ovarrf( self.outfile, "Mean edge safety factor", "(q)", physics_variables.q ) @@ -2712,7 +3462,7 @@ def outplas(self): * betath * physics_variables.rminor * physics_variables.bt - / physics_variables.plascur, + / physics_variables.plasma_current, "OP ", ) @@ -4118,35 +4868,35 @@ def outplas(self): if physics_module.err243 == 1: error_handling.report_error(243) - if current_drive_variables.bscfmax < 0.0e0: + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: po.ocmmnt( self.outfile, " (User-specified bootstrap current fraction used)" ) - elif physics_variables.ibss == 1: + elif physics_variables.i_bootstrap_current == 1: po.ocmmnt( self.outfile, " (ITER 1989 bootstrap current fraction model used)" ) - elif physics_variables.ibss == 2: + elif physics_variables.i_bootstrap_current == 2: po.ocmmnt( self.outfile, " (Nevins et al bootstrap current fraction model used)", ) - elif physics_variables.ibss == 3: + elif physics_variables.i_bootstrap_current == 3: po.ocmmnt( self.outfile, " (Wilson bootstrap current fraction model used)" ) - elif physics_variables.ibss == 4: + elif physics_variables.i_bootstrap_current == 4: po.ocmmnt( self.outfile, " (Sauter et al bootstrap current fraction model used)", ) - elif physics_variables.ibss == 5: + elif physics_variables.i_bootstrap_current == 5: po.ocmmnt( self.outfile, " (Sakai et al bootstrap current fraction model used)", ) - if physics_variables.idia == 0: + if physics_variables.i_diamagnetic_current == 0: po.ocmmnt( self.outfile, " (Diamagnetic current fraction not calculated)" ) @@ -4154,20 +4904,20 @@ def outplas(self): if current_drive_variables.diacf_scene > 0.01e0: error_handling.report_error(244) - elif physics_variables.idia == 1: + elif physics_variables.i_diamagnetic_current == 1: po.ocmmnt( self.outfile, " (Hender diamagnetic current fraction scaling used)" ) - elif physics_variables.idia == 2: + elif physics_variables.i_diamagnetic_current == 2: po.ocmmnt( self.outfile, " (SCENE diamagnetic current fraction scaling used)" ) - if physics_variables.ips == 0: + if physics_variables.i_pfirsch_schluter_current == 0: po.ocmmnt( self.outfile, " Pfirsch-Schluter current fraction not calculated" ) - elif physics_variables.ips == 1: + elif physics_variables.i_pfirsch_schluter_current == 1: po.ocmmnt( self.outfile, " (SCENE Pfirsch-Schluter current fraction scaling used)", @@ -4176,22 +4926,22 @@ def outplas(self): po.ovarrf( self.outfile, "Bootstrap fraction (enforced)", - "(bootipf.)", - current_drive_variables.bootipf, + "(bootstrap_current_fraction.)", + current_drive_variables.bootstrap_current_fraction, "OP ", ) po.ovarrf( self.outfile, "Diamagnetic fraction (enforced)", - "(diaipf.)", - current_drive_variables.diaipf, + "(diamagnetic_current_fraction.)", + current_drive_variables.diamagnetic_current_fraction, "OP ", ) po.ovarrf( self.outfile, "Pfirsch-Schlueter fraction (enforced)", - "(psipf.)", - current_drive_variables.psipf, + "(ps_current_fraction.)", + current_drive_variables.ps_current_fraction, "OP ", ) @@ -4199,9 +4949,9 @@ def outplas(self): self.outfile, "Loop voltage during burn (V)", "(vburn)", - physics_variables.plascur + physics_variables.plasma_current * physics_variables.rplas - * physics_variables.facoh, + * physics_variables.inductive_current_fraction, "OP ", ) po.ovarre( @@ -4330,7 +5080,7 @@ def igmarcal(self): physics_variables.kappa95, physics_variables.pchargemw, current_drive_variables.pinjmw, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.pcoreradpv, physics_variables.rmajor, physics_variables.rminor, @@ -4354,71 +5104,108 @@ def igmarcal(self): @staticmethod def bootstrap_fraction_iter89( - aspect, beta, bt, cboot, plascur, q95, q0, rmajor, vol - ): - """Original ITER calculation of bootstrap-driven fraction - of the plasma current. - author: P J Knight, CCFE, Culham Science Centre - aspect : input real : plasma aspect ratio - beta : input real : plasma total beta - bt : input real : toroidal field on axis (T) - cboot : input real : bootstrap current fraction multiplier - plascur : input real : plasma current (A) - q95 : input real : safety factor at 95% surface - q0 : input real : central safety factor - rmajor : input real : plasma major radius (m) - vol : input real : plasma volume (m3) - This routine performs the original ITER calculation of the - plasma current bootstrap fraction. - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + aspect: float, + beta: float, + bt: float, + plasma_current: float, + q95: float, + q0: float, + rmajor: float, + vol: float, + ) -> float: """ - xbs = min(10, q95 / q0) - cbs = cboot * (1.32 - 0.235 * xbs + 0.0185 * xbs**2) - bpbs = ( - constants.rmu0 - * plascur - / (2 * np.pi * np.sqrt(vol / (2 * np.pi**2 * rmajor))) - ) - betapbs = beta * bt**2 / bpbs**2 + Calculate the bootstrap-driven fraction of the plasma current. + + Args: + aspect (float): Plasma aspect ratio. + beta (float): Plasma total beta. + bt (float): Toroidal field on axis (T). + plasma_current (float): Plasma current (A). + q95 (float): Safety factor at 95% surface. + q0 (float): Central safety factor. + rmajor (float): Plasma major radius (m). + vol (float): Plasma volume (m3). + + Returns: + float: The bootstrap-driven fraction of the plasma current. + + This function performs the original ITER calculation of the plasma current bootstrap fraction. + + Reference: + - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + """ + + # Calculate the bootstrap current coefficient + c_bs = 1.32 - 0.235 * (q95 / q0) + 0.0185 * (q95 / q0) ** 2 + + # Calculate the average minor radius + average_a = np.sqrt(vol / (2 * np.pi**2 * rmajor)) + + b_pa = (plasma_current / 1e6) / (5 * average_a) - if betapbs <= 0.0: # only possible if beta <= 0.0 + # Calculate the poloidal beta for bootstrap current + betapbs = beta * (bt / b_pa) ** 2 + + # Ensure betapbs is positive + if betapbs <= 0.0: return 0.0 - return cbs * (betapbs / np.sqrt(aspect)) ** 1.3 + + # Calculate and return the bootstrap current fraction + return c_bs * (betapbs / np.sqrt(aspect)) ** 1.3 @staticmethod def bootstrap_fraction_wilson( - alphaj, alphap, alphat, betpth, q0, qpsi, rmajor, rminor - ): - """Bootstrap current fraction from Wilson et al scaling - author: P J Knight, CCFE, Culham Science Centre - alphaj : input real : current profile index - alphap : input real : pressure profile index - alphat : input real : temperature profile index - beta : input real : total beta - betpth : input real : thermal component of poloidal beta - q0 : input real : safety factor on axis - qpsi : input real : edge safety factor - rmajor : input real : major radius (m) - rminor : input real : minor radius (m) - This function calculates the bootstrap current fraction - using the numerically fitted algorithm written by Howard Wilson. - AEA FUS 172: Physics Assessment for the European Reactor Study - H. R. Wilson, Nuclear Fusion 32 (1992) 257 + alphaj: float, + alphap: float, + alphat: float, + betpth: float, + q0: float, + q95: float, + rmajor: float, + rminor: float, + ) -> float: + """ + Bootstrap current fraction from Wilson et al scaling + + Args: + alphaj (float): Current profile index. + alphap (float): Pressure profile index. + alphat (float): Temperature profile index. + betpth (float): Thermal component of poloidal beta. + q0 (float): Safety factor on axis. + q95 (float): Edge safety factor. + rmajor (float): Major radius (m). + rminor (float): Minor radius (m). + + Returns: + float: The bootstrap current fraction. + + This function calculates the bootstrap current fraction using the numerically fitted algorithm written by Howard Wilson. + + Reference: + - AEA FUS 172: Physics Assessment for the European Reactor Study, 1989 + - H. R. Wilson, Nuclear Fusion 32 (1992) 257 """ + term1 = np.log(0.5) - term2 = np.log(q0 / qpsi) + term2 = np.log(q0 / q95) + + # Re-arranging of parabolic profile to be equal to (r/a)^2 where the profile value is half of the core value termp = 1.0 - 0.5 ** (1.0 / alphap) termt = 1.0 - 0.5 ** (1.0 / alphat) termj = 1.0 - 0.5 ** (1.0 / alphaj) - alfpnw = term1 / np.log(np.log((q0 + (qpsi - q0) * termp) / qpsi) / term2) - alftnw = term1 / np.log(np.log((q0 + (qpsi - q0) * termt) / qpsi) / term2) - aj = term1 / np.log(np.log((q0 + (qpsi - q0) * termj) / qpsi) / term2) + # Assuming a parabolic safety factor profile of the form q = q0 + (q95 - q0) * (r/a)^2 + # Substitute (r/a)^2 term from temperature,pressure and current profiles into q profile when values is 50% of core value + # Take natural log of q profile over q95 and q0 to get the profile index - # Crude check for NaN errors or other illegal values... + alfpnw = term1 / np.log(np.log((q0 + (q95 - q0) * termp) / q95) / term2) + alftnw = term1 / np.log(np.log((q0 + (q95 - q0) * termt) / q95) / term2) + aj = term1 / np.log(np.log((q0 + (q95 - q0) * termj) / q95) / term2) + # Crude check for NaN errors or other illegal values. if np.isnan(aj) or np.isnan(alfpnw) or np.isnan(alftnw) or aj < 0: error_handling.fdiags[0] = aj error_handling.fdiags[1] = alfpnw @@ -4432,12 +5219,14 @@ def bootstrap_fraction_wilson( z = 1.0 # Inverse aspect ratio: r2 = maximum plasma radius, r1 = minimum - + # This is the definition used in the paper r2 = rmajor + rminor r1 = rmajor - rminor eps1 = (r2 - r1) / (r2 + r1) # Coefficients fitted using least squares techniques + + # Square root of current profile index term saj = np.sqrt(aj) a = np.array( @@ -4481,38 +5270,55 @@ def bootstrap_fraction_wilson( @staticmethod def bootstrap_fraction_nevins( - alphan, - alphat, - betat, - bt, - dene, - plascur, - q95, - q0, - rmajor, - rminor, - ten, - zeff, - ): - """Bootstrap current fraction from Nevins et al scaling - author: P J Knight, CCFE, Culham Science Centre - alphan : input real : density profile index - alphat : input real : temperature profile index - betat : input real : total plasma beta (with respect to the toroidal - field) - bt : input real : toroidal field on axis (T) - dene : input real : electron density (/m3) - plascur: input real : plasma current (A) - q0 : input real : central safety factor - q95 : input real : safety factor at 95% surface - rmajor : input real : plasma major radius (m) - rminor : input real : plasma minor radius (m) - ten : input real : density weighted average plasma temperature (keV) - zeff : input real : plasma effective charge - This function calculates the bootstrap current fraction, - using the Nevins et al method, 4/11/90. + alphan: float, + alphat: float, + betat: float, + bt: float, + dene: float, + plasma_current: float, + q95: float, + q0: float, + rmajor: float, + rminor: float, + te: float, + zeff: float, + ) -> float: + """ + Calculate the bootstrap current fraction from Nevins et al scaling. + + Args: + alphan (float): Density profile index. + alphat (float): Temperature profile index. + betat (float): Toroidal plasma beta. + bt (float): Toroidal field on axis (T). + dene (float): Electron density (/m3). + plasma_current (float): Plasma current (A). + q0 (float): Central safety factor. + q95 (float): Safety factor at 95% surface. + rmajor (float): Plasma major radius (m). + rminor (float): Plasma minor radius (m). + te (float): Volume averaged plasma temperature (keV). + zeff (float): Plasma effective charge. + + Returns: + float: The bootstrap current fraction. + + This function calculates the bootstrap current fraction using the Nevins et al method. + + Reference: See appendix of: + Keii Gi, Makoto Nakamura, Kenji Tobita, Yasushi Ono, + Bootstrap current fraction scaling for a tokamak reactor design study, + Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, + ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2014.07.009. + + Nevins, W. M. "Summary report: ITER specialists’ meeting on heating and current drive." + ITER-TN-PH-8-4, June 1988. 1988. + """ - # Calculate peak electron beta + # Calculate peak electron beta at plasma centre, this is not the form used in the paper + # The paper assumes parabolic profiles for calculating core values with the profile indexes. + # We instead use the directly calculated electron density and temperature values at the core. + # So that it is compatible with all profiles betae0 = ( physics_variables.ne0 @@ -4522,13 +5328,13 @@ def bootstrap_fraction_nevins( / (bt**2 / (2.0 * constants.rmu0)) ) - # Call integration routine + # Call integration routine using definite integral routine from scipy ainteg, _ = integrate.quad( - lambda y: bsinteg( + lambda y: _nevins_integral( y, dene, - ten, + te, bt, rminor, rmajor, @@ -4538,80 +5344,112 @@ def bootstrap_fraction_nevins( q0, q95, betat, - constants.echarge, - constants.rmu0, ), - 0, - 0.999, - epsabs=0.001, - epsrel=0.001, + 0, # Lower bound + 1.0, # Upper bound ) # Calculate bootstrap current and fraction aibs = 2.5 * betae0 * rmajor * bt * q95 * ainteg - return 1.0e6 * aibs / plascur + return 1.0e6 * aibs / plasma_current @staticmethod - def bootstrap_fraction_sauter(plasma_profile): - """Bootstrap current fraction from Sauter et al scaling - author: P J Knight, CCFE, Culham Science Centre - None - This function calculates the bootstrap current fraction - using the Sauter, Angioni and Lin-Liu scaling. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 + def bootstrap_fraction_sauter(plasma_profile: float) -> float: + """Calculate the bootstrap current fraction from the Sauter et al scaling. + + Args: + plasma_profile (PlasmaProfile): The plasma profile object. + + Returns: + float: The bootstrap current fraction. + + This function calculates the bootstrap current fraction using the Sauter, Angioni, and Lin-Liu scaling. + + Reference: + - O. Sauter, C. Angioni, Y. R. Lin-Liu; + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. + Phys. Plasmas 1 July 1999; 6 (7): 2834–2839. https://doi.org/10.1063/1.873240 + - O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime + [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 + + Note: + The code was supplied by Emiliano Fable, IPP Garching (private communication). """ + + # Number of radial data points along the profile NR = plasma_profile.profile_size + # Radial points from 0 to 1 seperated by 1/profile_size roa = plasma_profile.neprofile.profile_x + + # Local circularised minor radius rho = np.sqrt(physics_variables.xarea / np.pi) * roa + # Square root of local aspect ratio sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) + # Calculate electron and ion density profiles ne = plasma_profile.neprofile.profile_y * 1e-19 ni = (physics_variables.dnitot / physics_variables.dene) * ne + # Calculate electron and ion temperature profiles tempe = plasma_profile.teprofile.profile_y tempi = (physics_variables.ti / physics_variables.te) * tempe + # Flat Zeff profile assumed - zef = np.full_like(tempi, physics_variables.zeff) + # Return tempi like array object filled with zeff + zeff = np.full_like(tempi, physics_variables.zeff) - # mu = 1/safety factor + # inverse_q = 1/safety factor # Parabolic q profile assumed - mu = 1 / ( + inverse_q = 1 / ( physics_variables.q0 + (physics_variables.q - physics_variables.q0) * roa**2 ) - amain = np.full_like(mu, physics_variables.afuel) - zmain = np.full_like(mu, 1.0 + physics_variables.fhe3) + # Create new array of average mass of fuel portion of ions + amain = np.full_like(inverse_q, physics_variables.afuel) + # Create new array of average main ion charge + zmain = np.full_like(inverse_q, 1.0 + physics_variables.fhe3) + + # Prevent division by zero if ne[NR - 1] == 0.0: ne[NR - 1] = 1e-4 * ne[NR - 2] ni[NR - 1] = 1e-4 * ni[NR - 2] + # Prevent division by zero if tempe[NR - 1] == 0.0: tempe[NR - 1] = 1e-4 * tempe[NR - 2] tempi[NR - 1] = 1e-4 * tempi[NR - 2] # Calculate total bootstrap current (MA) by summing along profiles - # Looping from 2 because dcsa etc should return 0 @ j == 1 - nr_rng = np.arange(2, NR) - nr_rng_1 = nr_rng - 1 - drho = rho[nr_rng] - rho[nr_rng_1] - da = 2 * np.pi * rho[nr_rng_1] * drho # area of annulus - dlogte_drho = (np.log(tempe[nr_rng]) - np.log(tempe[nr_rng_1])) / drho - dlogti_drho = (np.log(tempi[nr_rng]) - np.log(tempi[nr_rng_1])) / drho - dlogne_drho = (np.log(ne[nr_rng]) - np.log(ne[nr_rng_1])) / drho + # Looping from 2 because _calculate_l31_coefficient() etc should return 0 @ j == 1 + radial_elements = np.arange(2, NR) + + # Change in localised minor radius to be used as delta term in derivative + drho = rho[radial_elements] - rho[radial_elements - 1] + + # Area of annulus, assuming circular plasma cross-section + da = 2 * np.pi * rho[radial_elements - 1] * drho # area of annulus + + # Create the partial derivatives + dlogte_drho = ( + np.log(tempe[radial_elements]) - np.log(tempe[radial_elements - 1]) + ) / drho + dlogti_drho = ( + np.log(tempi[radial_elements]) - np.log(tempi[radial_elements - 1]) + ) / drho + dlogne_drho = ( + np.log(ne[radial_elements]) - np.log(ne[radial_elements - 1]) + ) / drho + jboot = ( 0.5 * ( - dcsa( - nr_rng, + _calculate_l31_coefficient( + radial_elements, NR, physics_variables.rmajor, physics_variables.bt, @@ -4620,14 +5458,14 @@ def bootstrap_fraction_sauter(plasma_profile): ni, tempe, tempi, - mu, + inverse_q, rho, - zef, + zeff, sqeps, ) * dlogne_drho - + hcsa( - nr_rng, + + _calculate_l31_32_coefficient( + radial_elements, NR, physics_variables.rmajor, physics_variables.bt, @@ -4636,19 +5474,19 @@ def bootstrap_fraction_sauter(plasma_profile): ni, tempe, tempi, - mu, + inverse_q, rho, - zef, + zeff, sqeps, ) * dlogte_drho - + xcsa( - nr_rng, + + _calculate_l34_alpha_31_coefficient( + radial_elements, NR, physics_variables.rmajor, physics_variables.bt, physics_variables.triang, - mu, + inverse_q, sqeps, tempi, tempe, @@ -4657,7 +5495,7 @@ def bootstrap_fraction_sauter(plasma_profile): ni, ne, rho, - zef, + zeff, ) * dlogti_drho ) @@ -4665,12 +5503,12 @@ def bootstrap_fraction_sauter(plasma_profile): * ( -physics_variables.bt / (0.2 * np.pi * physics_variables.rmajor) - * rho[nr_rng_1] - * mu[nr_rng_1] + * rho[radial_elements - 1] + * inverse_q[radial_elements - 1] ) ) # A/m2 - return np.sum(da * jboot, axis=0) / physics_variables.plascur + return np.sum(da * jboot, axis=0) / physics_variables.plasma_current @staticmethod def bootstrap_fraction_sakai( @@ -4697,7 +5535,7 @@ def bootstrap_fraction_sakai( 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 profile assumed for the alphan and alpat indexes is only a parabolic 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. @@ -4711,11 +5549,11 @@ def bootstrap_fraction_sakai( 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 + # So the diamganetic current should not be calculated with this. i_diamagnetic_current = 0 return ( 10 ** (0.951 * eps - 0.948) * betap ** (1.226 * eps + 1.584) - * rli ** (-0.184 - 0.282) + * rli ** (-0.184 * eps - 0.282) * (q95 / q0) ** (-0.042 * eps - 0.02) * alphan ** (0.13 * eps + 0.05) * alphat ** (0.502 * eps - 0.273) @@ -4768,7 +5606,7 @@ def fhz(self, hhh): physics_variables.kappa95, physics_variables.pchargemw, current_drive_variables.pinjmw, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.pcoreradpv, physics_variables.rmajor, physics_variables.rminor, @@ -4825,7 +5663,7 @@ def pcond( kappa95, pchargemw, pinjmw, - plascur, + plasma_current, pcoreradpv, rmajor, rminor, @@ -4858,7 +5696,7 @@ def pcond( kappaa : output real : plasma elongation calculated using area ratio pchargemw : input real : non-alpha charged particle fusion power (MW) pinjmw : input real : auxiliary power to ions and electrons (MW) - plascur : input real : plasma current (A) + plasma_current : input real : plasma current (A) pcoreradpv: input real : total core radiation power (MW/m3) q : input real : edge safety factor (tokamaks), or rotational transform iotabar (stellarators) @@ -4939,7 +5777,7 @@ def pcond( n20 = dene / 1.0e20 # Plasma current in MA - pcur = plascur / 1.0e6 + pcur = plasma_current / 1.0e6 # Separatrix kappa defined with X-section for general use kappaa = xarea / (np.pi * rminor * rminor) @@ -5570,15 +6408,15 @@ def pcond( elif isc == 42: # High density relevant confinement scaling # P.T. Lang et al. 2012, IAEA conference proceeding EX/P4-01 - # q should be q95: incorrect if icurr = 2 (ST current scaling) + # q should be q95: incorrect if i_plasma_current = 2 (ST current scaling) qratio = q / qstar # Greenwald density in m^-3 - nGW = 1.0e14 * plascur / (np.pi * rminor * rminor) + nGW = 1.0e14 * plasma_current / (np.pi * rminor * rminor) nratio = dnla / nGW tauee = ( hfact * 6.94e-7 - * plascur**1.3678e0 + * plasma_current**1.3678e0 * bt**0.12e0 * dnla**0.032236e0 * (powerht * 1.0e6) ** (-0.74e0) @@ -5595,7 +6433,7 @@ def pcond( tauee = ( hfact * 0.014e0 - * (plascur / 1.0e6) ** 0.68e0 + * (plasma_current / 1.0e6) ** 0.68e0 * bt**0.77e0 * dnla20**0.02e0 * powerht ** (-0.29e0) @@ -5605,7 +6443,7 @@ def pcond( tauee = ( hfact * 0.014e0 - * (plascur / 1.0e6) ** 0.60e0 + * (plasma_current / 1.0e6) ** 0.60e0 * bt**0.70e0 * dnla20 ** (-0.03e0) * powerht ** (-0.33e0) @@ -5615,7 +6453,7 @@ def pcond( tauee = ( hfact * 0.014e0 - * (plascur / 1.0e6) ** 0.76e0 + * (plasma_current / 1.0e6) ** 0.76e0 * bt**0.84e0 * dnla20**0.07 * powerht ** (-0.25e0) diff --git a/process/plasma_geometry.py b/process/plasma_geometry.py index bd606e404e..4bb9bb1077 100644 --- a/process/plasma_geometry.py +++ b/process/plasma_geometry.py @@ -233,8 +233,8 @@ def geomty(self): ) physics_variables.sareao = xso - # icurr = 8 specifies use of the Sauter geometry as well as plasma current. - if physics_variables.icurr == 8: + # i_plasma_current = 8 specifies use of the Sauter geometry as well as plasma current. + if physics_variables.i_plasma_current == 8: ( physics_variables.pperim, physics_variables.sf, diff --git a/process/pulse.py b/process/pulse.py index 561b2ce487..6727c007ac 100755 --- a/process/pulse.py +++ b/process/pulse.py @@ -92,7 +92,7 @@ def tohswg(self, output: bool) -> None: # Maximum rate of change of plasma current (A/s) # - now a function of the plasma current itself (previously just 0.5e6) - ipdot = 0.0455e0 * physics_variables.plascur + ipdot = 0.0455e0 * physics_variables.plasma_current # Minimum plasma current ramp-up time (s) # - corrected (bus resistance is not a function of pfcoil_variables.turns) @@ -149,9 +149,9 @@ def burn(self, output: bool): # Loop voltage during flat-top (including MHD sawtooth enhancement) vburn = ( - physics_variables.plascur + physics_variables.plasma_current * physics_variables.rplas - * physics_variables.facoh + * physics_variables.inductive_current_fraction * physics_variables.csawth ) diff --git a/process/stellarator.py b/process/stellarator.py index 494f0bef77..62ff4a483b 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -221,7 +221,7 @@ def stigma(self): physics_variables.kappa95, physics_variables.pchargemw, current_drive_variables.pinjmw, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.pcoreradpv, physics_variables.rmajor, physics_variables.rminor, @@ -4203,7 +4203,7 @@ def stphys(self, output): physics_variables.kappa95, physics_variables.pchargemw, current_drive_variables.pinjmw, - physics_variables.plascur, + physics_variables.plasma_current, physics_variables.pcoreradpv, physics_variables.rmajor, physics_variables.rminor, @@ -4242,7 +4242,7 @@ def stphys(self, output): physics_variables.deni, physics_variables.fusionrate, physics_variables.alpharate, - physics_variables.plascur, + physics_variables.plasma_current, sbar, physics_variables.dnalp, physics_variables.taueff, diff --git a/process/structure.py b/process/structure.py index f4cac89f40..f8fc547a9f 100644 --- a/process/structure.py +++ b/process/structure.py @@ -47,7 +47,7 @@ def run(self, output: bool = False) -> None: stv.coldmass, stv.gsmass, ) = self.structure( - pv.plascur, + pv.plasma_current, pv.rmajor, pv.rminor, pv.kappa, diff --git a/process/utilities/errorlist.json b/process/utilities/errorlist.json index 20a2453f27..4b168b848e 100644 --- a/process/utilities/errorlist.json +++ b/process/utilities/errorlist.json @@ -193,12 +193,12 @@ { "no": 37, "level": 2, - "message": "CHECK: Usual current scaling for TARTs (icurr=2 or 9) is not being used" + "message": "CHECK: Usual current scaling for TARTs (i_plasma_current=2 or 9) is not being used" }, { "no": 38, "level": 3, - "message": "CHECK: Invalid boostrap current law for ST, do not use ibss = 1" + "message": "CHECK: Invalid boostrap current law for ST, do not use i_bootstrap_current = 1" }, { "no": 39, @@ -208,7 +208,7 @@ { "no": 40, "level": 2, - "message": "CHECK: icurr=2 is not a valid option for a non-TART device" + "message": "CHECK: i_plasma_current=2 is not a valid option for a non-TART device" }, { "no": 41, @@ -383,7 +383,7 @@ { "no": 75, "level": 3, - "message": "PHYSICS: Illegal value of ibss" + "message": "PHYSICS: Illegal value of i_bootstrap_current" }, { "no": 76, @@ -393,7 +393,7 @@ { "no": 77, "level": 3, - "message": "CULCUR: Illegal value for icurr" + "message": "CALCULATE_PLASMA_CURRENT: Illegal value for i_plasma_current" }, { "no": 78, @@ -1118,7 +1118,7 @@ { "no": 222, "level": 3, - "message": "CHECK: Lang 2012 confinement scaling cannot be used for icurr=2 due to wrong q" + "message": "CHECK: Lang 2012 confinement scaling cannot be used for i_plasma_current=2 due to wrong q" }, { "no": 223, @@ -1228,7 +1228,7 @@ { "no": 244, "level": 2, - "message": "PHYSICS: Diamagnetic fraction is more than 1%, but not calculated. Consider using idia=2 and ips=1." + "message": "PHYSICS: Diamagnetic fraction is more than 1%, but not calculated. Consider using i_diamagnetic_current=2 and i_pfirsch_schluter_current=1." }, { "no": 245, @@ -1428,7 +1428,7 @@ { "no": 284, "level": 3, - "message": "CHECK: idia = 0 should be used with the Sakai plasma current scaling." + "message": "CHECK: i_diamagnetic_current = 0 should be used with the Sakai plasma current scaling." }, { "no": 285, diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90 index a15bdd8501..d85a8f3c45 100755 --- a/source/fortran/constraint_equations.f90 +++ b/source/fortran/constraint_equations.f90 @@ -1985,7 +1985,7 @@ subroutine constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) !! Logic change during pre-factoring: err, symbol, units will be assigned only if present. !! fq : input real : f-value for edge safety factor !! q : safety factor 'near' plasma edge: equal to q95 - !! (unless icurr = 2 (ST current scaling), in which case q = mean edge safety factor qbar) + !! (unless i_plasma_current = 2 (ST current scaling), in which case q = mean edge safety factor qbar) !! qlim : input real : lower limit for edge safety factor !! itart : input integer : switch for spherical tokamak (ST) models: