diff --git a/documentation/source/physics-models/detailed_physics.md b/documentation/source/physics-models/detailed_physics.md index bc1844a899..b928feb8f5 100644 --- a/documentation/source/physics-models/detailed_physics.md +++ b/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} +$$ diff --git a/documentation/source/physics-models/plasma_confinement.md b/documentation/source/physics-models/plasma_confinement.md index 80124f1f89..24277c770d 100644 --- a/documentation/source/physics-models/plasma_confinement.md +++ b/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. diff --git a/documentation/source/physics-models/plasma_power_balance.md b/documentation/source/physics-models/plasma_power_balance.md index cdab37b71c..6cd5e9f588 100644 --- a/documentation/source/physics-models/plasma_power_balance.md +++ b/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. - -
-![Power balance](../images/Overall-power-flow.png){ width="100%"} -
Figure 1: Machine build for D-shaped major components
-
- -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). \ No newline at end of file +# 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. \ No newline at end of file diff --git a/process/core/init.py b/process/core/init.py index 5aba77d345..4a5d108d6c 100644 --- a/process/core/init.py +++ b/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") diff --git a/process/core/input.py b/process/core/input.py index 1d7c132b6a..270dd3ba3a 100644 --- a/process/core/input.py +++ b/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) + ), } diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index e8b2fea071..73bb5fe489 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -2818,17 +2818,17 @@ def plot_main_plasma_information( f"$\\mathbf{{Confinement:}}$\n \n" f"Confinement scaling law: {mfile.get('tauelaw', scan=scan)}\n" f"Confinement $H$ factor: {mfile.get('hfact', scan=scan):.4f}\n" - f"Energy confinement time from scaling: {mfile.get('t_energy_confinement', scan=scan):.4f} s\n" + f"$\\tau_{{\\text{{E}}}}$: {mfile.get('t_energy_confinement', scan=scan):.4f} s | $\\tau_{{\\text{{e}}}}$: {mfile.get('t_electron_energy_confinement', scan=scan):.4f} s | $\\tau_{{\\text{{i}}}}$: {mfile.get('t_ion_energy_confinement', scan=scan):.4f} s | $\\tau_{{\\text{{i}}}}/\\tau_{{\\text{{e}}}}$: {mfile.get('f_t_fuel_ion_electron_energy_confinement', scan=scan):.4f} \n" f"Fusion double product: {mfile.get('ntau', scan=scan):.4e} s/m³\n" f"Lawson Triple product: {mfile.get('nttau', scan=scan):.4e} keV·s/m³\n" f"Transport loss power assumed in scaling law: {mfile.get('p_plasma_loss_mw', scan=scan):.4f} MW\n" - f"Plasma thermal energy (inc. $\\alpha$), $W$: {mfile.get('e_plasma_beta', scan=scan) / 1e9:.4f} GJ\n" - f"Alpha particle confinement time: {mfile.get('t_alpha_confinement', scan=scan):.4f} s | $\\tau_{{\\alpha}}/\\tau_{{e}}$: {mfile.get('f_t_alpha_energy_confinement', scan=scan):.4f}" + f"$W_{{\\beta}}$: {mfile.get('e_plasma_beta', scan=scan) / 1e9:.4f} GJ | $W_{{th}}$: {mfile.get('e_plasma_thermal_total', scan=scan) / 1e9:.4f} GJ | $W_{{e}}$: {mfile.get('e_plasma_electrons_thermal', scan=scan) / 1e9:.4f} GJ | $W_{{i}}$: {mfile.get('e_plasma_ions_thermal', scan=scan) / 1e9:.4f} GJ\n" + f"Particle confinement: $\\tau_{{\\alpha}}$: {mfile.get('t_alpha_confinement', scan=scan):.4f} s | $\\tau_{{\\alpha}}/\\tau_{{e}}$: {mfile.get('f_t_alpha_energy_confinement', scan=scan):.4f}" ) axis.text( 0.025, - 0.57, + 0.575, textstr_confinement, fontsize=9, verticalalignment="top", @@ -11811,8 +11811,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m3/sec\n\n" f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m3\n" f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m3\n\n" - f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m3\n" - f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n" + f"Power per unit volume transferred to electrons: {mfile.get('pden_alpha_heating_electrons_mw', scan=scan):.4e} MW/m3\n" + f"Power per unit volume transferred to ions: {mfile.get('pden_alpha_heating_ions_mw', scan=scan):.4e} MW/m3\n\n" ) axis.text( @@ -14090,6 +14090,63 @@ def plot_resistivity_profile(axis: plt.Axes, mfile_data: MFile, scan: int) -> No axis.legend() +def plot_plasma_equilibration_time_profile( + axis: plt.Axes, mfile_data: MFile, scan: int +) -> None: + """Plot the plasma equilibration time on the given axis.""" + t_plasma_electron_deuteron_equilibration_profile = [ + mfile_data.data[f"t_plasma_electron_deuteron_equilibration_profile{i}"].get_scan( + scan + ) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + t_plasma_electron_triton_equilibration_profile = [ + mfile_data.data[f"t_plasma_electron_triton_equilibration_profile{i}"].get_scan( + scan + ) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + t_plasma_electron_alpha_thermal_equilibration_profile = [ + mfile_data.data[ + f"t_plasma_electron_alpha_thermal_equilibration_profile{i}" + ].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + + axis.plot( + np.linspace(0, 1, len(t_plasma_electron_deuteron_equilibration_profile)), + t_plasma_electron_deuteron_equilibration_profile, + color="pink", + linestyle="-", + label=r"$\tau_{e-D,eq}$", + ) + + axis.plot( + np.linspace(0, 1, len(t_plasma_electron_triton_equilibration_profile)), + t_plasma_electron_triton_equilibration_profile, + color="green", + linestyle="-", + label=r"$\tau_{e-T,eq}$", + ) + + axis.plot( + np.linspace(0, 1, len(t_plasma_electron_alpha_thermal_equilibration_profile)), + t_plasma_electron_alpha_thermal_equilibration_profile, + color="red", + linestyle="-", + label=r"$\tau_{e-\alpha,thermal,eq}$", + ) + + axis.set_yscale("log") + axis.set_ylabel("Equilibration Time [s]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + def plot_equality_constraint_equations(axis: plt.Axes, m_file_data: MFile, scan: int): """Plot the equality constraints for a solution and their normalised residuals @@ -14785,6 +14842,13 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int f"$\\langle\\eta_{{Spitzer}}\\rangle$: {mfile.get('res_plasma_fuel_spitzer_vol_avg', scan=scan):.4e} $\\Omega\\mathrm{{m}}$" ) + textstr_equilibriation = ( + f"$\\mathbf{{Equilibration \\ Times:}}$\n\n" + f"$\\langle\\tau_{{e-D,Equil}}\\rangle$: {mfile.get('t_plasma_electron_deuteron_equilibration_vol_avg', scan=scan):.4e} s\n" + f"$\\langle\\tau_{{e-T,Equil}}\\rangle$: {mfile.get('t_plasma_electron_triton_equilibration_vol_avg', scan=scan):.4e} s\n" + f"$\\langle\\tau_{{e-\\alpha,Equil}}\\rangle$: {mfile.get('t_plasma_electron_alpha_thermal_equilibration_vol_avg', scan=scan):.4e} s" + ) + light_yellow_box = { "boxstyle": "round", "facecolor": "lightyellow", @@ -14825,6 +14889,17 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int bbox=light_yellow_box, ) + axis.text( + 0.65, + 0.45, + textstr_equilibriation, + fontsize=9, + verticalalignment="top", + horizontalalignment="left", + transform=fig.transFigure, + bbox=light_yellow_box, + ) + light_cyan_box = { "boxstyle": "round", "facecolor": "lightcyan", @@ -16114,6 +16189,10 @@ def _add_page(name: str | None = None): plot_resistivity_profile(pages["detailed_params"].add_subplot(232), m_file, scan) + plot_plasma_equilibration_time_profile( + pages["detailed_params"].add_subplot(233), m_file, scan + ) + plot_detailed_plasma_parameters( pages["detailed_params"].add_subplot(233), fig=pages["detailed_params"], diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 3c2406507e..b40dc8c8b5 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -279,24 +279,25 @@ def constraint_equation_2(constraint_registration, data): """ # pscaling: total transport power per volume (MW/m3) - pscaling = ( - data.physics.pden_electron_transport_loss_mw - + data.physics.pden_ion_transport_loss_mw - ) + p_plasma_loss = ( + ( + data.physics.eden_plasma_electrons_thermal_vol_avg + + data.physics.eden_plasma_ions_thermal_vol_avg + ) + / 1e6 + ) / data.physics.t_energy_confinement # Total power lost is scaling power plus radiation: if data.physics.i_rad_loss == 0: - pnumerator = pscaling + data.physics.pden_plasma_rad_mw + p_plasma_loss += data.physics.pden_plasma_rad_mw elif data.physics.i_rad_loss == 1: - pnumerator = pscaling + data.physics.pden_plasma_core_rad_mw - else: - pnumerator = pscaling + p_plasma_loss += data.physics.pden_plasma_core_rad_mw # if plasma not ignited include injected power if ( PlasmaIgnitionModel(data.physics.i_plasma_ignited) == PlasmaIgnitionModel.NON_IGNITED ): - pdenom = ( + p_plasma_heating = ( data.physics.f_p_alpha_plasma_deposited * data.physics.pden_alpha_total_mw + data.physics.pden_non_alpha_charged_mw + data.physics.pden_plasma_ohmic_mw @@ -304,13 +305,13 @@ def constraint_equation_2(constraint_registration, data): ) else: # if plasma ignited - pdenom = ( + p_plasma_heating = ( data.physics.f_p_alpha_plasma_deposited * data.physics.pden_alpha_total_mw + data.physics.pden_non_alpha_charged_mw + data.physics.pden_plasma_ohmic_mw ) - return eq(pnumerator, pdenom, constraint_registration) + return eq(p_plasma_loss, p_plasma_heating, constraint_registration) @ConstraintManager.register_constraint(3, "MW/m³", "=") @@ -322,11 +323,11 @@ def constraint_equation_3(constraint_registration, data): pden_ion_transport_loss_mw: ion transport power per volume (MW/m3) - pden_ion_electron_equilibration_mw: + pden_ion_electron_equilibration_vol_avg_mw: ion/electron equilibration power per volume (MW/m3) f_p_alpha_plasma_deposited: fraction of alpha power deposited in plasma - f_pden_alpha_ions_mw: + pden_alpha_heating_ions_mw: alpha power per volume to ions (MW/m3) p_hcd_injected_ions_mw: auxiliary injected power to ions (MW) @@ -338,26 +339,31 @@ def constraint_equation_3(constraint_registration, data): PlasmaIgnitionModel(data.physics.i_plasma_ignited) == PlasmaIgnitionModel.NON_IGNITED ): + p_ion_heating = ( + data.physics.pden_alpha_heating_ions_mw + + (data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma) + + data.physics.pden_ion_electron_equilibration_vol_avg_mw + + data.physics.pden_non_alpha_charged_mw + ) + + p_ion_loss = data.physics.pden_ion_transport_loss_mw return eq( - ( - data.physics.pden_ion_transport_loss_mw - + data.physics.pden_ion_electron_equilibration_mw - ), - ( - data.physics.f_p_alpha_plasma_deposited - * data.physics.f_pden_alpha_ions_mw - + data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma - ), + (p_ion_loss), + (p_ion_heating), constraint_registration, ) # Plasma ignited + p_ion_heating = ( + data.physics.pden_alpha_heating_ions_mw + + data.physics.pden_ion_electron_equilibration_vol_avg_mw + + data.physics.pden_non_alpha_charged_mw + ) + + p_ion_loss = data.physics.pden_ion_transport_loss_mw return eq( - ( - data.physics.pden_ion_transport_loss_mw - + data.physics.pden_ion_electron_equilibration_mw - ), - (data.physics.f_p_alpha_plasma_deposited * data.physics.f_pden_alpha_ions_mw), + (p_ion_loss), + (p_ion_heating), constraint_registration, ) @@ -385,9 +391,9 @@ def constraint_equation_4(constraint_registration, data): total core radiation power per volume (MW/m3) f_p_alpha_plasma_deposited: fraction of alpha power deposited in plasma - f_pden_alpha_electron_mw: + pden_alpha_heating_electrons_mw: alpha power per volume to electrons (MW/m3) - pden_ion_electron_equilibration_mw: + pden_ion_electron_equilibration_vol_avg_mw: ion/electron equilibration power per volume (MW/m3) p_hcd_injected_electrons_mw: auxiliary injected power to electrons (MW) @@ -395,35 +401,33 @@ def constraint_equation_4(constraint_registration, data): """ # pscaling: total transport power per volume (MW/m3) - pscaling = data.physics.pden_electron_transport_loss_mw + p_electron_loss = ( + data.physics.pden_electron_transport_loss_mw + + data.physics.pden_ion_electron_equilibration_vol_avg_mw + ) # Total power lost is scaling power plus radiation: if data.physics.i_rad_loss == 0: - pnumerator = pscaling + data.physics.pden_plasma_rad_mw + p_electron_loss += data.physics.pden_plasma_rad_mw elif data.physics.i_rad_loss == 1: - pnumerator = pscaling + data.physics.pden_plasma_core_rad_mw - else: - pnumerator = pscaling + p_electron_loss += data.physics.pden_plasma_core_rad_mw # if plasma not ignited include injected power if ( PlasmaIgnitionModel(data.physics.i_plasma_ignited) == PlasmaIgnitionModel.NON_IGNITED ): - pdenom = ( - data.physics.f_p_alpha_plasma_deposited - * data.physics.f_pden_alpha_electron_mw - + data.physics.pden_ion_electron_equilibration_mw - + data.current_drive.p_hcd_injected_electrons_mw / data.physics.vol_plasma + p_electron_heating = ( + data.physics.pden_alpha_heating_electrons_mw + + (data.current_drive.p_hcd_injected_electrons_mw / data.physics.vol_plasma) + + data.physics.pden_plasma_ohmic_mw ) else: # if plasma ignited - pdenom = ( - data.physics.f_p_alpha_plasma_deposited - * data.physics.f_pden_alpha_electron_mw - + data.physics.pden_ion_electron_equilibration_mw + p_electron_heating = ( + data.physics.pden_alpha_heating_electrons_mw + + data.physics.pden_plasma_ohmic_mw ) - - return eq(pnumerator, pdenom, constraint_registration) + return eq(p_electron_loss, p_electron_heating, constraint_registration) @ConstraintManager.register_constraint(5, "/m³", "<=") diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index ac79dc2789..5ab1998a7f 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -235,6 +235,10 @@ class IterationVariable: 175: IterationVariable("kappa", "physics", 0.00, 10.00), 176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30), 177: IterationVariable("f_a_tf_turn_cable_space_extra_void", "tfcoil", 0.01, 1.0), + 178: IterationVariable("f_temp_plasma_ion_electron", "physics", 0.0, 2.0), + 179: IterationVariable( + "f_t_fuel_ion_electron_energy_confinement", "physics", 0.1, 5.0 + ), } diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 28e37dee79..596b821dfc 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -684,14 +684,14 @@ class PhysicsData: f_c_plasma_inductive: float = 0.0 """fraction of plasma current produced inductively""" - f_alpha_electron: float = 0.0 - """fraction of alpha energy to electrons""" + f_p_alpha_total_electron: float = 0.0 + """Fraction of total alpha power used to heat the electrons""" f_p_alpha_plasma_deposited: float = 0.95 """Fraction of alpha power deposited in plasma. Default of 0.95 taken from https://doi.org/10.1088/0029-5515/39/12/305.""" - f_alpha_ion: float = 0.0 - """fraction of alpha power to ions""" + f_p_alpha_total_ions: float = 0.0 + """Fraction of total alpha power used to heat the ions""" f_plasma_fuel_deuterium: float = 0.5 """Plasma deuterium fuel fraction""" @@ -1020,19 +1020,25 @@ class PhysicsData: """branching ratio for DD -> T""" pden_plasma_alpha_mw: float = 0.0 - """Alpha power per volume just from plasma [MW/m3]""" + """Alpha power per volume just from plasma [MW/m³]""" pden_alpha_total_mw: float = 0.0 - """Alpha power per volume from plasma and beams [MW/m3]""" + """Alpha power per volume from plasma and beams [MW/m³]""" - f_pden_alpha_electron_mw: float = 0.0 - """Alpha power per volume to electrons [MW/m3]""" + pden_alpha_heating_electrons_mw: float = 0.0 + """Alpha power heating per unit volume to electrons [MW/m³]""" + + p_alpha_heating_electrons_mw: float = 0.0 + """Alpha power heating to electrons [MW]""" p_fw_alpha_mw: float = 0.0 """alpha power escaping plasma and reaching first wall (MW)""" - f_pden_alpha_ions_mw: float = 0.0 - """alpha power per volume to ions (MW/m3)""" + pden_alpha_heating_ions_mw: float = 0.0 + """Alpha power heating per unit volume to ions [MW/m³]""" + + p_alpha_heating_ions_mw: float = 0.0 + """Alpha power heating to ions [MW]""" p_plasma_alpha_mw: float = 0.0 """Alpha power from only the plasma (MW)""" @@ -1109,8 +1115,11 @@ class PhysicsData: pflux_fw_rad_mw: float = 0.0 """Nominal mean radiation load on inside surface of reactor (MW/m2)""" - pden_ion_electron_equilibration_mw: float = 0.0 - """ion/electron equilibration power per volume (MW/m3)""" + pden_ion_electron_equilibration_vol_avg_mw: float = 0.5 + """Ion/electron volume averaged equilibration power per unit volume [MW/m³]""" + + p_ion_electron_equilibration_vol_avg_mw: float = 80.0 + """Ion/electron volume averaged equilibration power [MW]""" plasma_current: float = 0.0 """Plasma current (Iₚ) [A]""" @@ -1349,6 +1358,9 @@ class PhysicsData: f_t_alpha_energy_confinement: float = 0.0 """Alpha particle to energy confinement time ratio""" + f_t_fuel_ion_electron_energy_confinement: float = 1.0 + """Fuel ion to electron energy confinement time ratio""" + temp_plasma_electron_vol_avg_kev: float = 12.9 """Plasma volume averaged electron temperature (⟨Tₑ⟩) [keV] (`iteration variable 4`)""" @@ -1581,6 +1593,30 @@ class PhysicsData: t_plasma_electron_alpha_thermal_collision_vol_avg: float = 0.0 """Volume averaged electron-alpha collision time in plasma (s)""" + t_plasma_electron_deuteron_equilibration_vol_avg: float = 3.0 + """Volume averaged electron-deuteron equilibration time in plasma (s)""" + + t_plasma_electron_deuteron_equilibration_profile: list[float] = field( + default_factory=list + ) + """Profile of electron-deuteron equilibration time in plasma (s)""" + + t_plasma_electron_triton_equilibration_vol_avg: float = 4.5 + """Volume averaged electron-triton equilibration time in plasma (s)""" + + t_plasma_electron_triton_equilibration_profile: list[float] = field( + default_factory=list + ) + """Profile of electron-triton equilibration time in plasma (s)""" + + t_plasma_electron_alpha_thermal_equilibration_vol_avg: float = 4.5 + """Volume averaged electron-thermal alphas equilibration time in plasma (s)""" + + t_plasma_electron_alpha_thermal_equilibration_profile: list[float] = field( + default_factory=list + ) + """Profile of electron-thermal alphas equilibration time in plasma (s)""" + freq_plasma_electron_electron_collision_profile: list[float] = field( default_factory=list ) diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py index 2e9dd2ae8a..d7c8da094d 100644 --- a/process/models/physics/confinement_time.py +++ b/process/models/physics/confinement_time.py @@ -220,7 +220,7 @@ def calculate_confinement_time( if ( model == ConfinementTimeModel.USER_INPUT ): # t_electron_energy_confinement is an input - t_electron_confinement = self.data.physics.tauee_in + t_global_confinement_scaling = self.data.physics.tauee_in # ======================================================================== @@ -228,7 +228,7 @@ def calculate_confinement_time( if ( model == ConfinementTimeModel.NEO_ALCATOR ): # t_electron_energy_confinement is an input - t_electron_confinement = self.neo_alcator_confinement_time( + t_global_confinement_scaling = self.neo_alcator_confinement_time( n20, rminor, rmajor, qstar ) @@ -236,13 +236,15 @@ def calculate_confinement_time( # "Mirnov"-like scaling (H-mode) elif model == ConfinementTimeModel.MIRNOV: # Mirnov scaling (H-mode) - t_electron_confinement = self.mirnov_confinement_time(rminor, kappa95, pcur) + t_global_confinement_scaling = self.mirnov_confinement_time( + rminor, kappa95, pcur + ) # ======================================================================== # Merezhkin-Mukhovatov (MM) OH/L-mode scaling elif model == ConfinementTimeModel.MEREZHKIN_MUHKOVATOV: - t_electron_confinement = self.merezhkin_muhkovatov_confinement_time( + t_global_confinement_scaling = self.merezhkin_muhkovatov_confinement_time( rmajor, rminor, kappa95, @@ -256,7 +258,7 @@ def calculate_confinement_time( # Shimomura (S) optimized H-mode scaling elif model == ConfinementTimeModel.SHIMOMURA: - t_electron_confinement = self.shimomura_confinement_time( + t_global_confinement_scaling = self.shimomura_confinement_time( rmajor, rminor, b_plasma_toroidal_on_axis, kappa95, m_fuel_amu ) @@ -264,7 +266,7 @@ def calculate_confinement_time( # Kaye-Goldston scaling (L-mode) elif model == ConfinementTimeModel.KAYE_GOLDSTON: - t_electron_confinement = self.kaye_goldston_confinement_time( + t_global_confinement_scaling = self.kaye_goldston_confinement_time( pcur, rmajor, rminor, @@ -279,7 +281,7 @@ def calculate_confinement_time( # ITER Power scaling - ITER 89-P (L-mode) elif model == ConfinementTimeModel.ITER_89P: - t_electron_confinement = self.iter_89p_confinement_time( + t_global_confinement_scaling = self.iter_89p_confinement_time( pcur, rmajor, rminor, @@ -294,7 +296,7 @@ def calculate_confinement_time( # ITER Offset linear scaling - ITER 89-O (L-mode) elif model == ConfinementTimeModel.ITER_89_0: - t_electron_confinement = self.iter_89_0_confinement_time( + t_global_confinement_scaling = self.iter_89_0_confinement_time( pcur, rmajor, rminor, @@ -308,7 +310,7 @@ def calculate_confinement_time( # Rebut-Lallia offset linear scaling (L-mode) elif model == ConfinementTimeModel.REBUT_LALLIA: - t_electron_confinement = self.rebut_lallia_confinement_time( + t_global_confinement_scaling = self.rebut_lallia_confinement_time( rminor, rmajor, kappa, @@ -324,7 +326,7 @@ def calculate_confinement_time( # Goldston scaling (L-mode) elif model == ConfinementTimeModel.GOLDSTON: # Goldston scaling (L-mode) - t_electron_confinement = self.goldston_confinement_time( + t_global_confinement_scaling = self.goldston_confinement_time( pcur, rmajor, rminor, kappa95, m_fuel_amu, p_plasma_loss_mw ) @@ -332,7 +334,7 @@ def calculate_confinement_time( # T-10 scaling (L-mode) elif model == ConfinementTimeModel.T_10: - t_electron_confinement = self.t10_confinement_time( + t_global_confinement_scaling = self.t10_confinement_time( dnla20, rmajor, qstar, @@ -348,7 +350,7 @@ def calculate_confinement_time( # JAERI / Odajima-Shimomura L-mode scaling elif model == ConfinementTimeModel.JAERI: # JAERI scaling - t_electron_confinement = self.jaeri_confinement_time( + t_global_confinement_scaling = self.jaeri_confinement_time( kappa95, rminor, m_fuel_amu, @@ -365,7 +367,7 @@ def calculate_confinement_time( # Kaye "big" L-mode scaling (based only on big tokamak data) elif model == ConfinementTimeModel.KAYE_BIG: - t_electron_confinement = self.kaye_big_confinement_time( + t_global_confinement_scaling = self.kaye_big_confinement_time( rmajor, rminor, b_plasma_toroidal_on_axis, @@ -380,7 +382,7 @@ def calculate_confinement_time( # ITER H90-P H-mode scaling elif model == ConfinementTimeModel.ITER_H90_P: - t_electron_confinement = self.iter_h90_p_confinement_time( + t_global_confinement_scaling = self.iter_h90_p_confinement_time( pcur, rmajor, rminor, @@ -395,7 +397,7 @@ def calculate_confinement_time( # Minimum of ITER 89-P and ITER 89-O elif model == ConfinementTimeModel.MINIMUM_OF_ITER_89P_AND_ITER_89_0: - t_electron_confinement = min( + t_global_confinement_scaling = min( self.iter_89p_confinement_time( pcur, rmajor, @@ -422,7 +424,7 @@ def calculate_confinement_time( # Riedel scaling (L-mode) elif model == ConfinementTimeModel.RIEDEL_L: - t_electron_confinement = self.riedel_l_confinement_time( + t_global_confinement_scaling = self.riedel_l_confinement_time( pcur, rmajor, rminor, @@ -436,7 +438,7 @@ def calculate_confinement_time( # Christiansen et al scaling (L-mode) elif model == ConfinementTimeModel.CHRISTIANSEN: - t_electron_confinement = self.christiansen_confinement_time( + t_global_confinement_scaling = self.christiansen_confinement_time( pcur, rmajor, rminor, @@ -451,7 +453,7 @@ def calculate_confinement_time( # Lackner-Gottardi scaling (L-mode) elif model == ConfinementTimeModel.LACKNER_GOTTARDI: - t_electron_confinement = self.lackner_gottardi_confinement_time( + t_global_confinement_scaling = self.lackner_gottardi_confinement_time( pcur, rmajor, rminor, @@ -465,7 +467,7 @@ def calculate_confinement_time( # Neo-Kaye scaling (L-mode) elif model == ConfinementTimeModel.NEO_KAYE: - t_electron_confinement = self.neo_kaye_confinement_time( + t_global_confinement_scaling = self.neo_kaye_confinement_time( pcur, rmajor, rminor, @@ -479,7 +481,7 @@ def calculate_confinement_time( # Riedel scaling (H-mode) elif model == ConfinementTimeModel.RIEDEL_H: - t_electron_confinement = self.riedel_h_confinement_time( + t_global_confinement_scaling = self.riedel_h_confinement_time( pcur, rmajor, rminor, @@ -494,7 +496,7 @@ def calculate_confinement_time( # Amended version of ITER H90-P law elif model == ConfinementTimeModel.ITER_H90_P_AMENDED: - t_electron_confinement = self.iter_h90_p_amended_confinement_time( + t_global_confinement_scaling = self.iter_h90_p_amended_confinement_time( pcur, b_plasma_toroidal_on_axis, m_fuel_amu, @@ -507,7 +509,7 @@ def calculate_confinement_time( # Sudo et al. scaling (stellarators/heliotron) elif model == ConfinementTimeModel.SUDO_ET_AL: - t_electron_confinement = self.sudo_et_al_confinement_time( + t_global_confinement_scaling = self.sudo_et_al_confinement_time( rmajor, rminor, dnla20, @@ -519,7 +521,7 @@ def calculate_confinement_time( # Gyro-reduced Bohm scaling elif model == ConfinementTimeModel.GYRO_REDUCED_BOHM: - t_electron_confinement = self.gyro_reduced_bohm_confinement_time( + t_global_confinement_scaling = self.gyro_reduced_bohm_confinement_time( b_plasma_toroidal_on_axis, dnla20, p_plasma_loss_mw, @@ -531,20 +533,22 @@ def calculate_confinement_time( # Lackner-Gottardi stellarator scaling elif model == ConfinementTimeModel.LACKNER_GOTTARDI_STELLARATOR: - t_electron_confinement = self.lackner_gottardi_stellarator_confinement_time( - rmajor, - rminor, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - q95, + t_global_confinement_scaling = ( + self.lackner_gottardi_stellarator_confinement_time( + rmajor, + rminor, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + q95, + ) ) # ========================================================================== # ITER_93 ELM-free H-mode scaling elif model == ConfinementTimeModel.ITER_93H: - t_electron_confinement = self.iter_93h_confinement_time( + t_global_confinement_scaling = self.iter_93h_confinement_time( pcur, b_plasma_toroidal_on_axis, p_plasma_loss_mw, @@ -563,7 +567,7 @@ def calculate_confinement_time( # ELM-free: ITERH-97P elif model == ConfinementTimeModel.ITER_H97P: - t_electron_confinement = self.iter_h97p_confinement_time( + t_global_confinement_scaling = self.iter_h97p_confinement_time( pcur, b_plasma_toroidal_on_axis, p_plasma_loss_mw, @@ -578,7 +582,7 @@ def calculate_confinement_time( # ELMy: ITERH-97P(y) elif model == ConfinementTimeModel.ITER_H97P_ELMY: - t_electron_confinement = self.iter_h97p_elmy_confinement_time( + t_global_confinement_scaling = self.iter_h97p_elmy_confinement_time( pcur, b_plasma_toroidal_on_axis, p_plasma_loss_mw, @@ -593,7 +597,7 @@ def calculate_confinement_time( # ITER-96P (= ITER-97L) L-mode scaling elif model == ConfinementTimeModel.ITER_96P: - t_electron_confinement = self.iter_96p_confinement_time( + t_global_confinement_scaling = self.iter_96p_confinement_time( pcur, b_plasma_toroidal_on_axis, kappa95, @@ -609,7 +613,7 @@ def calculate_confinement_time( # Valovic modified ELMy-H mode scaling # WARNING: No reference found for this scaling. This may not be its real name elif model == ConfinementTimeModel.VALOVIC_ELMY: - t_electron_confinement = self.valovic_elmy_confinement_time( + t_global_confinement_scaling = self.valovic_elmy_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -625,7 +629,7 @@ def calculate_confinement_time( # Kaye PPPL Workshop April 1998 L-mode scaling # WARNING: No reference found for this scaling. This may not be its real name elif model == ConfinementTimeModel.KAYE: - t_electron_confinement = self.kaye_confinement_time( + t_global_confinement_scaling = self.kaye_confinement_time( pcur, b_plasma_toroidal_on_axis, kappa, @@ -641,7 +645,7 @@ def calculate_confinement_time( # ITERH-PB98P(y), ELMy H-mode scaling # WARNING: No reference found for this scaling. This may not be its real name elif model == ConfinementTimeModel.ITER_PB98P_Y: - t_electron_confinement = self.iter_pb98py_confinement_time( + t_global_confinement_scaling = self.iter_pb98py_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -656,7 +660,7 @@ def calculate_confinement_time( # IPB98(y), ELMy H-mode scaling elif model == ConfinementTimeModel.IPB98_Y: - t_electron_confinement = self.iter_ipb98y_confinement_time( + t_global_confinement_scaling = self.iter_ipb98y_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -671,7 +675,7 @@ def calculate_confinement_time( # IPB98(y,1), ELMy H-mode scaling elif model == ConfinementTimeModel.ITER_IPB98Y1: - t_electron_confinement = self.iter_ipb98y1_confinement_time( + t_global_confinement_scaling = self.iter_ipb98y1_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -686,7 +690,7 @@ def calculate_confinement_time( # IPB98(y,2), ELMy H-mode scaling elif model == ConfinementTimeModel.ITER_IPB98Y2: - t_electron_confinement = self.iter_ipb98y2_confinement_time( + t_global_confinement_scaling = self.iter_ipb98y2_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -701,7 +705,7 @@ def calculate_confinement_time( # IPB98(y,3), ELMy H-mode scaling elif model == ConfinementTimeModel.ITER_IPB98Y3: - t_electron_confinement = self.iter_ipb98y3_confinement_time( + t_global_confinement_scaling = self.iter_ipb98y3_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -716,7 +720,7 @@ def calculate_confinement_time( # IPB98(y,4), ELMy H-mode scaling elif model == ConfinementTimeModel.ITER_IPB98Y4: - t_electron_confinement = self.iter_ipb98y4_confinement_time( + t_global_confinement_scaling = self.iter_ipb98y4_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -733,7 +737,7 @@ def calculate_confinement_time( elif model == ConfinementTimeModel.ISS95_STELLARATOR: # dummy argument q95 is actual argument iotabar for stellarators iotabar = q95 - t_electron_confinement = self.iss95_stellarator_confinement_time( + t_global_confinement_scaling = self.iss95_stellarator_confinement_time( rminor, rmajor, dnla19, @@ -748,7 +752,7 @@ def calculate_confinement_time( elif model == ConfinementTimeModel.ISS04_STELLARATOR: # dummy argument q95 is actual argument iotabar for stellarators iotabar = q95 - t_electron_confinement = self.iss04_stellarator_confinement_time( + t_global_confinement_scaling = self.iss04_stellarator_confinement_time( rminor, rmajor, dnla19, @@ -761,7 +765,7 @@ def calculate_confinement_time( # DS03 beta-independent H-mode scaling elif model == ConfinementTimeModel.DS03: - t_electron_confinement = self.ds03_confinement_time( + t_global_confinement_scaling = self.ds03_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -776,7 +780,7 @@ def calculate_confinement_time( # Murari "Non-power law" scaling elif model == ConfinementTimeModel.MURARI: - t_electron_confinement = self.murari_confinement_time( + t_global_confinement_scaling = self.murari_confinement_time( pcur, rmajor, self.data.physics.kappa_ipb, @@ -789,7 +793,7 @@ def calculate_confinement_time( # Petty08, beta independent dimensionless scaling elif model == ConfinementTimeModel.PETTY08: - t_electron_confinement = self.petty08_confinement_time( + t_global_confinement_scaling = self.petty08_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -803,7 +807,7 @@ def calculate_confinement_time( # Lang high density relevant confinement scaling elif model == ConfinementTimeModel.LANG_HIGH_DENSITY: - t_electron_confinement = self.lang_high_density_confinement_time( + t_global_confinement_scaling = self.lang_high_density_confinement_time( plasma_current, b_plasma_toroidal_on_axis, nd_plasma_electron_line, @@ -821,7 +825,7 @@ def calculate_confinement_time( # Hubbard 2017 I-mode confinement time scaling - nominal elif model == ConfinementTimeModel.HUBBARD_NOMINAL: - t_electron_confinement = self.hubbard_nominal_confinement_time( + t_global_confinement_scaling = self.hubbard_nominal_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla20, @@ -832,7 +836,7 @@ def calculate_confinement_time( # Hubbard 2017 I-mode confinement time scaling - lower elif model == ConfinementTimeModel.HUBBARD_LOWER: - t_electron_confinement = self.hubbard_lower_confinement_time( + t_global_confinement_scaling = self.hubbard_lower_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla20, @@ -843,7 +847,7 @@ def calculate_confinement_time( # Hubbard 2017 I-mode confinement time scaling - upper elif model == ConfinementTimeModel.HUBBARD_UPPER: - t_electron_confinement = self.hubbard_upper_confinement_time( + t_global_confinement_scaling = self.hubbard_upper_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla20, @@ -854,7 +858,7 @@ def calculate_confinement_time( # Menard NSTX, ELMy H-mode scaling elif model == ConfinementTimeModel.MENARD_NSTX: - t_electron_confinement = self.menard_nstx_confinement_time( + t_global_confinement_scaling = self.menard_nstx_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -869,22 +873,24 @@ def calculate_confinement_time( # Menard NSTX-Petty08 Hybrid elif model == ConfinementTimeModel.MENARD_NSTX_PETTY08_HYBRID: - t_electron_confinement = self.menard_nstx_petty08_hybrid_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - self.data.physics.kappa_ipb, - aspect, - m_fuel_amu, + t_global_confinement_scaling = ( + self.menard_nstx_petty08_hybrid_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + self.data.physics.kappa_ipb, + aspect, + m_fuel_amu, + ) ) # ========================================================================== # NSTX gyro-Bohm (Buxton) elif model == ConfinementTimeModel.NSTX_GYRO_BOHM: - t_electron_confinement = self.nstx_gyro_bohm_confinement_time( + t_global_confinement_scaling = self.nstx_gyro_bohm_confinement_time( pcur, b_plasma_toroidal_on_axis, p_plasma_loss_mw, @@ -896,7 +902,7 @@ def calculate_confinement_time( # ITPA20 H-mode scaling elif model == ConfinementTimeModel.ITPA20: - t_electron_confinement = self.itpa20_confinement_time( + t_global_confinement_scaling = self.itpa20_confinement_time( pcur, b_plasma_toroidal_on_axis, dnla19, @@ -912,7 +918,7 @@ def calculate_confinement_time( # ITPA20-IL confinement time scaling elif model == ConfinementTimeModel.ITPA20_IL: - t_electron_confinement = self.itpa20_il_confinement_time( + t_global_confinement_scaling = self.itpa20_il_confinement_time( pcur, b_plasma_toroidal_on_axis, p_plasma_loss_mw, @@ -931,12 +937,6 @@ def calculate_confinement_time( i_confinement_time=i_confinement_time, ) - # Apply H-factor correction to chosen scaling - t_electron_energy_confinement = hfact * t_electron_confinement - - # Ion energy confinement time - t_ion_energy_confinement = t_electron_energy_confinement - # Calculate H* non-radiation corrected H factor # Note: we will assume the IPB-98y2 scaling. if ( @@ -972,11 +972,22 @@ def calculate_confinement_time( elif self.data.physics.i_rad_loss == ConfinementRadiationLossModel.NO_RADIATION: hstar = hfact - # Calculation of the transport power loss terms - # Transport losses in Watts/m3 are 3/2 * n.e.T / tau , with T in eV - # (here, temp_plasma_ion_density_weighted_kev and - # temp_plasma_electron_density_weighted_kev are in keV, and - # pden_electron_transport_loss_mw and pden_ion_transport_loss_mw are in MW/m3) + t_energy_confinement = t_global_confinement_scaling * hfact + + # Apply H-factor correction to chosen scaling + t_electron_energy_confinement = ( + self.calculate_electron_species_consistent_energy_confinement_time( + self.data.physics.f_t_fuel_ion_electron_energy_confinement, + eden_plasma_ions_thermal_vol_avg, + eden_plasma_electrons_thermal_vol_avg, + t_energy_confinement, + ) + ) + # Ion energy confinement time + t_ion_energy_confinement = ( + t_electron_energy_confinement + * self.data.physics.f_t_fuel_ion_electron_energy_confinement + ) # The transport losses is just the electron and ion thermal energies divided by # the confinement time. @@ -987,14 +998,6 @@ def calculate_confinement_time( eden_plasma_electrons_thermal_vol_avg / t_electron_energy_confinement ) / 1e6 # Convert from W/m³ to MW/m³ - ratio = eden_plasma_ions_thermal_vol_avg / eden_plasma_electrons_thermal_vol_avg - - # Global energy confinement time - - t_energy_confinement = (ratio + 1.0e0) / ( - ratio / t_ion_energy_confinement + 1.0e0 / t_electron_energy_confinement - ) - # For comparison directly calculate the confinement time from the stored energy # calculated from the total plasma beta and the loss power used above. self.data.physics.t_energy_confinement_beta = ( @@ -1041,6 +1044,48 @@ def calculate_double_and_triple_product( return ntau, nTtau + @staticmethod + def calculate_electron_species_consistent_energy_confinement_time( + f_t_fuel_ion_electron_energy_confinement: float, + eden_plasma_ions_thermal_vol_avg: float, + eden_plasma_electrons_thermal_vol_avg: float, + t_plasma_global_energy_confinement: float, + ) -> float: + """Calculate the electron energy confinement time based on the global plasma + energy confinement time and the species ratio. + + Parameters + ---------- + f_t_fuel_ion_electron_energy_confinement : float + Ratio of fuel ion to electron energy confinement times. + eden_plasma_ions_thermal_vol_avg : float + Volume averaged thermal energy density of ions [J/m³]. + eden_plasma_electrons_thermal_vol_avg : float + Volume averaged thermal energy density of electrons [J/m³]. + t_plasma_global_energy_confinement : float + Global plasma energy confinement time [s]. + + Returns + ------- + float + Electron energy confinement time [s]. + + Notes + ----- + - The variable `f_t_fuel_ion_electron_energy_confinement` is used to relate + the energy confinement times of fuel ions and electrons. R = τ_i / τ_e. + """ + return t_plasma_global_energy_confinement * ( + ( + eden_plasma_electrons_thermal_vol_avg + + ( + eden_plasma_ions_thermal_vol_avg + / f_t_fuel_ion_electron_energy_confinement + ) + ) + / (eden_plasma_electrons_thermal_vol_avg + eden_plasma_ions_thermal_vol_avg) + ) + def find_other_h_factors(self, i_confinement_time: int) -> float: """Function to find H-factor for the equivalent confinement time in other scalings. @@ -1194,6 +1239,13 @@ def output_confinement_time_info(self): po.ocmmnt(self.outfile, "----------------------------") po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Ratio of fuel ion to electron energy confinement times", + "(f_t_fuel_ion_electron_energy_confinement)", + self.data.physics.f_t_fuel_ion_electron_energy_confinement, + "OP ", + ) po.ovarre( self.outfile, "Ion energy confinement time, from scaling (s)", diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 53d44266ff..995ec44a59 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -796,7 +796,7 @@ def bosch_hale_reactivity( def set_fusion_powers( f_alpha_electron: float, - f_alpha_ion: float, + f_p_alpha_total_ions: float, p_beam_alpha_mw: float, pden_non_alpha_charged_mw: float, pden_plasma_neutron_mw: float, @@ -811,7 +811,7 @@ def set_fusion_powers( ---------- f_alpha_electron : float - f_alpha_ion : + f_p_alpha_total_ions : float p_beam_alpha_mw : float @@ -840,9 +840,9 @@ def set_fusion_powers( - p_non_alpha_charged_mw (float): Other total charged particle fusion power [MW]. - pden_alpha_total_mw (float): Alpha power per unit volume, from beams and plasma [MW/m³]. - - f_pden_alpha_electron_mw (float): Alpha power per unit volume to + - pden_alpha_heating_electrons_mw (float): Alpha power per unit volume to electrons [MW/m³]. - - f_pden_alpha_ions_mw (float): Alpha power per unit volume to ions [MW/m³]. + - pden_alpha_heating_ions_mw (float): Alpha power per unit volume to ions [MW/m³]. - p_charged_particle_mw (float): Charged particle fusion power [MW]. - p_fusion_total_mw (float): Total fusion power [MW]. @@ -896,8 +896,10 @@ def set_fusion_powers( # Alpha power to electrons and ions (used with electron # and ion power balance equations only) # No consideration of pden_non_alpha_charged_mw here. - f_pden_alpha_ions_mw = f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_alpha_ion - f_pden_alpha_electron_mw = ( + pden_alpha_heating_ions_mw = ( + f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_p_alpha_total_ions + ) + pden_alpha_heating_electrons_mw = ( f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_alpha_electron ) @@ -909,8 +911,8 @@ def set_fusion_powers( p_neutron_total_mw, p_non_alpha_charged_mw, pden_alpha_total_mw, - f_pden_alpha_electron_mw, - f_pden_alpha_ions_mw, + pden_alpha_heating_electrons_mw, + pden_alpha_heating_ions_mw, p_charged_particle_mw, p_fusion_total_mw, ) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 413e545d4c..f6a4046f0f 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -98,57 +98,42 @@ def calculate_cylindrical_safety_factor( @nb.jit(nopython=True, cache=True) -def rether( - alphan, - alphat, - nd_plasma_electrons_vol_avg, - dlamie, - te, - temp_plasma_ion_vol_avg_kev, - n_charge_plasma_effective_mass_weighted_vol_avg, -): +def calculate_ion_electron_equilibration_power( + nd_electrons: float | np.ndarray, + temp_plasma_electron_kev: float | np.ndarray, + temp_plasma_ion_kev: float | np.ndarray, + t_plasma_ion_electron_equilibration: float | np.ndarray, +) -> float | np.ndarray: """Routine to find the equilibration power between the - ions and electrons - This routine calculates the equilibration power between the - ions and electrons. - Unknown origin + ions and electrons (Pₑᵢ) in the plasma. Parameters ---------- - alphan : - density profile index - alphat : - temperature profile index - nd_plasma_electrons_vol_avg : - electron density (/m3) - dlamie : - ion-electron coulomb logarithm - te : - electron temperature (keV) - temp_plasma_ion_vol_avg_kev : - ion temperature (keV) - n_charge_plasma_effective_mass_weighted_vol_avg : - mass weighted plasma effective charge + nd_electrons : float | np.ndarray + Electron density [/m³] + temp_plasma_electron_kev : float | np.ndarray + Electron temperature [keV] + temp_plasma_ion_kev : float | np.ndarray + Ion temperature [keV] + t_plasma_ion_electron_equilibration : float | np.ndarray + Ion-electron equilibration time [s] Returns ------- - pden_ion_electron_equilibration_mw : - ion/electron equilibration power (MW/m3) + pden_ion_electron_equilibration : + ion/electron equilibration power (Pₑᵢ) [W/m³] """ - profie = (1.0 + alphan) ** 2 / ( - (2.0 * alphan - 0.5 * alphat + 1.0) * np.sqrt(1.0 + alphat) - ) - conie = ( - 2.42165e-41 - * dlamie - * nd_plasma_electrons_vol_avg**2 - * n_charge_plasma_effective_mass_weighted_vol_avg - * profie + return ( + (3 / 2) + * nd_electrons + * constants.KILOELECTRON_VOLT + * ( + (temp_plasma_electron_kev - temp_plasma_ion_kev) + / t_plasma_ion_electron_equilibration + ) ) - return conie * (temp_plasma_ion_vol_avg_kev - te) / (te**1.5) - # ----------------------------------------------------- # Diamagnetic and Pfirsch-Schlüter Current Calculations @@ -676,13 +661,13 @@ def run(self): self.data.physics.p_neutron_total_mw, self.data.physics.p_non_alpha_charged_mw, self.data.physics.pden_alpha_total_mw, - self.data.physics.f_pden_alpha_electron_mw, - self.data.physics.f_pden_alpha_ions_mw, + self.data.physics.pden_alpha_heating_electrons_mw, + self.data.physics.pden_alpha_heating_ions_mw, self.data.physics.p_charged_particle_mw, self.data.physics.p_fusion_total_mw, ) = reactions.set_fusion_powers( - self.data.physics.f_alpha_electron, - self.data.physics.f_alpha_ion, + self.data.physics.f_p_alpha_total_electron, + self.data.physics.f_p_alpha_total_ions, self.data.physics.p_beam_alpha_mw, self.data.physics.pden_non_alpha_charged_mw, self.data.physics.pden_plasma_neutron_mw, @@ -691,6 +676,15 @@ def run(self): self.data.physics.f_p_alpha_plasma_deposited, ) + self.data.physics.p_alpha_heating_electrons_mw = ( + self.data.physics.pden_alpha_heating_electrons_mw + * self.data.physics.vol_plasma + ) + + self.data.physics.p_alpha_heating_ions_mw = ( + self.data.physics.pden_alpha_heating_ions_mw * self.data.physics.vol_plasma + ) + self.data.physics.beta_fast_alpha = self.beta.fast_alpha_beta( b_plasma_poloidal_average=self.data.physics.b_plasma_surface_poloidal_average, b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis, @@ -707,14 +701,25 @@ def run(self): # Calculate ion/electron equilibration power - self.data.physics.pden_ion_electron_equilibration_mw = rether( - self.data.physics.alphan, - self.data.physics.alphat, - self.data.physics.nd_plasma_electrons_vol_avg, - self.data.physics.dlamie, - self.data.physics.temp_plasma_electron_vol_avg_kev, - self.data.physics.temp_plasma_ion_vol_avg_kev, - self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg, + pden_ion_electron_equilibration_vol_avg = calculate_ion_electron_equilibration_power( + nd_electrons=self.data.physics.nd_plasma_electrons_vol_avg, + temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_density_weighted_kev, + temp_plasma_ion_kev=self.data.physics.temp_plasma_ion_density_weighted_kev, + t_plasma_ion_electron_equilibration=1.0 + / ( + 1.0 / self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg + + 1.0 / self.data.physics.t_plasma_electron_triton_equilibration_vol_avg + + 1.0 + / self.data.physics.t_plasma_electron_alpha_thermal_equilibration_vol_avg + ), + ) + self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = ( + pden_ion_electron_equilibration_vol_avg / 1.0e6 + ) + + self.data.physics.p_ion_electron_equilibration_vol_avg_mw = ( + self.data.physics.pden_ion_electron_equilibration_vol_avg_mw + * self.data.physics.vol_plasma ) # Calculate radiation power @@ -1359,10 +1364,12 @@ def plasma_composition(self): else: pc = self.data.physics.f_temp_plasma_electron_density_vol_avg - self.data.physics.f_alpha_electron = 0.88155 * np.exp( + self.data.physics.f_p_alpha_total_electron = 0.88155 * np.exp( -self.data.physics.temp_plasma_electron_vol_avg_kev * pc / 67.4036 ) - self.data.physics.f_alpha_ion = 1.0 - self.data.physics.f_alpha_electron + self.data.physics.f_p_alpha_total_ions = ( + 1.0 - self.data.physics.f_p_alpha_total_electron + ) # ====================================================================== @@ -2004,15 +2011,29 @@ def outplas(self): po.ovarre( self.outfile, "Alpha power per unit volume transferred to electrons (MW/m³)", - "(f_pden_alpha_electron_mw)", - self.data.physics.f_pden_alpha_electron_mw, + "(pden_alpha_heating_electrons_mw)", + self.data.physics.pden_alpha_heating_electrons_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Alpha power transferred to electrons [MW]", + "(p_alpha_heating_electrons_mw)", + self.data.physics.p_alpha_heating_electrons_mw, "OP ", ) po.ovarre( self.outfile, "Alpha power per unit volume transferred to ions (MW/m³)", - "(f_pden_alpha_ions_mw)", - self.data.physics.f_pden_alpha_ions_mw, + "(pden_alpha_heating_ions_mw)", + self.data.physics.pden_alpha_heating_ions_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Alpha power transferred to ions [MW]", + "(p_alpha_heating_ions_mw)", + self.data.physics.p_alpha_heating_ions_mw, "OP ", ) @@ -2271,73 +2292,10 @@ def outplas(self): ) po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Total heating power given to the plasma (Pₕₑₐₜ) [MW]", - "(p_plasma_heating_total_mw)", - self.data.physics.p_plasma_heating_total_mw, - "OP ", - ) - po.ovarre( - self.outfile, - "Ohmic heating power (MW)", - "(p_plasma_ohmic_mw)", - self.data.physics.p_plasma_ohmic_mw, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of alpha power deposited in plasma", - "(f_p_alpha_plasma_deposited)", - self.data.physics.f_p_alpha_plasma_deposited, - "IP", - ) - po.ovarre( - self.outfile, - "Fraction of alpha power to electrons", - "(f_alpha_electron)", - self.data.physics.f_alpha_electron, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of alpha power to ions", - "(f_alpha_ion)", - self.data.physics.f_alpha_ion, - "OP ", - ) - po.ovarre( - self.outfile, - "Ion transport (MW)", - "(p_ion_transport_loss_mw)", - self.data.physics.p_ion_transport_loss_mw, - "OP ", - ) - po.ovarre( - self.outfile, - "Electron transport (MW)", - "(p_electron_transport_loss_mw)", - self.data.physics.p_electron_transport_loss_mw, - "OP ", - ) - po.ovarre( - self.outfile, - "Injection power to ions (MW)", - "(p_hcd_injected_ions_mw)", - self.data.current_drive.p_hcd_injected_ions_mw, - "OP ", - ) - po.ovarre( - self.outfile, - "Injection power to electrons (MW)", - "(p_hcd_injected_electrons_mw)", - self.data.current_drive.p_hcd_injected_electrons_mw, - "OP ", - ) - if ( - PlasmaIgnitionModel(self.data.physics.i_plasma_ignited) - == PlasmaIgnitionModel.IGNITED - ): + + self.output_plasma_power_balance() + + if self.data.physics.i_plasma_ignited == 1: po.ocmmnt(self.outfile, " (Injected power only used for start-up phase)") self.exhaust.output() @@ -2441,6 +2399,203 @@ def outplas(self): self.data.reinke.fzactual, ) + def output_plasma_power_balance(self) -> None: + """Output information about plasma power balance.""" + po.oheadr(self.outfile, "Plasma Power Balance") + + po.ocmmnt(self.outfile, "Global power balance") + po.ocmmnt(self.outfile, "----------------------------") + po.oblnkl(self.outfile) + + po.ocmmnt(self.outfile, "Heating powers:") + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Ohmic heating power (MW)", + "(p_plasma_ohmic_mw)", + self.data.physics.p_plasma_ohmic_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Total alpha heating power (MW)", + "(p_alpha_total_mw)", + self.data.physics.p_alpha_total_mw, + "OP ", + ) + po.ovarrf( + self.outfile, + "Fraction of alpha power deposited in plasma", + "(f_p_alpha_plasma_deposited)", + self.data.physics.f_p_alpha_plasma_deposited, + "IP", + ) + po.ovarrf( + self.outfile, + "Total power injected into the plasma (MW)", + "(p_hcd_injected_total_mw)", + self.data.current_drive.p_hcd_injected_total_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Non-alpha heating power (MW)", + "(p_non_alpha_charged_mw)", + self.data.physics.p_non_alpha_charged_mw, + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Total heating power given to the plasma (Pₕₑₐₜ) [MW]", + "(p_plasma_heating_total_mw)", + self.data.physics.p_plasma_heating_total_mw, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "Loss powers:") + po.oblnkl(self.outfile) + + po.ovarre( + self.outfile, + "Total transport loss power (Pₗ) [MW]", + "(p_plasma_loss_mw)", + self.data.physics.p_plasma_loss_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Total radiated power from the plasma (Pᵧ) [MW]", + "(p_plasma_rad_mw)", + self.data.physics.p_plasma_rad_mw, + "OP ", + ) + + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + po.osubhd(self.outfile, "Electron power balance :") + po.ocmmnt(self.outfile, "Heating powers:") + po.oblnkl(self.outfile) + + po.ovarre( + self.outfile, + "Ohmic heating power to electrons (MW)", + "(p_plasma_ohmic_mw)", + self.data.physics.p_plasma_ohmic_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Alpha heating power to electrons (MW)", + "(p_alpha_heating_electrons_mw)", + self.data.physics.p_alpha_heating_electrons_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Injection power to electrons (MW)", + "(p_hcd_injected_electrons_mw)", + self.data.current_drive.p_hcd_injected_electrons_mw, + "OP ", + ) + + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "Loss powers:") + po.oblnkl(self.outfile) + + po.ovarre( + self.outfile, + "Electron transport (MW)", + "(p_electron_transport_loss_mw)", + self.data.physics.p_electron_transport_loss_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Ion-electron equilibration power (MW)", + "(p_ion_electron_equilibration_vol_avg_mw)", + self.data.physics.p_ion_electron_equilibration_vol_avg_mw, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma radiation power (MW)", + "(p_plasma_rad_mw)", + self.data.physics.p_plasma_rad_mw, + "OP ", + ) + po.oblnkl(self.outfile) + p_electron_balance = ( + self.data.physics.p_plasma_ohmic_mw + + self.data.physics.p_alpha_heating_electrons_mw + + self.data.current_drive.p_hcd_injected_electrons_mw + - self.data.physics.p_electron_transport_loss_mw + - self.data.physics.p_ion_electron_equilibration_vol_avg_mw + - self.data.physics.p_plasma_rad_mw + ) + po.ocmmnt(self.outfile, "----------------------------") + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Electron power balance (MW)", + "(p_electron_balance_mw)", + p_electron_balance, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + po.osubhd(self.outfile, "Ions power balance :") + + po.ocmmnt(self.outfile, "Heating powers:") + po.oblnkl(self.outfile) + po.ovarrf( + self.outfile, + "Fraction of alpha power to ions", + "(f_p_alpha_total_ions)", + self.data.physics.f_p_alpha_total_ions, + "OP ", + ) + po.ovarre( + self.outfile, + "Alpha heating power to ions (MW)", + "(p_alpha_heating_ions_mw)", + self.data.physics.p_alpha_heating_ions_mw, + "OP ", + ) + + po.ovarre( + self.outfile, + "Injection power to ions (MW)", + "(p_hcd_injected_ions_mw)", + self.data.current_drive.p_hcd_injected_ions_mw, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "Loss powers:") + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Ion transport (MW)", + "(p_ion_transport_loss_mw)", + self.data.physics.p_ion_transport_loss_mw, + "OP ", + ) + p_ion_balance = ( + self.data.physics.p_alpha_heating_ions_mw + + self.data.current_drive.p_hcd_injected_ions_mw + - self.data.physics.p_ion_transport_loss_mw + - self.data.physics.p_ion_electron_equilibration_vol_avg_mw + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + po.ovarre( + self.outfile, + "Ion power balance (MW)", + "(p_ion_balance_mw)", + p_ion_balance, + "OP ", + ) + def output_temperature_density_profile_info(self) -> None: """Output information about plasma temperature and density profiles.""" po.oheadr(self.outfile, "Plasma Density and Temperature Profiles") @@ -5495,6 +5650,100 @@ def run(self): temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, ) + # ================================ + # Ion-electron equilibration times + # ================================ + + self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg = self.calculate_equilibriation_time( + temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev, + temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev + * self.data.physics.f_temp_plasma_ion_electron, + nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg + * self.data.physics.f_plasma_fuel_deuterium + * ( + self.data.physics.nd_plasma_fuel_ions_vol_avg + / self.data.physics.nd_plasma_electrons_vol_avg + ), + plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_deuteron_vol_avg, + m_ion=constants.DEUTERON_MASS, + n_charge_ion=1, + ) + + self.data.physics.t_plasma_electron_deuteron_equilibration_profile = self.calculate_equilibriation_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y + * self.data.physics.f_temp_plasma_ion_electron, + nd_plasma_ions=( + self.plasma_profile.neprofile.profile_y + * self.data.physics.f_plasma_fuel_deuterium + * ( + self.data.physics.nd_plasma_fuel_ions_vol_avg + / self.data.physics.nd_plasma_electrons_vol_avg + ) + ), + plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_deuteron_profile, + m_ion=constants.DEUTERON_MASS, + n_charge_ion=1, + ) + + self.data.physics.t_plasma_electron_triton_equilibration_vol_avg = self.calculate_equilibriation_time( + temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev, + temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev + * self.data.physics.f_temp_plasma_ion_electron, + nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg + * self.data.physics.f_plasma_fuel_tritium + * ( + self.data.physics.nd_plasma_fuel_ions_vol_avg + / self.data.physics.nd_plasma_electrons_vol_avg + ), + plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_triton_vol_avg, + m_ion=constants.TRITON_MASS, + n_charge_ion=1, + ) + + self.data.physics.t_plasma_electron_triton_equilibration_profile = self.calculate_equilibriation_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y + * self.data.physics.f_temp_plasma_ion_electron, + nd_plasma_ions=( + self.plasma_profile.neprofile.profile_y + * self.data.physics.f_plasma_fuel_tritium + * ( + self.data.physics.nd_plasma_fuel_ions_vol_avg + / self.data.physics.nd_plasma_electrons_vol_avg + ) + ), + plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_triton_profile, + m_ion=constants.TRITON_MASS, + n_charge_ion=1, + ) + + self.data.physics.t_plasma_electron_alpha_thermal_equilibration_vol_avg = self.calculate_equilibriation_time( + temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev, + temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev + * self.data.physics.f_temp_plasma_ion_electron, + nd_plasma_ions=self.data.physics.nd_plasma_alphas_thermal_vol_avg, + plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_alpha_thermal_vol_avg, + m_ion=constants.ALPHA_MASS, + n_charge_ion=2, + ) + + self.data.physics.t_plasma_electron_alpha_thermal_equilibration_profile = self.calculate_equilibriation_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y + * self.data.physics.f_temp_plasma_ion_electron, + nd_plasma_ions=( + self.plasma_profile.neprofile.profile_y + * ( + self.data.physics.nd_plasma_alphas_thermal_vol_avg + / self.data.physics.nd_plasma_electrons_vol_avg + ) + ), + plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_alpha_thermal_profile, + m_ion=constants.ALPHA_MASS, + n_charge_ion=2, + ) + @staticmethod @nb.njit(cache=True) def calculate_debye_length( @@ -5937,6 +6186,67 @@ def calculate_spitzer_resistivity( / (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT) ** 1.5 ) + @staticmethod + @nb.njit(cache=True) + def calculate_equilibriation_time( + temp_plasma_electron_kev: float | np.ndarray, + temp_plasma_ion_kev: float | np.ndarray, + nd_plasma_ions: float | np.ndarray, + plasma_coulomb_log_electron_ion: float | np.ndarray, + m_ion: float, + n_charge_ion: float = 1.0, + ): + """ + Calculate the equilibration time (τ_eq) between electrons and ions in a plasma. + + Parameters + ---------- + temp_plasma_electron_kev : float | np.ndarray + Electron temperature in keV. + temp_plasma_ion_kev : float | np.ndarray + Ion temperature in keV. + nd_plasma_ions : float | np.ndarray + Ion density (m^-3). + plasma_coulomb_log_electron_ion : float | np.ndarray + Coulomb logarithm for electron-ion collisions. + m_ion : float + Ion mass (kg). + n_charge_ion : float, optional + Charge number (Z) of the ion. Default is 1.0. + + Returns + ------- + float | np.ndarray + Equilibration time (s). + + Notes + ----- + - The equilibration time is the characteristic time for energy exchange between + electrons and ions in a plasma, leading to thermal equilibrium. + - It is the characteristic timescale required for two different particle species + to share heat and reach thermal equilibrium. + """ + return ( + 3 + * (4 * np.pi * constants.EPSILON0) ** 2.0 + * m_ion + * (constants.ELECTRON_MASS) + * ( + ( + (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT) + / constants.ELECTRON_MASS + ) + + ((temp_plasma_ion_kev * constants.KILOELECTRON_VOLT) / m_ion) + ) + ** 1.5 + ) / ( + (8 * np.sqrt(2 * np.pi)) + * nd_plasma_ions + * n_charge_ion**2 + * plasma_coulomb_log_electron_ion + * constants.ELECTRON_CHARGE**4 + ) + def output_detailed_physics(self): """Outputs detailed physics variables to file.""" po.oheadr(self.outfile, "Detailed Plasma") @@ -6495,3 +6805,58 @@ def output_detailed_physics(self): f"(res_plasma_fuel_spitzer_profile{i})", self.data.physics.res_plasma_fuel_spitzer_profile[i], ) + + po.osubhd(self.outfile, "Equilibration Times:") + + po.ovarre( + self.outfile, + "Volume averaged electron-deuteron equilibration time (τ_eq) (s)", + "(t_plasma_electron_deuteron_equilibration_vol_avg)", + self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg, + ) + + for i in range( + len(self.data.physics.t_plasma_electron_deuteron_equilibration_profile) + ): + po.ovarre( + self.mfile, + f"Electron-deuteron equilibration time at point {i}", + f"(t_plasma_electron_deuteron_equilibration_profile{i})", + self.data.physics.t_plasma_electron_deuteron_equilibration_profile[i], + ) + + po.ovarre( + self.outfile, + "Volume averaged electron-triton equilibration time (τ_eq) (s)", + "(t_plasma_electron_triton_equilibration_vol_avg)", + self.data.physics.t_plasma_electron_triton_equilibration_vol_avg, + ) + + for i in range( + len(self.data.physics.t_plasma_electron_triton_equilibration_profile) + ): + po.ovarre( + self.mfile, + f"Electron-triton equilibration time at point {i}", + f"(t_plasma_electron_triton_equilibration_profile{i})", + self.data.physics.t_plasma_electron_triton_equilibration_profile[i], + ) + + po.ovarre( + self.outfile, + "Volume averaged electron-alpha equilibration time (τ_eq) (s)", + "(t_plasma_electron_alpha_thermal_equilibration_vol_avg)", + self.data.physics.t_plasma_electron_alpha_thermal_equilibration_vol_avg, + ) + + for i in range( + len(self.data.physics.t_plasma_electron_alpha_thermal_equilibration_profile) + ): + po.ovarre( + self.mfile, + f"Electron-alpha equilibration time at point {i}", + f"(t_plasma_electron_alpha_thermal_equilibration_profile{i})", + self.data.physics.t_plasma_electron_alpha_thermal_equilibration_profile[ + i + ], + ) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index 4d23487128..b99325ed68 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -15,7 +15,10 @@ from process.core.model import Model from process.data_structure.physics_variables import PlasmaIgnitionModel from process.models.engineering.pumping import CoolantType -from process.models.physics.physics import Physics, rether +from process.models.physics.physics import ( + Physics, + calculate_ion_electron_equilibration_power, +) from process.models.power import PumpingPowerModelTypes from process.models.stellarator.build import st_build from process.models.stellarator.coils.calculate import st_coil @@ -2047,13 +2050,13 @@ def st_phys(self, output): self.data.physics.p_neutron_total_mw, self.data.physics.p_non_alpha_charged_mw, self.data.physics.pden_alpha_total_mw, - self.data.physics.f_pden_alpha_electron_mw, - self.data.physics.f_pden_alpha_ions_mw, + self.data.physics.pden_alpha_heating_electrons_mw, + self.data.physics.pden_alpha_heating_ions_mw, self.data.physics.p_charged_particle_mw, self.data.physics.p_fusion_total_mw, ) = reactions.set_fusion_powers( - self.data.physics.f_alpha_electron, - self.data.physics.f_alpha_ion, + self.data.physics.f_p_alpha_total_electron, + self.data.physics.f_p_alpha_total_ions, self.data.physics.p_beam_alpha_mw, self.data.physics.pden_non_alpha_charged_mw, self.data.physics.pden_plasma_neutron_mw, @@ -2104,14 +2107,16 @@ def st_phys(self, output): # Calculate ion/electron equilibration power - self.data.physics.pden_ion_electron_equilibration_mw = rether( - self.data.physics.alphan, - self.data.physics.alphat, - self.data.physics.nd_plasma_electrons_vol_avg, - self.data.physics.dlamie, - self.data.physics.temp_plasma_electron_vol_avg_kev, - self.data.physics.temp_plasma_ion_vol_avg_kev, - self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg, + self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = ( + calculate_ion_electron_equilibration_power( + self.data.physics.alphan, + self.data.physics.alphat, + self.data.physics.nd_plasma_electrons_vol_avg, + self.data.physics.dlamie, + self.data.physics.temp_plasma_electron_vol_avg_kev, + self.data.physics.temp_plasma_ion_vol_avg_kev, + self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg, + ) ) # Calculate radiation power diff --git a/tests/unit/models/physics/test_fusion_reactions.py b/tests/unit/models/physics/test_fusion_reactions.py index b64821e975..9bbc1f1a3d 100644 --- a/tests/unit/models/physics/test_fusion_reactions.py +++ b/tests/unit/models/physics/test_fusion_reactions.py @@ -12,9 +12,9 @@ class SetFusionPowersParam(NamedTuple): f_p_alpha_plasma_deposited: Any = None - f_alpha_electron: Any = None + f_p_alpha_total_electron: Any = None - f_alpha_ion: Any = None + f_p_alpha_total_ions: Any = None p_beam_alpha_mw: Any = None @@ -50,8 +50,8 @@ class SetFusionPowersParam(NamedTuple): [ SetFusionPowersParam( f_p_alpha_plasma_deposited=0.95, - f_alpha_electron=0.68, - f_alpha_ion=0.32, + f_p_alpha_total_electron=0.68, + f_p_alpha_total_ions=0.32, p_beam_alpha_mw=0, pden_non_alpha_charged_mw=0.00066, vol_plasma=2426.25, @@ -69,8 +69,8 @@ class SetFusionPowersParam(NamedTuple): ), SetFusionPowersParam( f_p_alpha_plasma_deposited=0.95, - f_alpha_electron=0.68, - f_alpha_ion=0.32, + f_p_alpha_total_electron=0.68, + f_p_alpha_total_ions=0.32, p_beam_alpha_mw=100.5, pden_non_alpha_charged_mw=0.00066, vol_plasma=2426.25, @@ -110,13 +110,13 @@ def test_set_fusion_powers(setfusionpowersparam): p_neutron_total_mw, p_non_alpha_charged_mw, pden_alpha_total_mw, - f_pden_alpha_electron_mw, - f_pden_alpha_ions_mw, + pden_alpha_heating_electrons_mw, + pden_alpha_heating_ions_mw, p_charged_particle_mw, p_fusion_total_mw, ) = reactions.set_fusion_powers( - f_alpha_electron=setfusionpowersparam.f_alpha_electron, - f_alpha_ion=setfusionpowersparam.f_alpha_ion, + f_p_alpha_total_electron=setfusionpowersparam.f_p_alpha_total_electron, + f_p_alpha_total_ions=setfusionpowersparam.f_p_alpha_total_ions, p_beam_alpha_mw=setfusionpowersparam.p_beam_alpha_mw, pden_non_alpha_charged_mw=setfusionpowersparam.pden_non_alpha_charged_mw, pden_plasma_neutron_mw=setfusionpowersparam.pden_plasma_neutron_mw, @@ -140,10 +140,10 @@ def test_set_fusion_powers(setfusionpowersparam): assert p_non_alpha_charged_mw == pytest.approx( setfusionpowersparam.expected_non_alpha_charged_power ) - assert f_pden_alpha_electron_mw == pytest.approx( + assert pden_alpha_heating_electrons_mw == pytest.approx( setfusionpowersparam.expected_alpha_power_electron_density ) - assert f_pden_alpha_ions_mw == pytest.approx( + assert pden_alpha_heating_ions_mw == pytest.approx( setfusionpowersparam.expected_alpha_power_ion_density ) assert p_charged_particle_mw == pytest.approx( diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index ccef3e81db..a241a81b56 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -14,9 +14,9 @@ PlasmaBeta, PlasmaInductance, calculate_cylindrical_safety_factor, + calculate_ion_electron_equilibration_power, ps_fraction_scene, res_diff_time, - rether, ) @@ -1235,7 +1235,7 @@ class PlasmaCompositionParam(NamedTuple): i_plasma_ignited: Any = None - f_alpha_electron: Any = None + f_p_alpha_total_electron: Any = None m_fuel_amu: Any = None @@ -1255,7 +1255,7 @@ class PlasmaCompositionParam(NamedTuple): f_nd_plasma_oxygen_electron: Any = None - f_alpha_ion: Any = None + f_p_alpha_total_ions: Any = None f_nd_alpha_thermal_electron: Any = None @@ -1379,7 +1379,7 @@ class PlasmaCompositionParam(NamedTuple): ]), alphat=1.45, i_plasma_ignited=PlasmaIgnitionModel.NON_IGNITED, - f_alpha_electron=0, + f_p_alpha_total_electron=0, m_fuel_amu=0, f_plasma_fuel_tritium=0.5, nd_plasma_fuel_ions_vol_avg=0, @@ -1389,7 +1389,7 @@ class PlasmaCompositionParam(NamedTuple): n_charge_plasma_effective_mass_weighted_vol_avg=0, f_nd_plasma_carbon_electron=0, f_nd_plasma_oxygen_electron=0, - f_alpha_ion=0, + f_p_alpha_total_ions=0, f_nd_alpha_thermal_electron=0.10000000000000001, f_nd_beam_electron=0, n_charge_plasma_effective_vol_avg=0, @@ -1485,7 +1485,7 @@ class PlasmaCompositionParam(NamedTuple): ).transpose(), alphat=1.45, i_plasma_ignited=PlasmaIgnitionModel.NON_IGNITED, - f_alpha_electron=0.6845930883190634, + f_p_alpha_total_electron=0.6845930883190634, m_fuel_amu=2.5, f_plasma_fuel_tritium=0.5, nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19, @@ -1495,7 +1495,7 @@ class PlasmaCompositionParam(NamedTuple): n_charge_plasma_effective_mass_weighted_vol_avg=0.43046641789338563, f_nd_plasma_carbon_electron=0, f_nd_plasma_oxygen_electron=0, - f_alpha_ion=0.3154069116809366, + f_p_alpha_total_ions=0.3154069116809366, f_nd_alpha_thermal_electron=0.10000000000000001, f_nd_beam_electron=0, n_charge_plasma_effective_vol_avg=2.0909945616489103, @@ -1581,7 +1581,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics): for field in [ "alphat", "i_plasma_ignited", - "f_alpha_electron", + "f_p_alpha_total_electron", "m_fuel_amu", "f_plasma_fuel_tritium", "nd_plasma_fuel_ions_vol_avg", @@ -1591,7 +1591,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics): "n_charge_plasma_effective_mass_weighted_vol_avg", "f_nd_plasma_carbon_electron", "f_nd_plasma_oxygen_electron", - "f_alpha_ion", + "f_p_alpha_total_ions", "f_nd_alpha_thermal_electron", "f_nd_beam_electron", "n_charge_plasma_effective_vol_avg", @@ -1621,7 +1621,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics): plasmacompositionparam.expected_impurity_arr_frac ) - assert physics.data.physics.f_alpha_electron == pytest.approx( + assert physics.data.physics.f_p_alpha_total_electron == pytest.approx( plasmacompositionparam.expected_f_alpha_electron ) @@ -1646,7 +1646,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics): == pytest.approx(plasmacompositionparam.expected_zeffai) ) - assert physics.data.physics.f_alpha_ion == pytest.approx( + assert physics.data.physics.f_p_alpha_total_ions == pytest.approx( plasmacompositionparam.expected_f_alpha_ion ) @@ -1970,10 +1970,13 @@ def test_phyaux(phyauxparam, monkeypatch, physics): assert t_alpha_confinement == pytest.approx(phyauxparam.expected_t_alpha_confinement) -def test_rether(): - assert rether( - 1.0, 1.45, 7.5e19, 17.81065204, 12.0, 13.0, 0.43258985 - ) == pytest.approx(0.028360489673597476) +def test_calculate_ion_electron_equilibration_power(): + assert calculate_ion_electron_equilibration_power( + nd_electrons=7.5e19, + temp_plasma_electron_kev=17.81065204, + temp_plasma_ion_kev=12.0, + t_plasma_ion_electron_equilibration=0.43258985, + ) == pytest.approx(242109.29342520377) class PohmParam(NamedTuple):