🧲 Update vscalc() function and add to documentation

documentation/proc-pages/development/standards.md

@@ -254,7 +254,7 @@ This should be used for units of $\text{kg} \cdot \text{m}^{-2}\text{s}^{-1}$
 
 ##### Inductances
 
-- Inductances should start with the `h_` prefix
+- Inductances should start with the `ind_` prefix
 
 ---------------------
 
@@ -313,6 +313,15 @@ This should be used for units of $\text{kg} \cdot \text{m}^{-2}\text{s}^{-1}$
 
 ---------------------
 
+##### Magnetic flux
+
+- Magnetic fluxes can start with the `web_` prefix representing Webers.
+
+- Since magnetic flux units are more commonly used in inductive current drive it may be more appropriate
+    to use the `vs_` prefix instead representing a $\text{Vs}$.
+
+---------------------
+
 ##### Frequencies
 
 - Frequencies should start with the `freq_`

documentation/proc-pages/eng-models/central-solenoid.md

@@ -78,7 +78,7 @@ in the case of the CS coils).
     7.  The stress is monotonic (since the hoop stress is always positive), and its minimum value 
         is the residual stress (input).  Default: 240 MPa.  
     8.  The mean stress is taken into account using the [Walker](https://en.wikipedia.org/wiki/Crack_growth_equation#Walker_equation) 
-        modification of the Paris equation, with coefficient $\gamma$=0.436  
+        modification of the Paris equation, with coefficient $\ejima_coeff$=0.436  
     9.  Failure occurs when the crack dimension a equals the conduit thickness, or dimension c reaches 
         the conductor width.  
     10. No safety factor is used for the number of cycles.  

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

@@ -123,15 +123,15 @@ be calculated. The following switch options are available below:
 
 This can be activated by stating `iprofile = 0` in the input file.
 
-`alphaj`, `rli` and `beta_norm_max` are inputs.
+`alphaj`, `ind_plasma_internal_norm` and `beta_norm_max` are inputs.
 
 ---------
 
 #### Wesson Relation
 
 This can be activated by stating `iprofile = 1` in the input file.
 
-`alphaj`, `rli` and `beta_norm_max` are calculated consistently. 
+`alphaj`, `ind_plasma_internal_norm` and `beta_norm_max` are calculated consistently. 
 
 `beta_norm_max` is calculated using:  
 
@@ -149,7 +149,7 @@ This is only recommended for high aspect ratio tokamaks[^3].
 
 This can be activated by stating `iprofile = 2` in the input file.
 
-`alphaj` and `rli` are inputs. `beta_norm_max` calculated using:
+`alphaj` and `ind_plasma_internal_norm` are inputs. `beta_norm_max` calculated using:
 
 $$
 g=2.7(1+5\epsilon^{3.5})
@@ -224,7 +224,7 @@ $$
 
 This can be activated by stating `iprofile = 3` in the input file.
 
-`alphaj` and `rli` are inputs. `beta_norm_max` calculated using[^4]:
+`alphaj` and `ind_plasma_internal_norm` are inputs. `beta_norm_max` calculated using[^4]:
 
 $$
 g=3.12+3.5\epsilon^{1.7}
@@ -299,23 +299,23 @@ $$
 
 This can be activated by stating `iprofile = 4` in the input file.
 
-`alphaj` and `beta_norm_max` are inputs. `rli` calculated from elongation [^4]. This is only recommended for spherical tokamaks.
+`alphaj` and `beta_norm_max` are inputs. `ind_plasma_internal_norm` calculated from elongation [^4]. This is only recommended for spherical tokamaks.
 
 ---------
 
 #### Menard Beta & Inductance Relation
 
 This can be activated by stating `iprofile = 5` in the input file.
 
-`alphaj` is an input.  `rli` calculated from elongation and `beta_norm_max` calculated using $g=3.12+3.5\epsilon^{1.7}$ [^4]. This is only recommended for spherical tokamaks.
+`alphaj` is an input.  `ind_plasma_internal_norm` calculated from elongation and `beta_norm_max` calculated using $g=3.12+3.5\epsilon^{1.7}$ [^4]. This is only recommended for spherical tokamaks.
 
 ---------
 
 #### Tholerus Relation
 
 This can be activated by stating `iprofile = 6` in the input file.
 
-`alphaj` and `c_beta` are inputs.  `rli` calculated from elongation and `beta_norm_max` calculated using 
+`alphaj` and `c_beta` are inputs.  `ind_plasma_internal_norm` calculated from elongation and `beta_norm_max` calculated using 
 
 $$
 C_{\beta}\approx\frac{(g-3.7)F_p}{12.5-3.5 F_p}
@@ -327,7 +327,7 @@ where $F_p$ is the pressure peaking, $F_p = p_{\text{ax}} / \langle p \rangle$ a
 
 ---------
 
-Further details on the calculation of `alphaj` and `rli` is given in [Plasma Current](./plasma_current.md).
+Further details on the calculation of `alphaj` and `ind_plasma_internal_norm` is given in [Plasma Current](./plasma_current.md).
 
 ----------------------
 

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

@@ -1,7 +1,163 @@
 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.
+This calculated fraction (`inductive_current_fraction`) is then used in the `calculate_volt_second_requirements()` 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.
+    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.
+
+
+--------------------
+
+## Volt-second requirements | `calculate_volt_second_requirements()`
+
+The plasma requires a constants magnetic flux change in order to keep inductively driving a current in itself.
+
+By Faraday's law of induction, any change in flux through a circuit induces an electromotive force (EMF, $V$),
+in the circuit, proportional to the rate of change of flux.
+
+$$
+V(t) = - \frac{\mathrm{d}}{\mathrm{d}t} \Phi(t)
+$$
+
+Inductance, $L$ is the ratio between the induced voltage from the flux change and the rate of change of current that produced the flux change.
+
+$$
+V(t) = L \ \frac{\mathrm{d}I}{\mathrm{d}t} 
+$$
+
+
+The flux requirements are defined by the sum of the pulse ramp up and flat top / burn phases.
+
+-----------
+
+### Current ramp phase
+
+#### Resistive Component
+
+In the ramp up phase we need to take the plasma current from 0 Amps to the plasma current value, which is normally ten's of Mega Amps. Since the plasma has a non zero resistance, the current induced in the plasma will be dissipated due to resistive losses. This can be tricky to calculate as the plasma resistance decreases as the plasma temperature increases towards our flat-top full plasma scenario state. The inductance of the plasma also varies during this ramp phase and so the amount of current driven for the same change in flux will vary too. 
+
+Thankfully a formulation of the resistive flux consumption during ramp phase is given by Ejima et.al [^1].
+
+$$
+\mathtt{vs\_res\_ramp} = \Phi_{\text{res,ramp}} = C_{\text{eji}}\mu_0I_{\text{p}}R
+$$
+
+where $C_{\text{ejima}}$ is the empirical Ejima coefficient defined by the user or at its default of 0.4.
+
+This relation is is done by analyzing a wide range of cross-sections,
+ranging from circular to doublet produced in the Doublet III machine. The plasma cross-section is
+slightly elongated, with $\kappa=1.2$. The toroidal field is 2.4 $\text{T}$. The current swing of the Ohmic-heating transformer is nominally from $-25 \text{kA}$ to $+80 \text{kA}$, corresponding to a flux swing of $\approx 2.6 \ \text{Vs}$.
+
+The initial one-turn loop voltage is around $50 \ \text{V}$, causing the plasma current to rise to $300 \ \text{kA}$ in about $80 \  \text{ms}$. The plasma current is then increased to higher flat-top currents at a steady rate of $2 \text{MA} \text{s}^{-1}$.
+
+The calculation of $\Phi_{\text{res,ramp}}$ is based on the fact that the ramp up takes of the order of the resistive current penetration time $\frac{\mu_0 a^2}{\rho_{\text{p}}}$, where $\rho_{\text{p}}$ is the resistivity of the plasma. The
+flux consumption is therefore independent of the resistivity and the minor radius.
+
+-------------
+
+#### Self-Inductance Component
+
+
+The internal inductance is defined as the part of the inductance obtained by integrating over the plasma volume.
+
+The internal component of the plasma self inductance flux consumption during the ramp up phase is given by:
+
+$$
+\mathtt{vs\_plasma\_internal} = \Phi_{\text{ind,internal,ramp}} =  I_{\text{p}} \times \underbrace{\left[\frac{\mu_0 R l_i}{2}\right]}_{\mathtt{ind\_plasma\_internal}}
+$$
+
+Here $l_i$ is known as the normalized internal inductance defined for circular cross section plasmas with minor radius $a$:
+
+$$
+l_i = \frac{\langle B_{\text{p}}^2 \rangle }{B_{\text{p}}^2(a)}
+$$
+
+You may also see the following approximation term used, $l_i(3)$[^2].
+
+$$
+l_i(3) = \frac{2V\langle B_{\text{p}}^2 \rangle }{\mu_0^2I^2R}
+$$
+
+which is equal to $l_i$ if the plasma has a perfect circular cross-section.
+
+The external inductance needs to be accounted for also as even though we assume the toroidal current density vanishes at the plasma edge, there still exists a vacuum poloidal magnetic field around the plasma.
+
+Hirshman et.al[^3] gives a formula for the external plasma inductance in the form:
+
+$$
+\mathtt{ind\_plasma\_external} = L_{\text{ext}} = \mu_0 R\frac{a(\epsilon)(1-\epsilon)}{1-\epsilon + b(\epsilon)\kappa}
+$$
+
+$$
+a(\epsilon)  = (1+1.18\sqrt{\epsilon}+2.05\epsilon)\ln{\left(\frac{8}{\epsilon}\right)} \\
+- (2.0+9.25\sqrt{\epsilon}-1.21\epsilon)
+$$
+
+$$
+b(\epsilon) = 0.73\sqrt{\epsilon}(1+2\epsilon^4 - 6\epsilon^5+3.7\epsilon^6)
+$$
+
+where $\epsilon$ is the plasma inverse aspect ratio and $\kappa$ is the separatrix elongation.
+
+
+The total plasma inductance is then calculated as:
+
+$$
+\mathtt{ind\_plasma\_total} = \mathtt{ind\_plasma\_external} + \mathtt{ind\_plasma\_internal}
+$$
+
+Therefore the total inductive flux consumption during ramp up is given by:
+
+$$
+\mathtt{vs\_self\_ind\_ramp} = \Phi_{\text{ind,ramp}} =  \mathtt{ind\_plasma\_total} \times I_{\text{p}}
+$$
+
+So the total resisitive and inductive flux consumption at current ramp up which is the total flux requires is given by:
+
+$$
+\mathtt{vs\_ramp\_required} = \mathtt{vs\_res\_ramp} + \mathtt{vs\_self\_ind\_ramp} \\
+\Phi_{\text{tot,ramp}} =  \Phi_{\text{res,ramp}} + \Phi_{\text{ind,ramp}}
+$$
+
+------------------
+
+### Steady current burn phase
+
+At plasma current flat-top there is no self inductance contribution as the plasma current is expected to be at a constant value. However there is still a resistive contribution.
+
+
+
+#### Resistive Component
+
+
+For the flat top resistive component we can just take the loop voltage value based on the plasmas resistivity and the known fraction of current driven inductively:
+
+$$
+\mathtt{v\_burn\_resistive} = V_{\text{loop}} = I_{\text{p}}  \rho_\text{p} f_{\text{ind}}
+$$
+
+where $\rho_\text{p}$ is the calculated [plasma resistivity](./plasma_resistive_heating.md) and $f_{\text{ind}}$ is the inductive current fraction.
+
+The total flux required is then simply found by multiplying the loop voltage above by the required duration of the burn phase:
+
+$$
+\mathtt{vs\_burn\_required} = \Phi_{\text{res,burn}} = I_{\text{p}}  \rho_\text{p} f_{\text{ind}} \times T_{\text{burn}} \\
+= \mathtt{v\_burn\_resistive} \times \mathtt{t\_burn}
+$$
+
+----------------
+
+Finally we can now find the minimum flux required for the full duration of the pulse:
+
+$$
+\mathtt{vs\_total\_required} = \Phi_{\text{tot}} = \Phi_{\text{tot,ramp}} + \Phi_{\text{res,burn}} \\
+= \mathtt{vs\_ramp\_required} + \mathtt{vs\_burn\_required}
+$$
+
+
+
+[^1]: S. Ejima, R. W. Callis, J. L. Luxon, R. D. Stambaugh, T. S. Taylor, and J. C. Wesley, “Volt-second analysis and consumption in Doublet III plasmas,” Nuclear Fusion, vol. 22, no. 10, pp. 1313-1319, Oct. 1982, doi: https://doi.org/10.1088/0029-5515/22/10/006.
+[^2]: G. L. Jackson et al., “ITER startup studies in the DIII-D tokamak,” Nuclear Fusion, vol. 48, no. 12, p. 125002, Nov. 2008, doi: https://doi.org/10.1088/0029-5515/48/12/125002.
+[^3]: S. P. Hirshman and G. H. Neilson, “External inductance of an axisymmetric plasma,” The Physics of Fluids, vol. 29, no. 3, pp. 790-793, Mar. 1986, doi: https://doi.org/10.1063/1.865934.
+‌

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

@@ -92,7 +92,7 @@ Instead of $q_a$, $q_{95}$ is used as in plasma configurations with divertors th
 
 ## 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.
+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{ind_plasma_internal_norm}$ ($l_{\mathtt{i}}$). This is done in 5 separate steps which are shown in the following numbered sections.
 
 
 $$\begin{aligned}
@@ -556,7 +556,7 @@ $$
 
 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 (`beta_norm_max`) using the 
+profile peaking factor $\alpha_J$ (`alphaj`),  the normalised internal inductance $l_i$ (`ind_plasma_internal_norm`) and beta limit $g$-factor (`beta_norm_max`) using the 
 safety factor on axis `q0` and the cylindrical safety factor $q_*$ (`qstar`):   
 
 $$\begin{aligned}

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

@@ -2,24 +2,47 @@
 
 The ohmic component of the plasma heating is given by that from the ITER 1989 Physics Design Guidelines[^1]
 
+-------------------
+
+## Plasma Ohmic Heating Power | `plasma_ohmic_heating()`
+
 Using the resistive loop voltage for a reference profile of parabolic shape with:
 
 $$
 \alpha_n \approx 0.5, \alpha_T \approx 1.0,  \alpha_J \approx 1.5
 $$
 
+We calculate the plasma resistance as:
 
 $$
-\Omega_{\text{plasma}} \approx 2.15 \times 10^{-3} Z_{\text{eff}}\langle \gamma_{\text{NC}} \rangle \frac{R_0}{\kappa a^2} \frac{1}{T_{10}^{1.5}}
+\mathtt{res\_plasma} = \Omega_{\text{plasma}} \approx 2.15 \times 10^{-3} Z_{\text{eff}}\langle \gamma_{\text{NC}} \rangle \frac{R_0}{\kappa a^2} \frac{1}{T_{10}^{1.5}}
 $$
 
-The neoclassical (avergae) resisitivity enhancement factor $\left(\langle \gamma_{\text{NC}} \rangle \right)$ is given by an empirical fit:
+where $Z_{\text{eff}}$ is the plasma effective charge and $T_{10}$ is the density-weighted temperature in units of 10 keV.
+
+The neoclassical (average) resistivity enhancement factor $\left(\langle \gamma_{\text{NC}} \rangle \right)$ is given by an empirical fit:
 
 $$
 \langle \gamma_{\text{NC}} \rangle = 4.3 -0.6A
 $$
 
-where $A$ is valid in the range of 2.5 - 4.0.
+where $A$ is valid in the range of 2.5 - 4.0. If $A < 2.5$ then $\langle \gamma_{\text{NC}}\rangle$ is et equal to 1.0
+
+---------------
+
+The ohmic heating power in MW is then simply found using Joules law:
+
+$$
+\mathtt{p\_plasma\_ohmic\_mw} = P_{\text{OH}} = 1\times 10^{-6} \left[f_{\text{ind}}I_{\text{p}}^2\Omega_{\text{plasma}}\right]
+$$
+
+where $f_{\text{ind}}$ is the fraction of plasma current driven by inductive means.
+
+Likewise, the ohmic heating per unit volume simply as:
+
+$$
+\mathtt{pden\_plasma\_ohmic\_mw} = \frac{\mathtt{p\_plasma\_ohmic\_mw}}{V_{\text{p}}}
+$$
 
 
 [^1]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989',

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

@@ -193,7 +193,7 @@ will be inaccurate for single-null plasmas, `i_single_null = 1`)
 - `i_plasma_geometry = 8` -- The input values for `kappa` and `triang` are used directly and the 95% flux 
   surface values are calculated using the FIESTA fit from `i_plasma_geometry = 7`.
 ---------------------------------------------------------------------
-- `i_plasma_geometry = 9` -- The input values for `triang` and `rli` are used, `kappa` and the 95% flux 
+- `i_plasma_geometry = 9` -- The input values for `triang` and `ind_plasma_internal_norm` are used, `kappa` and the 95% flux 
   surface values are calculated.
 
   $$

documentation/proc-pages/physics-models/pulsed-plant.md

@@ -5,8 +5,8 @@ necessary to operate the plant in a pulsed manner as the current swing in the ce
 cannot be continued indefinitely. `PROCESS` can perform a number of calculations relevant to a 
 pulsed power plant, as detailed below.
 
-Switch `lpulse` determines whether the power plant is assumed to be based on steady-state 
-(`lpulse = 0`) or pulsed (`lpulse = 1`) operation.  The current ramp calculations apply in both 
+Switch `i_pulsed_plant` determines whether the power plant is assumed to be based on steady-state 
+(`i_pulsed_plant = 0`) or pulsed (`i_pulsed_plant = 1`) operation.  The current ramp calculations apply in both 
 cases, as even a steady-state reactor has to be started up.
 
 ## Start-up power requirements
@@ -29,14 +29,14 @@ change equal to the maximum proposed in [^1], or it can be set by the user.  The
 constraint is likely to depend on whether the ramp-up is purely inductive or includes current drive, 
 but this is not taken ito account.
 
-In the steady-state scenario (`lpulse` = 0), the plasma current ramp-up time `t_current_ramp_up` is determined as follows. 
+In the steady-state scenario (`i_pulsed_plant` = 0), the plasma current ramp-up time `t_current_ramp_up` is determined as follows. 
 
 - If `tohsin` = 0, the rate of change of plasma current is 0.5 MA/s. The PF coil ramp time `t_precharge` 
   and shutdown time `t_ramp_down` are (arbitrarily) set equal to `t_current_ramp_up`. 
 - If `tohsin` $\neq$ 0, the plasma current ramp-up time `t_current_ramp_up` = `tohsin`, and the PF coil ramp 
   and shutdown times are input parameters.
 
-In the pulsed scenario, (`lpulse` = 1), the plasma current ramp-up time `t_current_ramp_up` is an input, and it 
+In the pulsed scenario, (`i_pulsed_plant` = 1), the plasma current ramp-up time `t_current_ramp_up` is an input, and it 
 can be set as an iteration variable (65). The ramp-up and shutdown time in the pulsed case are set 
 equal to `t_current_ramp_up`. To ensure that the plasma current ramp rate during start-up is prevented from being 
 too high, as governed by the requirement to maintain plasma stability by ensuring that the induced