Add proper ion & electron equilibration

documentation/source/physics-models/detailed_physics.md

@@ -160,3 +160,13 @@ For $T_\text{e}$ in eV
 $$
 \eta_{\text{spitzer}} = \frac{4 \sqrt{2 \pi}}{3 }\frac{Z_i e^2 \sqrt{m_{\text{e}}} \ln \Lambda_{\text{ei}}}{\left(4 \pi \epsilon_0 \right)^2 T_{\text{e}}^{\frac{3}{2}}}
 $$
+
+------------------
+
+### Ion-electron equilibration time | `calculate_equilibriation_time()`
+
+For species temperatures $T$, measured in $[\text{eV}]$
+
+$$
+\frac{3(4\pi\epsilon_0)^2m_{\text{i}} m_{\text{e}}}{8\sqrt{2\pi}n_{\text{i}} \text{e}^4 Z_{\text{i}}^2 \ln\Lambda_{\text{ei}}}\left(\frac{T_{\text{e}}}{m_{\text{e}}}+\frac{T_{\text{i}}}{m_{\text{i}}}\right)^\frac{3}{2}
+$$

documentation/source/physics-models/plasma_confinement.md

@@ -647,60 +647,105 @@ $$
 
 -------------------------
 
+### Species Confinement Times
+
+Once the global energy confinement time $\tau_{\text{E}}$ is calculated from the scaling we calculate the species and $H$ factor corrected confinement time.
+
+Using out $H$ factor we find the new corrected value for the global energy confinement time:
+
+$$
+\tau_{\text{E}} = H\tau_{\text{E,scaling}}
+$$
+
+We then find the electron confinement time via the consistency relation from [`calculate_electron_species_consistent_energy_confinement_time()`](#global-confinement-time-consistency--calculate_electron_species_consistent_energy_confinement_time)
+
+$$
+\tau_{\text{e}} = \tau_{\text{E}}\left(\frac{W_{\text{e}} +\frac{W_{\text{i}}}{R}}{W_{\text{e}} +W_{\text{i}}}\right)
+$$
+
+The ion energy confinement time is then simply defined as:
+
+$$
+\tau_{\text{i}} = R\tau_{\text{e}}
+$$
+
+
+----------------
+
 ### Transport Powers
 
-After the confinement time scaling with $H$-factor correction has been calculated, the ion and electron transport power densities are found. `PROCESS` assumes the scaling confinement time to be equal to the ion and electron energy confinement time.
+After the confinement time scaling with $H$-factor correction has been calculated along with the specific ion and electron energy confinement times, the ion and electron transport power densities are found.
 
 This is simply the volume averaged thermal energy of the electron and ions divided by the $H$-factor corrected confinement time from the chosen scaling.
 
 $$
-\mathtt{pden\_ion\_transport\_loss\_mw} = \frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}}
+\overbrace{\left(\frac{P_{\text{i,loss}}}{V}\right)}^{\texttt{pden_ion_transport_loss_mw}} = \frac{3}{2}\frac{\frac{W_{\text{i}}}{V}}{\tau_{\text{E}}}
 $$
 
 $$
-\mathtt{pden\_electron\_transport\_loss\_mw} = \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}}
+\overbrace{\left(\frac{P_{\text{e,loss}}}{V}\right)}^{\texttt{pden_electron_transport_loss_mw}} = \frac{3}{2}\frac{\frac{W_{\text{e}}}{V}}{\tau_{\text{E}}}
 $$
 
-Here $\langle T_{\text{i}} \rangle$ and $\langle T_{\text{e}} \rangle$ are the ion and electron density weighted temperatures respectively.
+Here $W_{\text{e}}$ and $W_{\text{i}}$ are the electron and ion stored thermal energies respectively and $V$ is the plasma volume.
+
+
 
-Calculate the density and density weighted ratio:
+----------
+
+## Global confinement time consistency | `calculate_electron_species_consistent_energy_confinement_time()`
+
+It is possible to set a relative ratio between the electron and ion confinement times $\overbrace{R}^{\texttt{f_t_fuel_ion_electron_energy_confinement}} = \frac{\tau_\text{i}}{\tau_\text{e}}$ that still ensures the global energy confinment time is consistent.
+
+It is derived as:
+
+The global energy confinement time is defined as:
 
 $$
-\frac{n_{\text{i}}}{n_{\text{e}}}\frac{\langle T_{\text{i}} \rangle_{\text{n}}}{\langle T_{\text{e}} \rangle_{\text{n}}}
+\tau_{\text{E}} = \frac{W_{\text{e}} +W_{\text{i}}}{P_{\text{e}} +P_{\text{i}}}
 $$
 
-The density weighted global energy confinement time is then found in terms of this ratio:
+where $W_{\text{e}},W_{\text{i}}$ is the stored thermal energy of the electrons and ions respectively and $P_{\text{e}},P_{\text{i}}$ is the electron and ion energy transport powers.
+
+The electron and ion energy confinement times alone are:
 
 $$
-\tau_{\text{E}} = \frac{\frac{n_{\text{i}}}{n_{\text{e}}}\frac{\langle T_{\text{i}} \rangle_{\text{n}}}{\langle T_{\text{e}} \rangle_{\text{n}}} + 1}{\left(\frac{\frac{n_{\text{i}}}{n_{\text{e}}}\frac{\langle T_{\text{i}} \rangle_{\text{n}}}{\langle T_{\text{e}} \rangle_{\text{n}}}}{\tau_{\text{i}}}+\frac{1}{\tau_{\text{e}}}\right)}
+\tau_{\text{e}} = \frac{W_{\text{e}}}{P_{\text{e}}}
 $$
 
-----------
-
-## Key Constraints
+$$
+\tau_{\text{i}} = \frac{W_{\text{i}}}{P_{\text{i}}}
+$$
 
-### Global plasma power balance
+Substituting the individual ion and electron terms into the global term we get:
 
-This constraint can be activated by stating `icc = 2` in the input file.
+$$
+\tau_{\text{E}} = \frac{W_{\text{e}} +W_{\text{i}}}{\frac{W_{\text{e}}}{\tau_{\text{e}}} +\frac{W_{\text{i}}}{\tau_{\text{i}}}}
+$$
 
-This constraint ensures self consistency between the the transport loss power used for the confinement scalings and the calculated confinement time in relation to the plasmas total thermal energy:
+Substituting in our known ratio $\left(R = \frac{\tau_\text{i}}{\tau_\text{e}}\right)$:
 
 $$
-P_{\text{L}} = \frac{W}{\tau_{\text{E}}}
+\tau_{\text{E}} = \frac{W_{\text{e}} +W_{\text{i}}}{\frac{W_{\text{e}}}{\tau_{\text{e}}} +\frac{W_{\text{e}}}{R\tau_{\text{e}}}}
 $$
 
+Factor out the $\frac{1}{\tau_{\text{e}}}$ from the denominator:
+
 $$
-\underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}} + \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}}}_{\frac{W}{\tau_{\text{E}}}}  = \underbrace{\frac{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}{V_{\text{P}}} - \frac{P_{\text{rad}}}{V_{\text{p}}}}_{P_{\text{L}}}
+\tau_{\text{E}} = \tau_{\text{e}}\left(\frac{W_{\text{e}} +W_{\text{i}}}{W_{\text{e}} +\frac{W_{\text{i}}}{R}}\right)
 $$
 
-The $\frac{3}{2}n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}$ value is simply the volume averaged ion thermal energy density where $\langle T_{\text{i}} \rangle_{\text{n}}$ is the density weighted temperature. The same goes for the $\frac{3}{2}n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{e}}$ electron thermal energy density term. $\tau_{\text{E}}$ is the confinement time calculated from the chosen confinement scaling via `i_confinement_time`. 
+Finally re-arrange the equation to isolate $\tau_{\text{e}}$:
 
-The constraint uses the loss power and thermal densities hence the inclusion of the $V_{\text{p}}$ plasma volume term. The constraint is adapted depending on the condition of `i_rad_loss` which governs the radiation contribution to the loss power definition, see the [radiation and energy confinement section](#effect-of-radiation-on-energy-confinement) for more info. The injected heating and current drive contribution $P_{\text{HCD}}$ is also included or excluded depending if the plasma is deemed to be ignited with the `i_plasma_ignited` switch.
+$$
+\tau_{\text{e}} = \tau_{\text{E}}\left(\frac{W_{\text{e}} +\frac{W_{\text{i}}}{R}}{W_{\text{e}} +W_{\text{i}}}\right)
+$$
 
-**It is highly recommended to always have this constraint on as it is a global consistency checker**
+Using the above we can specify plasmas where $\tau_{\text{e}} \neq \tau_{\text{i}} \neq \tau_{\text{E}}$ and still ensure that our global plasma and species energy transport is self consistent.
 
 ----------
 
+## Key Constraints
+
 ### Lower limit on alpha particle confinement time ratio
 
 This constraint can be activated by stating `icc = 62` in the input file.

documentation/source/physics-models/plasma_power_balance.md

@@ -1,20 +1,84 @@
-# Plasma Core Power Balance
-
-The figure below shows the flow of power as calculated by the code.
-
-<figure markdown>
-![Power balance](../images/Overall-power-flow.png){ width="100%"}
-<figcaption>Figure 1: Machine build for D-shaped major components</figcaption>
-</figure>
-
-The primary sources of power are the fusion reactions themselves, ohmic power
-due to resistive heating within the plasma, and any auxiliary power provided
-for heating and current drive. The power carried by the fusion-generated 
-neutrons is lost from the plasma, but is deposited in the surrounding material. 
-A fraction `f_p_alpha_plasma_deposited` of the alpha particle power is assumed to stay within the
-plasma core to contribute to the plasma power balance. The sum of this core
-alpha power, any power carried by non-alpha charged particles, the ohmic power
-and any injected power, is converted into charged particle transport power
-($P_{\mbox{loss}}$) plus core radiation power, as shown in the Figure. The core 
-power balance calculation is turned on using constraint equation no. 2 (which 
-should therefore always be used).
+# Plasma Power Balance
+
+There are several equations that need to be satisfied in order to ensure that the plasma is in 
+power equilibrium. This includes the total energy and power leaving the plasma as a whole and then also the rate of energy transfer between the ion and electron species themselves.
+
+!!! warning "Global and species power balance"
+
+    It is not possible to have `icc = 2` for the global plasma power balance on at the same time as the ion `icc = 3` and electron `icc = 4` power balance. This is due to a solver error where the close to convergence the equality Jacobian loses rank. This is because the solver sees only two independant equality conditions but is asked to enforce three.
+
+    Using the ion and electron power balance constraints is preffered and still enforces the global power balance consistency. 
+
+## Global plasma power balance
+
+This constraint can be activated by stating `icc = 2` in the input file.
+
+This constraint ensures self consistency between the the transport loss power used for the confinement scalings and the calculated confinement time in relation to the plasmas total thermal energy:
+
+$$
+P_{\text{L}} = \frac{W}{\tau_{\text{E}}}
+$$
+
+$$
+\underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{i}}} + \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{e}}}}_{\frac{W}{\tau_{\text{E}}}}  = \underbrace{\frac{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}{V_{\text{P}}} - \frac{P_{\text{rad}}}{V_{\text{p}}}}_{P_{\text{L}}}
+$$
+
+The $\frac{3}{2}n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}$ value is simply the volume averaged ion thermal energy density where $\langle T_{\text{i}} \rangle_{\text{n}}$ is the density weighted temperature. The same goes for the $\frac{3}{2}n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{e}}$ electron thermal energy density term. $\tau_{\text{E}}$ is the confinement time calculated from the chosen confinement scaling via `i_confinement_time`. 
+
+The constraint uses the loss power and thermal densities hence the inclusion of the $V_{\text{p}}$ plasma volume term. The constraint is adapted depending on the condition of `i_rad_loss` which governs the radiation contribution to the loss power definition, see the [radiation and energy confinement section](#effect-of-radiation-on-energy-confinement) for more info. The injected heating and current drive contribution $P_{\text{HCD}}$ is also included or excluded depending if the plasma is deemed to be ignited with the `i_plasma_ignited` switch.
+
+**It is highly recommended to always have this constraint on as it is a global consistency checker**
+
+----------
+
+## Plasma Species Power Balance
+
+While the plasma may be in global power equilibrium, the electrons and ions are continually transferring energy to each other via Coulomb collisions and certain heating and current drive systems only heat particular species. Therefore in order to be in true equilibrium the net power transfer between all species in the plasma must also be equal to zero.
+
+!!! tip "Solving for particle power balance"
+
+    It is highly recommended to have `f_temp_plasma_ion_electron` which represent the ion to electron temperature ratio $\frac{T_i}{T_e}$, as an iteration variable. This allows the solver to have another degree of freedom upon which to change the fusion power and the equilibration power.
+
+-------------------
+
+## Electron Power Balance
+
+This constraint can be activated by stating `icc = 3` in the input file.
+
+For the electrons this internal power blance must be satisifed:
+
+$$
+\underbrace{P_{\Omega} + P_{\text{HCD,e}}+ (f_{e}\times(f_{\alpha}P_{\alpha}))}_{\text{Power gain}} = \underbrace{P_{\text{rad}} - P_{\text{ei}} - \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{e}}}}_{\text{Power loss}}
+$$
+
+where $P_{\Omega}$ is the plasma ohmic heating power, $P_{\text{HCD,e}}$ is the external heating and current drive power that only goes to electrons,$(f_{e}\times(f_{\alpha}P_{\alpha}))$ is the fraction of the coupled alpha particle power that goes to the electrons, $P_{\text{rad}}$ is the total radiation power given off by the plasma, $P_{\text{ei}}$ is the [electron-ion equilibration power](#ion-electron-equilibration-power-density--calculate_ion_electron_equilibration_power).
+
+**It is highly recommended to always have this constraint on as it is a plasma equilibrium checker**
+
+------------
+
+## Ion Power Balance
+
+This constraint can be activated by stating `icc = 4` in the input file.
+
+For the ions this internal power blance must be satisifed:
+
+$$
+\underbrace{P_{\text{HCD,i}} + (f_{i}\times(f_{\alpha}P_{\alpha})) + P_{\text{ei}}}_{\text{Power gain}} = \underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{i}}}}_{\text{Power Loss}}
+$$
+
+where $P_{\text{HCD,i}}$ is the external heating and current drive power that only goes to ions,$(f_{e}\times(f_{\alpha}P_{\alpha}))$ is the fraction of the coupled alpha particle power that goes to the ions, $P_{\text{ei}}$ is the [electron-ion equilibration power](#ion-electron-equilibration-power-density--calculate_ion_electron_equilibration_power).
+
+**It is highly recommended to always have this constraint on as it is a plasma equilibrium checker**
+
+-----------
+
+## Ion-electron equilibration power density | `calculate_ion_electron_equilibration_power()`
+
+The equilibration power $P_{\text{ei}}$ refers to the rate of energy transfer between plasma species—typically ions and electrons—due to elastic Coulomb collisions. It is driven by the temperature difference between the species. The equation below represents the rate of energy transfer from the electrons to the ions due to these elastic Coulomb collisions caused by the temperature difference between the electrons and the ions:
+
+$$
+P_{\text{ei}} = \frac{3}{2} \frac{n_e \text{e} (T_e - T_i)}{\tau_{\text{eq}}}
+$$
+
+where $n_e$ is the electron density, $T_i$ and $T_e$ is the ion and electron temperatures in $[\text{eV}]$ and $\tau_{\text{eq}}$ is the ion-electron equilibration time.

process/core/init.py

@@ -294,7 +294,7 @@ def check_process(inputs, data):  # noqa: ARG001
         )
 
     # Deprecate constraints
-    for depcrecated_constraint in [3, 4, 10, 74, 42]:
+    for depcrecated_constraint in [10, 74, 42]:
         if (
             data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns]
             == depcrecated_constraint
@@ -1160,6 +1160,17 @@ def check_process(inputs, data):  # noqa: ARG001
             "turn off CS temperature margin constraint icc = 60 when using REBCO"
         )
 
+    # Cannot use temperature margin constraint with REBCO CS coils
+    if (
+        (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 2).any()
+        and (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 3).any()
+        and (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 4).any()
+    ):
+        raise ProcessValidationError(
+            "Cannot have global power balance on as well as individual particle power "
+            "balance, this will cause a solver error. Either use icc = 2 alone or icc = 3 and 4 alone"
+        )
+
     # Cold end of the cryocooler should be colder than the TF
     if data.tfcoil.temp_tf_cryo > data.tfcoil.tftmp:
         raise ProcessValidationError("temp_tf_cryo should be lower than tftmp")

process/core/input.py

@@ -1184,6 +1184,9 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "force_vmcon_inequality_tolerance": InputVariable(
         "numerics", float, range=(0.0, 1e10)
     ),
+    "f_t_fuel_ion_electron_energy_confinement": InputVariable(
+        "physics", float, range=(0.1, 5.0)
+    ),
 }