diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md index 6032e13737..d5a772bb23 100644 --- a/documentation/source/physics-models/profiles/plasma_profiles.md +++ b/documentation/source/physics-models/profiles/plasma_profiles.md @@ -636,6 +636,47 @@ The same function is run from the `i_plasma_pedestal == 0 ` profile case, found ----- + + +## Calculate profile volume average | `calculate_vol_avg_of_profile()` + +General method to calculate the volume averaged value of any normalised radial profile. + +### Derivation + +The integral $I$ of any quantity $G(\rho)$ over the volume of the plasma where $\rho$ is the normalised minor radius is: + +$$ +I=\int_{}^{} G(\rho) \ dV +$$ + +The volume element is the surface area of the toroidal surface times the thickness $dV=2\pi R_0 2\pi\kappa r dr $ + +$$ +I=\int_{0}^{1} G(\rho) 2\pi R_0 2\pi\kappa r \ dr +$$ + +Substituting the normalised minor radius coordinate $\rho=r/a$, and bringing the constants out in front, + +$$ +I=4 \pi^2 R_0 \kappa a^2 \int_{0}^{1} G(\rho) \rho \ d\rho +$$ + +However, if we calculate the volume $V$, the exact same integral is used where $G=1$, + +$$ +V=4 \pi^2 R_0 \kappa a^2 \int_{0}^{1} \rho \ d\rho = \frac{1}{2} 4 \pi^2 R_0 \kappa a^2 +$$ + +Therefore the general volume integral of $G(\rho)$, where $V$ is the volume, is + +$$ +I=2 V \int_{0}^{1} G(\rho) \rho \ d\rho +$$ + +If $G$ is a density then $2 \int_{0}^{1} G(\rho) \rho \ d\rho$ is the volume averaged density. + + ## Key Constraints -------- diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 30b2009ee6..e5f342dd46 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -10,6 +10,7 @@ from process.core.data_structure.base import DataStructure from process.data_structure.physics_variables import PhysicsData from process.models.physics.plasma_profiles import PlasmaProfile +from process.models.physics.profiles import calculate_vol_avg_of_profile logger = logging.getLogger(__name__) @@ -80,27 +81,37 @@ class FusionReactionRate: Attributes ---------- - plasma_profile (PlasmaProfile): The parameterized temperature and density - profiles of the plasma. - sigmav_dt_average (float): Average fusion reaction rate for D-T. - dhe3_power_density (float): Fusion power density produced by the D-3He reaction. - dd_power_density (float): Fusion power density produced by the D-D reactions. - dt_power_density (float): Fusion power density produced by the D-T reaction. - alpha_power_density (float): Power density of alpha particles produced. - pden_non_alpha_charged_mw (float): Power density of charged particles produced. - neutron_power_density (float): Power density of neutrons produced. - fusion_rate_density (float): Fusion reaction rate density. - alpha_rate_density (float): Alpha particle production rate density. - proton_rate_density (float): Proton production rate density. - f_dd_branching_trit (float): The rate of tritium producing D-D reactions to - 3He ones. + plasma_profile : PlasmaProfile + Parameterized temperature and density profiles of the plasma. + sigmav_dt_average : float + Volume-averaged D-T fusion reactivity, 〈σv〉ᵥ, for D-T. + dhe3_power_density : float + Fusion power density produced by the D-3He reaction. + dd_power_density : float + Fusion power density produced by the D-D reactions. + dt_power_density : float + Fusion power density produced by the D-T reaction. + alpha_power_density : float + Power density of alpha particles produced. + pden_non_alpha_charged_mw : float + Power density of charged particles produced. + neutron_power_density : float + Power density of neutrons produced. + fusion_rate_density : float + Fusion reaction rate density. + alpha_rate_density : float + Alpha particle production rate density. + proton_rate_density : float + Proton production rate density. + f_dd_branching_trit : float + Ratio of tritium-producing D-D reactions to 3He-producing D-D reactions. References ---------- - - H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections - and thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, - Apr. 1992, doi: https://doi.org/10.1088/0029-5515/32/4/i07. - """ + [1] H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections + and thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, + Apr. 1992, doi: https://doi.org/10.1088/0029-5515/32/4/i07. + """ # noqa: RUF002 def __init__(self, plasma_profile: PlasmaProfile, data: DataStructure): """ @@ -139,15 +150,15 @@ def deuterium_branching(self, ion_temperature: float) -> float: Notes ----- - For ion temperatures between 0.5 keV and 200 keV. - The deviation of the fit from the R-matrix branching ratio is always smaller + - For ion temperatures between 0.5 keV and 200 keV. + - The deviation of the fit from the R-matrix branching ratio is always smaller than 0.5%. References ---------- - - H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections - and thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, - Apr. 1992, doi: https://doi.org/10.1088/0029-5515/32/4/i07. + [1] H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections + and thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, + Apr. 1992, doi: https://doi.org/10.1088/0029-5515/32/4/i07. """ # Divide by 2 to get the branching ratio for the D-D reaction that produces # tritium as the output is just the ratio of the two normalized cross sections @@ -164,34 +175,35 @@ def dt_reaction(self): This method calculates the fusion reaction rate and power density for the deuterium-tritium (D-T) fusion reaction. It uses the Bosch-Hale parametrization - to compute the volumetric fusion reaction rate and integrates over + to compute the volumetric fusion reaction rate 〈σv〉 and integrates over the plasma cross-section to find the core plasma fusion power. The method updates the following attributes: - - self.sigmav_dt_average: Average fusion reaction rate for D-T. + - self.sigmav_dt_average: Volume averaged D-T fusion reactivity 〈σv〉ᵥ for D-T. - self.dt_power_density: Fusion power density produced by the D-T reaction. - self.alpha_power_density: Power density of alpha particles produced. - self.pden_non_alpha_charged_mw: Power density of charged particles - produced. + produced. - self.neutron_power_density: Power density of neutrons produced. - self.fusion_rate_density: Fusion reaction rate density. - self.alpha_rate_density: Alpha particle production rate density. - self.proton_rate_density: Proton production rate density. - - """ + """ # noqa: RUF002 # Initialize Bosch-Hale constants for the D-T reaction dt = BoschHaleConstants(**REACTION_CONSTANTS_DT) - self.data.physics.fusrat_plasma_dt_profile = ( - bosch_hale_reactivity( - ( - self.data.physics.temp_plasma_ion_vol_avg_kev - / self.data.physics.temp_plasma_electron_vol_avg_kev - ) - * self.plasma_profile.teprofile.profile_y, - dt, - ) + ion_temperature_profile = ( + self.data.physics.temp_plasma_ion_vol_avg_kev + / self.data.physics.temp_plasma_electron_vol_avg_kev + ) * self.plasma_profile.teprofile.profile_y + + sigv_profile = bosch_hale_reactivity( + ion_temperature_profile=ion_temperature_profile, reaction_constants=dt + ) + + reaction_rate_profile = ( + sigv_profile * self.data.physics.f_plasma_fuel_deuterium * self.data.physics.f_plasma_fuel_tritium * ( @@ -204,36 +216,19 @@ def dt_reaction(self): ** 2 ) - # Calculate the fusion reaction rate integral using Simpson's rule - sigmav = integrate.simpson( - fusion_rate_integral( - self.plasma_profile, dt, physics_data=self.data.physics - ), - x=self.plasma_profile.neprofile.profile_x, - dx=self.plasma_profile.neprofile.profile_dx, - ) + self.data.physics.fusrat_plasma_dt_profile = reaction_rate_profile - # Store the average fusion reaction rate - self.sigmav_dt_average = sigmav + reaction_rate_density = calculate_vol_avg_of_profile( + profile_x=self.plasma_profile.teprofile.profile_x, + profile_y=reaction_rate_profile, + ) # Reaction energy in MegaJoules [MJ] reaction_energy = constants.D_T_ENERGY / 1.0e6 - # Calculate the fusion power density produced [MW/m^3] - fusion_power_density = ( - sigmav - * reaction_energy - * ( - self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) - * ( - self.data.physics.f_plasma_fuel_tritium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) - ) + fusion_power_density = reaction_rate_density * reaction_energy - # Power densities for different particles [MW/m^3] + # Power densities for different particles [MW/m³] # Alpha particle gets approximately 20% of the fusion power alpha_power_density = ( 1.0 - constants.DT_NEUTRON_ENERGY_FRACTION @@ -243,7 +238,7 @@ def dt_reaction(self): constants.DT_NEUTRON_ENERGY_FRACTION * fusion_power_density ) - # Calculate the fusion rate density [reactions/m^3/second] + # Calculate the fusion rate density [reactions/m³/second] fusion_rate_density = fusion_power_density / reaction_energy alpha_rate_density = fusion_rate_density proton_rate_density = 0.0 @@ -266,7 +261,7 @@ def dhe3_reaction(self): This method calculates the fusion reaction rate and power density for the deuterium-helium-3 (D-3He) fusion reaction. It uses the Bosch-Hale - parametrization to compute the volumetric fusion reaction rate and + parametrization to compute the volumetric fusion reaction rate 〈σv〉 and integrates over the plasma cross-section to find the core plasma fusion power. The method updates the following attributes: @@ -281,28 +276,22 @@ def dhe3_reaction(self): - self.proton_rate_density: Proton production rate density. - """ + """ # noqa: RUF002 # Initialize Bosch-Hale constants for the D-3He reaction dhe3 = BoschHaleConstants(**REACTION_CONSTANTS_DHE3) - # Calculate the fusion reaction rate integral using Simpson's rule - sigmav = integrate.simpson( - fusion_rate_integral( - self.plasma_profile, dhe3, physics_data=self.data.physics - ), - x=self.plasma_profile.neprofile.profile_x, - dx=self.plasma_profile.neprofile.profile_dx, + ion_temperature_profile = ( + self.data.physics.temp_plasma_ion_vol_avg_kev + / self.data.physics.temp_plasma_electron_vol_avg_kev + ) * self.plasma_profile.teprofile.profile_y + + sigv_profile = bosch_hale_reactivity( + ion_temperature_profile=ion_temperature_profile, + reaction_constants=dhe3, ) - self.data.physics.fusrat_plasma_dhe3_profile = ( - bosch_hale_reactivity( - ( - self.data.physics.temp_plasma_ion_vol_avg_kev - / self.data.physics.temp_plasma_electron_vol_avg_kev - ) - * self.plasma_profile.teprofile.profile_y, - dhe3, - ) + reaction_rate_profile = ( + sigv_profile * self.data.physics.f_plasma_fuel_deuterium * self.data.physics.f_plasma_fuel_helium3 * ( @@ -315,24 +304,20 @@ def dhe3_reaction(self): ** 2 ) + self.data.physics.fusrat_plasma_dhe3_profile = reaction_rate_profile + + reaction_rate_density = calculate_vol_avg_of_profile( + profile_x=self.plasma_profile.teprofile.profile_x, + profile_y=reaction_rate_profile, + ) + # Reaction energy in MegaJoules [MJ] reaction_energy = constants.D_HELIUM_ENERGY / 1.0e6 - # Calculate the fusion power density produced [MW/m^3] - fusion_power_density = ( - sigmav - * reaction_energy - * ( - self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) - * ( - self.data.physics.f_plasma_fuel_helium3 - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) - ) + # Calculate the fusion power density produced [MW/m³] + fusion_power_density = reaction_rate_density * reaction_energy - # Power densities for different particles [MW/m^3] + # Power densities for different particles [MW/m³] # Alpha particle gets approximately 20% of the fusion power alpha_power_density = ( 1.0 - constants.DHELIUM_PROTON_ENERGY_FRACTION @@ -342,10 +327,10 @@ def dhe3_reaction(self): ) neutron_power_density = 0.0 - # Calculate the fusion rate density [reactions/m^3/second] + # Calculate the fusion rate density [reactions/m³/second] fusion_rate_density = fusion_power_density / reaction_energy alpha_rate_density = fusion_rate_density - proton_rate_density = fusion_rate_density # Proton production rate [m^3/second] + proton_rate_density = fusion_rate_density # Proton production rate [m³/second] # Update the cumulative D-3He power density self.dhe3_power_density = fusion_power_density @@ -384,26 +369,18 @@ def dd_helion_reaction(self): # Initialize Bosch-Hale constants for the D-D reaction dd1 = BoschHaleConstants(**REACTION_CONSTANTS_DD1) - # Calculate the fusion reaction rate integral using Simpson's rule - sigmav = integrate.simpson( - fusion_rate_integral( - self.plasma_profile, - dd1, - physics_data=self.data.physics, - ), - x=self.plasma_profile.neprofile.profile_x, - dx=self.plasma_profile.neprofile.profile_dx, + ion_temperature_profile = ( + self.data.physics.temp_plasma_ion_vol_avg_kev + / self.data.physics.temp_plasma_electron_vol_avg_kev + ) * self.plasma_profile.teprofile.profile_y + + sigv_profile = bosch_hale_reactivity( + ion_temperature_profile=ion_temperature_profile, + reaction_constants=dd1, ) - self.data.physics.fusrat_plasma_dd_helion_profile = ( - bosch_hale_reactivity( - ( - self.data.physics.temp_plasma_ion_vol_avg_kev - / self.data.physics.temp_plasma_electron_vol_avg_kev - ) - * self.plasma_profile.teprofile.profile_y, - dd1, - ) + reaction_rate_profile = ( + sigv_profile * self.data.physics.f_plasma_fuel_deuterium * self.data.physics.f_plasma_fuel_deuterium * ( @@ -416,6 +393,13 @@ def dd_helion_reaction(self): ** 2 ) + self.data.physics.fusrat_plasma_dd_helion_profile = reaction_rate_profile + + reaction_rate_density = calculate_vol_avg_of_profile( + profile_x=self.plasma_profile.teprofile.profile_x, + profile_y=reaction_rate_profile, + ) + # Reaction energy in MegaJoules [MJ] reaction_energy = constants.DD_HELIUM_ENERGY / 1.0e6 @@ -423,17 +407,7 @@ def dd_helion_reaction(self): # The power density is scaled by the branching ratio to simulate the different # product pathways fusion_power_density = ( - sigmav - * reaction_energy - * (1.0 - self.f_dd_branching_trit) - * ( - self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) - * ( - self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) + reaction_rate_density * reaction_energy * (1.0 - self.f_dd_branching_trit) ) # Power densities for different particles [MW/m^3] @@ -488,26 +462,18 @@ def dd_triton_reaction(self): # Initialize Bosch-Hale constants for the D-D reaction dd2 = BoschHaleConstants(**REACTION_CONSTANTS_DD2) - # Calculate the fusion reaction rate integral using Simpson's rule - sigmav = integrate.simpson( - fusion_rate_integral( - self.plasma_profile, - dd2, - physics_data=self.data.physics, - ), - x=self.plasma_profile.neprofile.profile_x, - dx=self.plasma_profile.neprofile.profile_dx, + ion_temperature_profile = ( + self.data.physics.temp_plasma_ion_vol_avg_kev + / self.data.physics.temp_plasma_electron_vol_avg_kev + ) * self.plasma_profile.teprofile.profile_y + + sigv_profile = bosch_hale_reactivity( + ion_temperature_profile=ion_temperature_profile, + reaction_constants=dd2, ) - self.data.physics.fusrat_plasma_dd_triton_profile = ( - bosch_hale_reactivity( - ( - self.data.physics.temp_plasma_ion_vol_avg_kev - / self.data.physics.temp_plasma_electron_vol_avg_kev - ) - * self.plasma_profile.teprofile.profile_y, - dd2, - ) + reaction_rate_profile = ( + sigv_profile * self.data.physics.f_plasma_fuel_deuterium * self.data.physics.f_plasma_fuel_deuterium * ( @@ -520,6 +486,13 @@ def dd_triton_reaction(self): ** 2 ) + self.data.physics.fusrat_plasma_dd_triton_profile = reaction_rate_profile + + reaction_rate_density = calculate_vol_avg_of_profile( + profile_x=self.plasma_profile.teprofile.profile_x, + profile_y=reaction_rate_profile, + ) + # Reaction energy in MegaJoules [MJ] reaction_energy = constants.DD_TRITON_ENERGY / 1.0e6 @@ -527,17 +500,7 @@ def dd_triton_reaction(self): # The power density is scaled by the branching ratio to simulate the different # product pathways fusion_power_density = ( - sigmav - * reaction_energy - * self.f_dd_branching_trit - * ( - self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) - * ( - self.data.physics.f_plasma_fuel_deuterium - * self.data.physics.nd_plasma_fuel_ions_vol_avg - ) + reaction_rate_density * reaction_energy * self.f_dd_branching_trit ) # Power densities for different particles [MW/m³] @@ -663,12 +626,22 @@ class BoschHaleConstants: cc7: float -def fusion_rate_integral( +def calculate_vol_avg_reactivity( plasma_profile: PlasmaProfile, reaction_constants: BoschHaleConstants, physics_data: PhysicsData, ) -> np.ndarray: - """Evaluate the integrand for the fusion power integration. + """Calculate the volume-averaged fusion reactivity 〈σv〉ᵥ for a given plasma profile + and Bosch-Hale reaction constants. + + This function computes the integrand for the fusion power by normalizing the density + profile with respect to the volume-averaged electron density and multiplying it by + the fusion reactivity. + The result is then integrated over the plasma volume to obtain the volume-averaged + reactivity. + + Parameters + ---------- Parameters ---------- @@ -686,10 +659,10 @@ def fusion_rate_integral( References ---------- - - H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and - thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, Apr. 1992, - doi: https://doi.org/10.1088/0029-5515/32/4/i07. - """ + [1] H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and + thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, Apr. 1992, + doi: https://doi.org/10.1088/0029-5515/32/4/i07. + """ # noqa: RUF002 # Since the electron temperature profile is only calculated directly, we scale the # ion temperature profile by the ratio of the volume averaged ion to electron # temperature @@ -713,15 +686,17 @@ def fusion_rate_integral( # Calculate a volume averaged fusion reaction integral that allows for fusion power # to be scaled with just the volume averaged ion density. - return ( - 2.0 * plasma_profile.teprofile.profile_x * sigv * density_profile_normalised**2 + + return calculate_vol_avg_of_profile( + profile_x=plasma_profile.teprofile.profile_x, + profile_y=sigv * density_profile_normalised**2, ) def bosch_hale_reactivity( ion_temperature_profile: np.ndarray, reaction_constants: BoschHaleConstants ) -> np.ndarray: - """Calculate the volumetric fusion reaction rate 〈sigmav〉 (m³/s) for one of four + """Calculate the velocity-space average reactivity 〈σv〉 [m³/s] for one of four nuclear reactions using the Bosch-Hale parametrization. The valid range of the fit is 0.2 keV < t < 100 keV except for D-3He where it is @@ -743,15 +718,15 @@ def bosch_hale_reactivity( Returns ------- : - np.ndarray: Volumetric fusion reaction rate 〈sigmav〉 in m^3/s for each point in + np.ndarray: Velocity-space average reactivity 〈σv〉 [m³/s] for each point in the ion temperature profile. References ---------- - - H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and - thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, Apr. 1992, - doi: https://doi.org/10.1088/0029-5515/32/4/i07. - """ + [1] H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and + thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611-631, Apr. 1992, + doi: https://doi.org/10.1088/0029-5515/32/4/i07. + """ # noqa: RUF002 theta1 = ( ion_temperature_profile * ( @@ -776,7 +751,7 @@ def bosch_hale_reactivity( xi = ((reaction_constants.bg**2) / (4.0 * theta)) ** (1 / 3) - # Volumetric reaction rate / reactivity 〈sigmav〉 (m³/s) + # Volumetric reaction rate / reactivity 〈σv〉 [m³/s] # noqa: RUF003 # Original form is in [cm³/s], so multiply by 1.0e-6 to convert to [m³/s] sigmav = ( 1.0e-6 diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 53342d9abc..abde62fdec 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -558,3 +558,56 @@ def set_physics_variables(self): / self.data.physics.temp_plasma_electron_vol_avg_kev * self.data.physics.temp_plasma_electron_on_axis_kev ) + + +def calculate_vol_avg_of_profile( + profile_x: np.ndarray, profile_y: np.ndarray, profile_dx: float | None = None +) -> float: + """Calculate the volume averaged value (⟨profile_y⟩) of a radially normalised + profile. + + Parameters + ---------- + profile_x : + The x-values of the profile. + profile_y : + The y-values of the profile. + profile_dx : + The spacing between consecutive x-values in the profile. + + Returns + ------- + float + The volume-averaged value (⟨profile_y⟩) of the profile. + + Raises + ------ + ValueError + If profile_x is not a 1D array, contains fewer than 2 points, + does not span from 0 to 1,or is not strictly increasing. + + + Notes + ----- + - The 2 factor in the calculation arises using both sides of the profile and the + radial normalisation of the profile. + + + """ + if profile_x.ndim != 1: + raise ValueError("profile_x must be a 1D array.") + + if profile_x.size < 2: + raise ValueError("profile_x must contain at least 2 points.") + + if not np.isclose(profile_x[0], 0.0) or not np.isclose(profile_x[-1], 1.0): + raise ValueError("profile_x must span from 0 to 1.") + + if np.any(np.diff(profile_x) <= 0): + raise ValueError("profile_x must be strictly increasing.") + + return 2.0 * sp.integrate.simpson( + profile_y * profile_x, + x=profile_x, + dx=profile_dx if profile_dx is not None else profile_x[1] - profile_x[0], + )