From 0d92ce41128b2e27a58602926701b8324307ab61 Mon Sep 17 00:00:00 2001 From: DarioV86 <85017346+DarioV86@users.noreply.github.com> Date: Wed, 11 Sep 2024 08:50:53 +0100 Subject: [PATCH] Fixing heat flux corrective factor (#3558) * Fixing DN heat flux corrective factor * Added no intersection region to SN * Added notes for the perpendicular transport at the mid-plane in analyse(). Fixed test --- .../advective_transport.py | 171 ++++++++++++++---- .../test_advective_transport.py | 12 +- 2 files changed, 142 insertions(+), 41 deletions(-) diff --git a/bluemira/radiation_transport/advective_transport.py b/bluemira/radiation_transport/advective_transport.py index d99233bee4..503646b94e 100644 --- a/bluemira/radiation_transport/advective_transport.py +++ b/bluemira/radiation_transport/advective_transport.py @@ -19,8 +19,9 @@ from bluemira.base.constants import EPS from bluemira.base.look_and_feel import bluemira_warn from bluemira.display.plotter import Zorder, plot_coordinates -from bluemira.geometry.coordinates import Coordinates +from bluemira.geometry.coordinates import Coordinates, coords_plane_intersect from bluemira.geometry.plane import BluemiraPlane +from bluemira.geometry.tools import make_polygon from bluemira.radiation_transport.error import AdvectionTransportError from bluemira.radiation_transport.flux_surfaces_maker import _clip_flux_surfaces @@ -126,8 +127,7 @@ def _check_params(self): " in either the lower or upper directions." ) - @staticmethod - def _process_first_wall(first_wall): + def _process_first_wall(self, first_wall): """ Force working first wall geometry to be closed and counter-clockwise. """ @@ -142,7 +142,11 @@ def _process_first_wall(first_wall): if not first_wall.closed: bluemira_warn("First wall should be a closed geometry. Closing it.") first_wall.close() - return first_wall + + int_intersection = coords_plane_intersect(first_wall, self._yz_plane)[0] + out_intersection = coords_plane_intersect(first_wall, self._yz_plane)[1] + + return first_wall, int_intersection, out_intersection @staticmethod def _get_arrays(flux_surfaces): @@ -217,6 +221,35 @@ def _clip_flux_surfaces(self, first_wall): ], ) + def _no_wall_intersection_region( + self, x_up_inter, z_up_inter, x_down_inter, z_down_inter, *, lfs=True + ): + """ + Get first wall mid-plane region with no flux line inetrsections. + """ + up_end_i = self.first_wall.argmin(np.array([x_up_inter[-1], 0, z_up_inter[-1]])) + down_end_i = self.first_wall.argmin( + np.array([x_down_inter[-1], 0, z_down_inter[-1]]) + ) + + reg_i = np.nonzero( + (self.first_wall.z < self.first_wall.z[up_end_i]) + & (self.first_wall.z >= self.first_wall.z[down_end_i]) + & ( + (self.first_wall.x > self._o_point.x) + if lfs + else (self.first_wall.x < self._o_point.x) + ) + )[0] + + x_reg_inter = self.first_wall.x[reg_i] + z_reg_inter = self.first_wall.z[reg_i] + + reg_wire = make_polygon(self.first_wall.T[reg_i]) + wire_length = reg_wire.length + + return x_reg_inter, z_reg_inter, wire_length + def analyse(self, first_wall: Coordinates): """ Perform the calculation to obtain charged particle heat fluxes on the @@ -235,8 +268,23 @@ def analyse(self, first_wall: Coordinates): The z coordinates of the flux surface intersections heat_flux: np.array The perpendicular heat fluxes at the intersection points [MW/m^2] + + Notes + ----- + The heat flux model assumes pure parallel transport and fudges + the perpendicular transport via the power decay length, lambda. + This approach, while is widely used, leads to no power deposited + on the wall at the mid-plane. + _analyse_SN and _analyse_DN assume, in the area in proximity of + the mid-plane, where the outermost flux tube is open, that the + remaining power of from the exponential decay is deposited on + the wall perpendicularly. """ - self.first_wall = self._process_first_wall(first_wall) + ( + self.first_wall, + self.imp_int, + self.omp_int, + ) = self._process_first_wall(first_wall) if self.eq.is_double_null: x, z, hf = self._analyse_DN() @@ -249,6 +297,7 @@ def analyse(self, first_wall: Coordinates): def _analyse_SN(self): """ Calculation for the case of single nulls. + """ self._make_flux_surfaces_ob() @@ -284,16 +333,42 @@ def _analyse_SN(self): heat_flux_lfs = self.params.f_lfs_lower_target * q_par_lfs * np.sin(alpha_lfs) heat_flux_hfs = self.params.f_hfs_lower_target * q_par_hfs * np.sin(alpha_hfs) - # Correct power (energy conservation) - q_omp_int = 2 * np.pi * np.sum(q_par_omp / (B_omp / Bp_omp) * self.dx_mp * x_omp) - f_correct_power = self.params.P_sep_particle / q_omp_int + # Find FW portion for perpendicular power + if self.first_wall.argmin( + np.array([x_hfs_inter[-1], 0, z_hfs_inter[-1]]) + ) != self.first_wall.argmin(np.array([x_lfs_inter[-1], 0, z_lfs_inter[-1]])): + x_out_inter, z_out_inter, outb_length = self._no_wall_intersection_region( + x_hfs_inter, z_hfs_inter, x_lfs_inter, z_lfs_inter, lfs=True + ) + else: + mid_i = self.first_wall.argmin( + np.array([x_lfs_inter[-1], 0, z_lfs_inter[-1]]) + ) + x_out_inter = self.first_wall.x[mid_i] + z_out_inter = self.first_wall.z[mid_i] + outb_length = 1 + + # Calculating missing power from parallel transport + q_omp_int = 2 * np.pi * np.sum(q_par_omp * Bp_omp / B_omp * self.dx_mp * x_omp) + miss_omp = self.params.P_sep_particle - q_omp_int + outb_surf = outb_length * 2 * np.pi * self.omp_int[0] + + # Calculating mid-outboard and mid-inboard heat flux + heat_flux_x_outb = miss_omp / outb_surf + if outb_length != 1: + heat_flux_x_outb = [heat_flux_x_outb] * len(x_out_inter) + return ( - np.append(x_lfs_inter, x_hfs_inter), - np.append(z_lfs_inter, z_hfs_inter), - f_correct_power * np.append(heat_flux_lfs, heat_flux_hfs), + np.concatenate([np.atleast_1d(x_out_inter), x_lfs_inter, x_hfs_inter]), + np.concatenate([np.atleast_1d(z_out_inter), z_lfs_inter, z_hfs_inter]), + np.concatenate([ + np.atleast_1d(heat_flux_x_outb), + heat_flux_lfs, + heat_flux_hfs, + ]), ) - def _analyse_DN(self): + def _analyse_DN(self): # noqa: PLR0914 """ Calculation for the case of double nulls. """ @@ -336,7 +411,8 @@ def _analyse_DN(self): Bt_imp = self.eq.Bt(x_imp) B_imp = np.hypot(Bp_imp, Bt_imp) - # Parallel power at the outboard and inboard midplane + # Parallel power set-up at the outboard and inboard midplane + # Note that the power is not split yet into hfs and lfs rates q_par_omp = self._q_par(x_omp, dx_omp, B_omp, Bp_omp) q_par_imp = self._q_par(x_imp, dx_imp, B_imp, Bp_imp, outboard=False) @@ -347,6 +423,7 @@ def _analyse_DN(self): Bp_hfs_up = self.eq.Bp(x_hfs_up_inter, z_hfs_up_inter) # Calculate parallel power at the intersections + # Each q_par_* stores full P_sep_particle # Note that flux expansion terms cancel down to this q_par_lfs_down = q_par_omp * Bp_lfs_down / B_omp q_par_lfs_up = q_par_omp * Bp_lfs_up / B_omp @@ -354,6 +431,7 @@ def _analyse_DN(self): q_par_hfs_up = q_par_imp * Bp_hfs_up / B_imp # Calculate perpendicular heat fluxes + # Here P_sep_particle actually gets distributed over the four targets heat_flux_lfs_down = ( self.params.f_lfs_lower_target * q_par_lfs_down * np.sin(alpha_lfs_down) ) @@ -367,44 +445,65 @@ def _analyse_DN(self): self.params.f_hfs_upper_target * q_par_hfs_up * np.sin(alpha_hfs_up) ) - # Correct power (energy conservation) - q_omp_int = 2 * np.pi * np.sum(q_par_omp * Bp_omp / B_omp * self.dx_mp * x_omp) - q_imp_int = 2 * np.pi * np.sum(q_par_imp * Bp_imp / B_imp * self.dx_mp * x_imp) - - total_power = self.params.P_sep_particle - f_outboard = self.params.f_lfs_lower_target + self.params.f_lfs_upper_target - f_inboard = self.params.f_hfs_lower_target + self.params.f_hfs_upper_target - f_correct_lfs_down = ( - total_power * self.params.f_lfs_lower_target / f_outboard - ) / q_omp_int - f_correct_lfs_up = ( - total_power * self.params.f_lfs_upper_target / f_outboard - ) / q_omp_int - f_correct_hfs_down = ( - total_power * self.params.f_hfs_lower_target / f_inboard - ) / q_imp_int - f_correct_hfs_up = ( - total_power * self.params.f_hfs_upper_target / f_inboard - ) / q_imp_int + # Find FW portion for perpendicular power + x_out_inter, z_out_inter, outb_length = self._no_wall_intersection_region( + x_lfs_up_inter, z_lfs_up_inter, x_lfs_down_inter, z_lfs_down_inter, lfs=True + ) + x_in_inter, z_in_inter, inb_length = self._no_wall_intersection_region( + x_hfs_up_inter, z_hfs_up_inter, x_hfs_down_inter, z_hfs_down_inter, lfs=False + ) + # Calculating missing power from parallel transport + q_omp_int = ( + 2 + * np.pi + * np.sum(q_par_omp * Bp_omp / B_omp * self.dx_mp * x_omp) + * (self.params.f_lfs_lower_target + self.params.f_lfs_upper_target) + ) + q_imp_int = ( + 2 + * np.pi + * np.sum(q_par_imp * Bp_imp / B_imp * self.dx_mp * x_imp) + * (self.params.f_hfs_lower_target + self.params.f_hfs_upper_target) + ) + miss_omp = ( + self.params.P_sep_particle + * (self.params.f_lfs_lower_target + self.params.f_lfs_upper_target) + ) - q_omp_int + miss_imp = ( + self.params.P_sep_particle + * (self.params.f_hfs_lower_target + self.params.f_hfs_upper_target) + ) - q_imp_int + outb_surf = outb_length * 2 * np.pi * self.omp_int[0] + inb_surf = inb_length * 2 * np.pi * self.imp_int[0] + + # Calculating mid-outboard and mid-inboard heat flux + heat_flux_x_outb = [miss_omp / outb_surf] * len(x_out_inter) + heat_flux_x_inb = [miss_imp / inb_surf] * len(x_in_inter) return ( np.concatenate([ + x_out_inter, x_lfs_down_inter, x_lfs_up_inter, + x_in_inter, x_hfs_down_inter, x_hfs_up_inter, ]), np.concatenate([ + z_out_inter, z_lfs_down_inter, z_lfs_up_inter, + z_in_inter, z_hfs_down_inter, z_hfs_up_inter, ]), np.concatenate([ - f_correct_lfs_down * heat_flux_lfs_down, - f_correct_lfs_up * heat_flux_lfs_up, - f_correct_hfs_down * heat_flux_hfs_down, - f_correct_hfs_up * heat_flux_hfs_up, + heat_flux_x_outb, + heat_flux_lfs_down, + heat_flux_lfs_up, + heat_flux_x_inb, + heat_flux_hfs_down, + heat_flux_hfs_up, ]), ) diff --git a/tests/radiation_transport/test_advective_transport.py b/tests/radiation_transport/test_advective_transport.py index 27a94b0969..4aea715fb3 100644 --- a/tests/radiation_transport/test_advective_transport.py +++ b/tests/radiation_transport/test_advective_transport.py @@ -108,20 +108,22 @@ def test_single_null_warnings(self, caplog): def test_recursion(self): assert np.isclose(np.max(self.hf), 5.379, rtol=1e-2) - assert np.argmax(self.hf) == 0 + assert np.argmax(self.hf) == 1 assert np.isclose(np.sum(self.hf), 399, rtol=1e-2) def test_n_intersections(self): """ Because it is a single null, we expect the same number of flux surfaces LFS and HFS. + The final number of intersection points may still be higher than the number of + flux surfaces as there are additional points in the mid-plane region """ len_ob_lfs = len(self.solver.flux_surfaces_ob_down) len_ob_hfs = len(self.solver.flux_surfaces_ob_up) assert len_ob_hfs == len_ob_lfs - assert len_ob_lfs + len_ob_hfs == len(self.x) - assert len(self.solver.flux_surfaces) == len(self.x) + assert len_ob_lfs + len_ob_hfs <= len(self.x) + assert len(self.solver.flux_surfaces) <= len(self.x) def test_integrals(self): n_fs = len(self.solver.flux_surfaces) @@ -235,8 +237,8 @@ def test_double_null_warnings(self, caplog): assert caplog.records[0].levelname == "WARNING" def test_recursion(self): - assert np.isclose(np.max(self.hf), 86.194, rtol=1e-2) - assert np.isclose(np.sum(self.hf), 830.6, rtol=1e-2) + assert np.isclose(np.max(self.hf), 130.992, rtol=1e-2) + assert np.isclose(np.sum(self.hf), 1254.656, rtol=1e-2) def test_analyse_DN(self, caplog): fw = deepcopy(self.solver.first_wall)