diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index 0cb370a7..64e4a29b 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -1150,12 +1150,19 @@ Recycling has been implemented at the target, the SOL edge and the PFR edge. Each is off by default and must be activated with a separate flag. Each can be assigned a separate recycle multiplier and recycle energy. -The chosen species must feature an outflow through the boundary - any cells +Configuring thermal recycling +^^^^^^^^^^^^^^^ + +A simple and commonly used way to model recycling is to assume it is fully thermal, +i.e. that every incident ion recombines into a neutral molecule and thermalises with the surface +before becoming re-emitted. Hermes-3 does not yet have a hydrogenic molecule model, and so +the molecules are assumed to instantly dissociate at the Franck-Condon dissociation temperature of 3.5eV. + +In order to set this up, the chosen species must feature an outflow through the boundary - any cells with an inflow have their recycling source set to zero. If a sheath boundary condition is enabled, then this is automatically satisfied at the target through the Bohm condition. If it is not enabled, then the target boundary must be set to `free_o2`, `free_o3` or `decaylength` to -allow an outflow. If SOL or PFR recycling is enabled, a `free_o2`, `free_o3` or `decaylength`on -their respective boundaries is required at all times. +allow an outflow. The recycling component has a `species` option, that is a list of species to recycle. For each of the species in that list, `recycling` will look in @@ -1193,6 +1200,40 @@ Each returning atom has an energy of 3.5eV: pfr_recycle_multiplier = 1 # Recycling fraction pfr_recycle_energy = 3.5 # Energy of recycled particles [eV] +Allowing for fast recycling +^^^^^^^^^^^^^^^ + +In reality, a fraction of incident ions will undergo specular reflection off the surface and +preserve a fraction of their energy. In the popular Monte-Carlo neutral code EIRENE, the +fast recycling fraction and the energy reflection factor are provided by the `TRIM database `_ +as a function of incident angle, surface material and incident particle energy. +Studies found that sheath acceleration can make the ion angle relatively consistent, e.g. 60 degrees; in (`Jae-Sun Park et al 2021 Nucl. Fusion 61 016021 `_). + +The recycled heat flux is: + +.. math:: + + \begin{aligned} + \Gamma_{E_{n}} &= R \times (R_{f} \alpha_{E} \Gamma_{E_{i}}^{sheath} + (1 - R_{f}) T_{R} \Gamma_{N_{i}})) \\ + \end{aligned} + +Where :math:`R` is the recycle multiplier, :math:`R_{f}` is the fast reflection fraction, :math:`\alpha_{E}` is the energy reflection factor, +:math:`\Gamma_{E_{i}}^{sheath}` is the incident heat flux from the sheath boundary condition, :math:`T_{R}` is the recycle energy and :math:`\Gamma_{N_{i}}` is the incident ion flux. + +:math:`R_{f}` and :math:`\alpha_{E}` can be set as in the below example. They can also be set to different values for the SOL and PFR by replacing +the word "target" with either "sol" or "pfr". + +.. code-block:: ini + + [d+] + recycle_as = d # Species to recycle as + + target_recycle = true + target_recycle_multiplier = 0.95 # Recycling fraction + target_recycle_energy = 3.5 # Energy of recycled particles [eV] + target_fast_recycle_energy_factor = 0.70 + target_fast_recycle_fraction = 0.80 + Neutral pump ^^^^^^^^^^^^^^^ @@ -1204,8 +1245,8 @@ The pump requires wall recycling to be enabled on the relevant wall region. The particle loss rate :math:`\Gamma_{N_{n}}` is the sum of the incident ions that are not recycled and the incident neutrals which are not reflected, both of which are controlled by the pump multiplier :math:`M_{p}` -which is set by the `pump_multiplier` option in the input file. The unrecycled ion flux :math:`\Gamma_{N_{i}}^{unrecycled}` is calculated exactly the same -as for edge recycling but with `pump_multiplier` replacing the `recycle_multiplier`. +which is set by the `pump_multiplier` option in the input file. The unrecycled ion flux :math:`\Gamma_{N_{i}}^{unrecycled}` is calculated using the recycling +model and allows for either thermal or fast recycling, but with the difference that the `pump_multiplier` replaces the `recycle_multiplier`. .. math:: diff --git a/examples/1D-hydrogen-equalT/BOUT.inp b/examples/1D-hydrogen-equalT/BOUT.inp index b164d07f..37eaf69c 100644 --- a/examples/1D-hydrogen-equalT/BOUT.inp +++ b/examples/1D-hydrogen-equalT/BOUT.inp @@ -95,6 +95,8 @@ density_controller_p = 5e2 diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] diff --git a/examples/1D-hydrogen/BOUT.inp b/examples/1D-hydrogen/BOUT.inp index 889602b3..51b0d3d6 100644 --- a/examples/1D-hydrogen/BOUT.inp +++ b/examples/1D-hydrogen/BOUT.inp @@ -98,6 +98,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] diff --git a/examples/1D-hydrogen/qn/BOUT.inp b/examples/1D-hydrogen/qn/BOUT.inp index e1b020d1..37d8ca25 100644 --- a/examples/1D-hydrogen/qn/BOUT.inp +++ b/examples/1D-hydrogen/qn/BOUT.inp @@ -100,6 +100,7 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] diff --git a/examples/1D-neon-source/BOUT.inp b/examples/1D-neon-source/BOUT.inp index ba0b3cc5..0204a8fc 100644 --- a/examples/1D-neon-source/BOUT.inp +++ b/examples/1D-neon-source/BOUT.inp @@ -103,6 +103,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] @@ -169,6 +171,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+] @@ -189,6 +193,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+2] @@ -209,6 +215,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+3] @@ -229,6 +237,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -251,8 +261,9 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction - diagnose = true [Nne+5] @@ -273,6 +284,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -295,6 +308,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -317,6 +332,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -339,6 +356,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -361,6 +380,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true diff --git a/examples/1D-neon/BOUT.inp b/examples/1D-neon/BOUT.inp index 3f4b4d19..2e0b9376 100644 --- a/examples/1D-neon/BOUT.inp +++ b/examples/1D-neon/BOUT.inp @@ -104,7 +104,9 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d -target_recycle_multiplier = 1.0 # Recycling fraction +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy +target_recycle_multiplier = 1 # Recycling fraction [Nd+] @@ -168,6 +170,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction [Nne+] @@ -188,6 +192,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction [Nne+2] @@ -208,6 +214,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction [Nne+3] @@ -228,6 +236,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -250,6 +260,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -272,6 +284,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -294,6 +308,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -316,6 +332,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -338,6 +356,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -360,6 +380,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true diff --git a/examples/1D-recycling/BOUT.inp b/examples/1D-recycling/BOUT.inp index d3ca2049..73666e46 100644 --- a/examples/1D-recycling/BOUT.inp +++ b/examples/1D-recycling/BOUT.inp @@ -101,6 +101,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/1D-recycling/cvode/BOUT.inp b/examples/1D-recycling/cvode/BOUT.inp index f0f11284..03713ba0 100644 --- a/examples/1D-recycling/cvode/BOUT.inp +++ b/examples/1D-recycling/cvode/BOUT.inp @@ -92,6 +92,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true recycle_as = d +target_recycle = true +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/solkit-comparison/single-temperature/BOUT.inp b/examples/solkit-comparison/single-temperature/BOUT.inp index 475b0efb..342e3f4e 100644 --- a/examples/solkit-comparison/single-temperature/BOUT.inp +++ b/examples/solkit-comparison/single-temperature/BOUT.inp @@ -100,7 +100,9 @@ AA = 2 diagnose = true # Output diagnostics for these components? recycle_as = d -target_recycle_multiplier = 1.0 # Recycling fraction +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy +target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/solkit-comparison/two-temperatures/BOUT.inp b/examples/solkit-comparison/two-temperatures/BOUT.inp index 6f429253..b7c91ee8 100644 --- a/examples/solkit-comparison/two-temperatures/BOUT.inp +++ b/examples/solkit-comparison/two-temperatures/BOUT.inp @@ -100,7 +100,9 @@ AA = 2 diagnose = true # Output diagnostics for these components? recycle_as = d -target_recycle_multiplier = 1.0 # Recycling fraction +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy +target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/tokamak/recycling-dthe-drifts/BOUT.inp b/examples/tokamak/recycling-dthe-drifts/BOUT.inp index 01cbe8bc..292892ce 100644 --- a/examples/tokamak/recycling-dthe-drifts/BOUT.inp +++ b/examples/tokamak/recycling-dthe-drifts/BOUT.inp @@ -53,6 +53,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nd+] @@ -89,6 +90,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = t +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nt+] @@ -125,6 +127,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = he +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nhe+] diff --git a/examples/tokamak/recycling-dthe/BOUT.inp b/examples/tokamak/recycling-dthe/BOUT.inp index ced75dd8..74fa2eed 100644 --- a/examples/tokamak/recycling-dthe/BOUT.inp +++ b/examples/tokamak/recycling-dthe/BOUT.inp @@ -53,6 +53,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nd+] @@ -89,6 +90,7 @@ anomalous_chi = 1.0 thermal_conduction = true recycle_as = t +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nt+] @@ -125,6 +127,7 @@ anomalous_chi = 1.0 thermal_conduction = true recycle_as = he +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true diff --git a/examples/tokamak/recycling-dthene/BOUT.inp b/examples/tokamak/recycling-dthene/BOUT.inp index 61d87c1f..21c65bb0 100644 --- a/examples/tokamak/recycling-dthene/BOUT.inp +++ b/examples/tokamak/recycling-dthene/BOUT.inp @@ -36,6 +36,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nd+] @@ -76,6 +77,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = t +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nt+] @@ -116,6 +118,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = he +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nhe+] @@ -156,6 +159,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = ne +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nne+] diff --git a/examples/tokamak/recycling/BOUT.inp b/examples/tokamak/recycling/BOUT.inp index 57a9b35d..310d0ee0 100644 --- a/examples/tokamak/recycling/BOUT.inp +++ b/examples/tokamak/recycling/BOUT.inp @@ -59,6 +59,7 @@ anomalous_chi = 4 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true diff --git a/include/recycling.hxx b/include/recycling.hxx index 696ba3d6..69c355fb 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -48,6 +48,8 @@ private: /// Combination of recycling fraction and species change e.g h+ -> h2 results in 0.5 multiplier BoutReal target_multiplier, sol_multiplier, pfr_multiplier, pump_multiplier; BoutReal target_energy, sol_energy, pfr_energy; ///< Energy of recycled particle (normalised to Tnorm) + BoutReal target_fast_recycle_fraction, pfr_fast_recycle_fraction, sol_fast_recycle_fraction; ///< Fraction of ions undergoing fast reflection + BoutReal target_fast_recycle_energy_factor, sol_fast_recycle_energy_factor, pfr_fast_recycle_energy_factor; ///< Fraction of energy retained by fast recycled neutrals }; std::vector channels; // Recycling channels @@ -56,16 +58,17 @@ private: bool diagnose; ///< Save additional post-processing variables? Field3D density_source, energy_source; ///< Recycling particle and energy sources for all locations - Field2D is_pump; ///< 1 = pump, 0 = no pump. Works only in SOL/PFR + Field3D energy_flow_ylow, energy_flow_xlow; ///< Cell edge fluxes used for calculating fast recycling energy source + Field3D particle_flow_xlow; ///< Radial wall particle fluxes for recycling calc. No need to get poloidal from here, it's calculated from sheath velocity + + Field2D is_pump; ///< 1 = pump, 0 = no pump. Works only in SOL/PFR. Provided by user in grid file. // Recycling particle and energy sources for the different sources of recycling // Note that SOL, PFR and pump are not applicable to 1D Field3D target_recycle_density_source, target_recycle_energy_source; Field3D wall_recycle_density_source, wall_recycle_energy_source; ///< Recycling particle and energy sources for pfr + sol recycling - Field3D pump_recycle_density_source, pump_recycle_energy_source; + Field3D pump_density_source, pump_energy_source; ///< Recycling particle and energy sources for pump recycling - Field3D radial_particle_outflow, radial_energy_outflow; ///< Radial fluxes coming from evolve_density and evolve_pressure used in recycling calc - }; namespace { diff --git a/src/recycling.cxx b/src/recycling.cxx index e81bf377..68cb7e48 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -3,6 +3,7 @@ #include // for trim, strsplit #include "../include/hermes_utils.hxx" // For indexAt +#include "../include/hermes_utils.hxx" // For indexAt #include #include #include @@ -22,13 +23,13 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { .as(), ','); - + // Neutral pump // Mark cells as having a pump by setting the Field2D is_pump to 1 in the grid file // Works only on SOL and PFR edges, where it locally modifies the recycle multiplier to the pump albedo is_pump = 0.0; mesh->get(is_pump, std::string("is_pump")); - + for (const auto& species : species_list) { std::string from = trim(species, " \t\r()"); // The species name in the list @@ -63,7 +64,6 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { from_options["pump_recycle_multiplier"] .doc("Multiply the pump boundary recycling flux by this factor (like albedo). Should be >=0 and <= 1") .withDefault(1.0); - BoutReal target_recycle_energy = from_options["target_recycle_energy"] .doc("Fixed energy of the recycled particles at target [eV]") @@ -80,6 +80,36 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { .withDefault(3.0) / Tnorm; // Normalise from eV + BoutReal target_fast_recycle_fraction = + from_options["target_fast_recycle_fraction"] + .doc("Fraction of ions undergoing fast reflection at target") + .withDefault(0); + + BoutReal pfr_fast_recycle_fraction = + from_options["pfr_fast_recycle_fraction"] + .doc("Fraction of ions undergoing fast reflection at pfr") + .withDefault(0); + + BoutReal sol_fast_recycle_fraction = + from_options["sol_fast_recycle_fraction"] + .doc("Fraction of ions undergoing fast reflection at sol") + .withDefault(0); + + BoutReal target_fast_recycle_energy_factor = + from_options["target_fast_recycle_energy_factor"] + .doc("Fraction of energy retained by fast recycled neutrals at target") + .withDefault(0); + + BoutReal sol_fast_recycle_energy_factor = + from_options["sol_fast_recycle_energy_factor"] + .doc("Fraction of energy retained by fast recycled neutrals at sol") + .withDefault(0); + + BoutReal pfr_fast_recycle_energy_factor = + from_options["pfr_fast_recycle_energy_factor"] + .doc("Fraction of energy retained by fast recycled neutrals at pfr") + .withDefault(0); + if ((target_recycle_multiplier < 0.0) or (target_recycle_multiplier > 1.0) or (sol_recycle_multiplier < 0.0) or (sol_recycle_multiplier > 1.0) or (pfr_recycle_multiplier < 0.0) or (pfr_recycle_multiplier > 1.0) @@ -91,7 +121,9 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { channels.push_back({ from, to, target_recycle_multiplier, sol_recycle_multiplier, pfr_recycle_multiplier, pump_recycle_multiplier, - target_recycle_energy, sol_recycle_energy, pfr_recycle_energy}); + target_recycle_energy, sol_recycle_energy, pfr_recycle_energy, + target_fast_recycle_fraction, pfr_fast_recycle_fraction, sol_fast_recycle_fraction, + target_fast_recycle_energy_factor, sol_fast_recycle_energy_factor, pfr_fast_recycle_energy_factor}); // Boolean flags for enabling recycling in different regions target_recycle = from_options["target_recycle"] @@ -108,7 +140,7 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { neutral_pump = from_options["neutral_pump"] .doc("Neutral pump enabled? Note, need location in grid file") - .withDefault(false); + .withDefault(false); } } @@ -129,9 +161,9 @@ void Recycling::transform(Options& state) { const Field3D N = get(species_from["density"]); const Field3D V = get(species_from["velocity"]); // Parallel flow velocity + const Field3D T = get(species_from["temperature"]); // Ion temperature Options& species_to = state["species"][channel.to]; - const Field3D Nn = get(species_to["density"]); const Field3D Pn = get(species_to["pressure"]); const Field3D Tn = get(species_to["temperature"]); @@ -146,9 +178,20 @@ void Recycling::transform(Options& state) { ? getNonFinal(species_to["energy_source"]) : 0.0; + // Recycling at the divertor target plates if (target_recycle) { + // Fast recycling needs to know how much energy the "from" species is losing to the boundary + if (species_from.isSet("energy_flow_ylow")) { + energy_flow_ylow = get(species_from["energy_flow_ylow"]); + } else { + energy_flow_ylow = 0; + if (channel.target_fast_recycle_fraction > 0) { + throw BoutException("Target fast recycle enabled but no sheath heat flow available in energy_flow_ylow! \nCurrently only sheath_boundary_simple is supported for fast recycling."); + } + } + target_recycle_density_source = 0; target_recycle_energy_source = 0; @@ -166,23 +209,30 @@ void Recycling::transform(Options& state) { flux = 0.0; } - // Flow of recycled species inwards + // Flow of recycled neutrals into domain [s-1] BoutReal flow = channel.target_multiplier * flux - * (J(r.ind, mesh->ystart) + J(r.ind, mesh->ystart - 1)) - / (sqrt(g_22(r.ind, mesh->ystart)) + sqrt(g_22(r.ind, mesh->ystart - 1))); - - // Add to density source - target_recycle_density_source(r.ind, mesh->ystart, jz) += flow - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - density_source(r.ind, mesh->ystart, jz) += flow - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - - // energy of recycled particles - target_recycle_energy_source(r.ind, mesh->ystart, jz) += channel.target_energy * flow - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - energy_source(r.ind, mesh->ystart, jz) += channel.target_energy * flow - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); + * (J(r.ind, mesh->ystart) + J(r.ind, mesh->ystart - 1)) / (sqrt(g_22(r.ind, mesh->ystart)) + sqrt(g_22(r.ind, mesh->ystart - 1))) + * 0.5*(dx(r.ind, mesh->ystart) + dx(r.ind, mesh->ystart - 1)) * 0.5*(dz(r.ind, mesh->ystart) + dz(r.ind, mesh->ystart - 1)); + + BoutReal volume = J(r.ind, mesh->ystart) * dx(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart) * dz(r.ind, mesh->ystart); + + // Calculate sources in the final cell [m^-3 s^-1] + target_recycle_density_source(r.ind, mesh->ystart, jz) += flow / volume; // For diagnostic + density_source(r.ind, mesh->ystart, jz) += flow / volume; // For use in solver + + // Energy of recycled particles + BoutReal ion_energy_flow = energy_flow_ylow(r.ind, mesh->ystart, jz) * -1; // This is ylow end so take first domain cell and flip sign + + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal recycle_energy_flow = + ion_energy_flow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part + + // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] + target_recycle_energy_source(r.ind, mesh->ystart, jz) += recycle_energy_flow / volume; + energy_source(r.ind, mesh->ystart, jz) += recycle_energy_flow / volume; } } @@ -201,37 +251,64 @@ void Recycling::transform(Options& state) { flux = 0.0; } - // Flow of neutrals inwards + // Flow of recycled neutrals into domain [s-1] BoutReal flow = - channel.target_multiplier * flux * (J(r.ind, mesh->yend) + J(r.ind, mesh->yend + 1)) - / (sqrt(g_22(r.ind, mesh->yend)) + sqrt(g_22(r.ind, mesh->yend + 1))); - - // Rate of change of neutrals in final cell - // Add to density source - target_recycle_density_source(r.ind, mesh->yend, jz) += flow - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - density_source(r.ind, mesh->yend, jz) += flow - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - - target_recycle_energy_source(r.ind, mesh->yend, jz) += channel.target_energy * flow - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - energy_source(r.ind, mesh->yend, jz) += channel.target_energy * flow - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); + channel.target_multiplier * flux + * (J(r.ind, mesh->yend) + J(r.ind, mesh->yend + 1)) / (sqrt(g_22(r.ind, mesh->yend)) + sqrt(g_22(r.ind, mesh->yend + 1))) + * 0.5*(dx(r.ind, mesh->yend) + dx(r.ind, mesh->yend + 1)) * 0.5*(dz(r.ind, mesh->yend) + dz(r.ind, mesh->yend + 1)); + + BoutReal volume = J(r.ind, mesh->yend) * dx(r.ind, mesh->yend) * dy(r.ind, mesh->yend) * dz(r.ind, mesh->yend); + + // Calculate sources in the final cell [m^-3 s^-1] + target_recycle_density_source(r.ind, mesh->yend, jz) += flow / volume; // For diagnostic + density_source(r.ind, mesh->yend, jz) += flow / volume; // For use in solver + + // Energy of recycled particles + BoutReal ion_energy_flow = energy_flow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell + + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal recycle_energy_flow = + ion_energy_flow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part + + + // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] + target_recycle_energy_source(r.ind, mesh->yend, jz) += recycle_energy_flow / volume; + energy_source(r.ind, mesh->yend, jz) += recycle_energy_flow / volume; } } } // Initialise counters of pump recycling fluxes - pump_recycle_density_source = 0; - pump_recycle_energy_source = 0; wall_recycle_density_source = 0; wall_recycle_energy_source = 0; + pump_density_source = 0; + pump_energy_source = 0; + + // Get edge particle and heat for the species being recycled + if (sol_recycle or pfr_recycle) { + + if (species_from.isSet("energy_flow_xlow")) { + energy_flow_xlow = get(species_from["energy_flow_xlow"]); + } else if ((channel.sol_fast_recycle_fraction > 0) or (channel.pfr_fast_recycle_fraction > 0)) { + throw BoutException("SOL/PFR fast recycle enabled but no cell edge heat flow available, check your wall BC choice"); + }; + + if (species_from.isSet("particle_flow_xlow")) { + particle_flow_xlow = get(species_from["particle_flow_xlow"]); + } else if ((channel.sol_fast_recycle_fraction > 0) or (channel.pfr_fast_recycle_fraction > 0)) { + throw BoutException("SOL/PFR fast recycle enabled but no cell edge particle flow available, check your wall BC choice"); + }; + + } // Recycling at the SOL edge (2D/3D only) if (sol_recycle) { // Flow out of domain is positive in the positive coordinate direction - radial_particle_outflow = get(species_from["particle_flow_xlow"]); + Field3D radial_particle_outflow = particle_flow_xlow; + Field3D radial_energy_outflow = energy_flow_xlow; if(mesh->lastX()){ // Only do this for the processor which has the edge region for(int iy=0; iy < mesh->LocalNy ; iy++){ @@ -242,12 +319,12 @@ void Recycling::transform(Options& state) { * dy(mesh->xend, iy) * dz(mesh->xend, iy); // If cell is a pump, overwrite multiplier with pump multiplier - BoutReal multiplier = channel.pfr_multiplier; + BoutReal multiplier = channel.sol_multiplier; if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { multiplier = channel.pump_multiplier; } - // Flow of recycled species back from the edge + // Flow of recycled species back from the edge due to recycling // SOL edge = LHS flow of inner guard cells on the high X side (mesh->xend+1) // Recycling source is 0 for each cell where the flow goes into instead of out of the domain BoutReal recycle_particle_flow = 0; @@ -255,14 +332,19 @@ void Recycling::transform(Options& state) { recycle_particle_flow = multiplier * radial_particle_outflow(mesh->xend+1, iy, iz); } - // Divide by volume to get source - BoutReal recycle_source = recycle_particle_flow / volume; + BoutReal ion_energy_flow = radial_energy_outflow(mesh->xend+1, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take guard cell - // Pump source terms - BoutReal esink = 0; - BoutReal psink = 0; + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal recycle_energy_flow = + ion_energy_flow * channel.sol_multiplier * channel.sol_fast_recycle_energy_factor * channel.sol_fast_recycle_fraction // Fast recycling part + + recycle_particle_flow * (1 - channel.sol_fast_recycle_fraction) * channel.sol_energy; // Thermal recycling part + + // Calculate neutral pump neutral sinks due to neutral impingement + // These are additional to particle sinks due to recycling + BoutReal pump_neutral_energy_sink = 0; + BoutReal pump_neutral_particle_sink = 0; - // Compute the neutral loss sink and combine with the recycled source to compute overall neutral sources if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { auto i = indexAt(Nn, mesh->xend, iy, iz); // Final domain cell @@ -293,32 +375,29 @@ void Recycling::transform(Options& state) { // Calculate particle and energy fluxes of neutrals hitting the pump // Assume thermal velocity greater than perpendicular velocity and use it for flux calc - BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] - psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + BoutReal pump_neutral_particle_flow = v_th * nnsheath * dasheath; // [s^-1] + pump_neutral_particle_sink = pump_neutral_particle_flow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] // Use gamma=2.0 as per Stangeby p.69, total energy of static Maxwellian - BoutReal eflow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] - esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] - - // Pump puts neutral particle and energy source in final domain cell - // Source accounts for recycled ions and the particle sink due to neutrals hitting the pump - // Note that the ions are explicitly transported out of the domain by anomalous_diffusion.cxx - // the pump picks this up and adds a recycling source based on this, but doesn't need an ion sink. - // Neutrals are not explicitly transported out by any component and must be taken out by the sink below. - // Pump multiplier controls both fraction of recycled ions and fraction of returned neutrals - pump_recycle_density_source(mesh->xend, iy, iz) += recycle_source - psink; - pump_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy - esink; + BoutReal pump_neutral_energy_flow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] + pump_neutral_energy_sink = pump_neutral_energy_flow / dv * (1 - multiplier); // heatsink [W m^-3] + + // Divide flows by volume to get sources + // Save to pump diagnostic + pump_density_source(mesh->xend, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + pump_energy_source(mesh->xend, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + } else { - wall_recycle_density_source(mesh->xend, iy, iz) += recycle_source; - wall_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy; + // Save to wall diagnostic (pump sinks are 0 if not on pump) + wall_recycle_density_source(mesh->xend, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + wall_recycle_energy_source(mesh->xend, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; } - // Add to density source which will be picked up by evolve_density.cxx - // Add to energy source which will be picked up by evolve_pressure.cxx - // psink and esink are additional sinks from the neutral pump which are 0 if disabled - density_source(mesh->xend, iy, iz) += recycle_source - psink; - energy_source(mesh->xend, iy, iz) += recycle_source * channel.pfr_energy - esink; + // Save to solver + density_source(mesh->xend, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + energy_source(mesh->xend, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + } } @@ -329,25 +408,25 @@ void Recycling::transform(Options& state) { if (pfr_recycle) { // PFR is flipped compared to edge: x=0 is at the PFR edge. Therefore outflow is in the negative coordinate direction. - radial_particle_outflow = get(species_from["particle_flow_xlow"]) * -1; + Field3D radial_particle_outflow = particle_flow_xlow * -1; + Field3D radial_energy_outflow = energy_flow_xlow * -1; if(mesh->firstX()){ // Only do this for the processor which has the core region if (!mesh->periodicY(mesh->xstart)) { // Only do this for the processor with a periodic Y, i.e. the PFR for(int iy=0; iy < mesh->LocalNy ; iy++){ for(int iz=0; iz < mesh->LocalNz; iz++){ - // Volume of cell adjacent to wall which will receive source BoutReal volume = J(mesh->xstart, iy) * dx(mesh->xstart, iy) * dy(mesh->xstart, iy) * dz(mesh->xstart, iy); // If cell is a pump, overwrite multiplier with pump multiplier BoutReal multiplier = channel.pfr_multiplier; - if ((is_pump(mesh->xstart, iy, iz) == 1.0) and (neutral_pump)) { + if ((is_pump(mesh->xstart, iy) == 1.0) and (neutral_pump)) { multiplier = channel.pump_multiplier; } - - // Flow of recycled species back from the edge + + // Flow of recycled species back from the edge due to recycling // PFR edge = LHS flow of the first domain cell on the low X side (mesh->xstart) // Recycling source is 0 for each cell where the flow goes into instead of out of the domain BoutReal recycle_particle_flow = 0; @@ -355,15 +434,22 @@ void Recycling::transform(Options& state) { recycle_particle_flow = multiplier * radial_particle_outflow(mesh->xstart, iy, iz); } - // Divide by volume to get source - BoutReal recycle_source = recycle_particle_flow / volume; + BoutReal ion_energy_flow = radial_energy_outflow(mesh->xstart, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take first domain cell + + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal recycle_energy_flow = + ion_energy_flow * channel.pfr_multiplier * channel.pfr_fast_recycle_energy_factor * channel.pfr_fast_recycle_fraction // Fast recycling part + + recycle_particle_flow * (1 - channel.pfr_fast_recycle_fraction) * channel.pfr_energy; // Thermal recycling part - // Pump source terms - BoutReal esink = 0; - BoutReal psink = 0; + // Calculate neutral pump neutral sinks due to neutral impingement + // These are additional to particle sinks due to recycling + BoutReal pump_neutral_energy_sink = 0; + BoutReal pump_neutral_particle_sink = 0; // Add to appropriate diagnostic field depending if pump or not if ((is_pump(mesh->xstart, iy) == 1.0) and (neutral_pump)) { + auto i = indexAt(Nn, mesh->xstart, iy, iz); // Final domain cell auto is = i.xp(); // Second to final domain cell auto ig = i.xm(); // First guard cell @@ -392,36 +478,35 @@ void Recycling::transform(Options& state) { // Calculate particle and energy fluxes of neutrals hitting the pump // Assume thermal velocity greater than perpendicular velocity and use it for flux calc - BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] - psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + BoutReal pump_neutral_particle_flow = v_th * nnsheath * dasheath; // [s^-1] + pump_neutral_particle_sink = pump_neutral_particle_flow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] // Use gamma=2.0 as per Stangeby p.69, total energy of static Maxwellian - BoutReal eflow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] - esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] - - // Pump puts neutral particle and energy source in final domain cell - // Source accounts for recycled ions and the particle sink due to neutrals hitting the pump - // Note that the ions are explicitly transported out of the domain by anomalous_diffusion.cxx - // the pump picks this up and adds a recycling source based on this, but doesn't need an ion sink. - // Neutrals are not explicitly transported out by any component and must be taken out by the sink below. - // Pump multiplier controls both fraction of recycled ions and fraction of returned neutrals - pump_recycle_density_source(mesh->xend, iy, iz) += recycle_source - psink; - pump_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy - esink; + BoutReal pump_neutral_energy_flow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] + pump_neutral_energy_sink = pump_neutral_energy_flow / dv * (1 - multiplier); // heatsink [W m^-3] + + // Divide flows by volume to get sources + // Save to pump diagnostic + pump_density_source(mesh->xstart, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + pump_energy_source(mesh->xstart, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + + } else { - wall_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; - wall_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; + // Save to wall diagnostic (pump sinks are 0 if not on pump) + wall_recycle_density_source(mesh->xstart, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + wall_recycle_energy_source(mesh->xstart, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + } - - // Add to density source which will be picked up by evolve_density.cxx - // Add to energy source which will be picked up by evolve_pressure.cxx - // psink and esink are additional sinks from the neutral pump which are 0 if disabled - density_source(mesh->xend, iy, iz) += recycle_source - psink; - energy_source(mesh->xend, iy, iz) += recycle_source * channel.pfr_energy - esink; + + // Save to solver + density_source(mesh->xstart, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + energy_source(mesh->xstart, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; } } } } + } // Put the updated sources back into the state @@ -482,11 +567,11 @@ void Recycling::outputVars(Options& state) { {"standard_name", "energy source"}, {"long_name", std::string("Wall recycling energy source of ") + channel.to}, {"source", "recycling"}}); - } + } // Neutral pump if (neutral_pump) { - set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump")}], pump_recycle_density_source, + set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump")}], pump_density_source, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, {"conversion", Nnorm * Omega_ci}, @@ -494,7 +579,7 @@ void Recycling::outputVars(Options& state) { {"long_name", std::string("Pump recycling particle source of ") + channel.to}, {"source", "recycling"}}); - set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump")}], pump_recycle_energy_source, + set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump")}], pump_energy_source, {{"time_dimension", "t"}, {"units", "W m^-3"}, {"conversion", Pnorm * Omega_ci}, diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index 62513b86..f987d787 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -279,7 +279,10 @@ void SheathBoundarySimple::transform(Options& state) { } } } - + + // Field to capture total sheath heat flux for diagnostics + Field3D electron_sheath_power_ylow = zeroFrom(Ne); + ////////////////////////////////////////////////////////////////// // Electrons @@ -327,13 +330,21 @@ void SheathBoundarySimple::transform(Options& state) { * nesheath * vesheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + BoutReal heatflow = q * (coord->J[i] + coord->J[im]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); electron_energy_source[i] += power; + + // Total heat flux for diagnostic purposes + q = gamma_e * tesheath * nesheath * vesheath; // Wm-2 + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W + + electron_sheath_power_ylow[i] += heatflow; // lower Y, so power placed in final domain cell + } } } @@ -379,13 +390,20 @@ void SheathBoundarySimple::transform(Options& state) { * nesheath * vesheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal heatflow = q * (coord->J[i] + coord->J[ip]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); electron_energy_source[i] -= power; + + // Total heat flux for diagnostic purposes + q = gamma_e * tesheath * nesheath * vesheath; // Wm-2 + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W + + electron_sheath_power_ylow[ip] -= heatflow; // upper Y, so power placed in first guard cell } } } @@ -399,6 +417,9 @@ void SheathBoundarySimple::transform(Options& state) { // Note: electron_energy_source includes any sources previously set in other components set(electrons["energy_source"], fromFieldAligned(electron_energy_source)); + // Add the total sheath power flux to the tracker of y power flows + add(electrons["energy_flow_ylow"], fromFieldAligned(electron_sheath_power_ylow)); + if (IS_SET_NOBOUNDARY(electrons["velocity"])) { setBoundary(electrons["velocity"], fromFieldAligned(Ve)); } @@ -452,6 +473,9 @@ void SheathBoundarySimple::transform(Options& state) { ? toFieldAligned(getNonFinal(species["energy_source"])) : zeroFrom(Ni); + // Field to capture total sheath heat flux for diagnostics + Field3D ion_sheath_power_ylow = zeroFrom(Ne); + if (lower_y) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { @@ -496,13 +520,20 @@ void SheathBoundarySimple::transform(Options& state) { * nisheath * visheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + BoutReal heatflow = q * (coord->J[i] + coord->J[im]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); energy_source[i] += power; + + // Calculation of total heat flux for diagnostic purposes + q = gamma_i * tisheath * nisheath * visheath; // Wm-2 + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W + + ion_sheath_power_ylow[i] += heatflow; // lower Y, so power placed in final domain cell } } } @@ -553,18 +584,25 @@ void SheathBoundarySimple::transform(Options& state) { * nisheath * visheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal heatflow = q * (coord->J[i] + coord->J[ip]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); ASSERT2(std::isfinite(power)); energy_source[i] -= power; // Note: Sign negative because power > 0 + + // Calculation of total heat flux for diagnostic purposes + q = gamma_i * tisheath * nisheath * visheath; + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W + + ion_sheath_power_ylow[ip] += heatflow; // Upper Y, so power placed in first guard cell } } + } - // Finished boundary conditions for this species // Put the modified fields back into the state. setBoundary(species["density"], fromFieldAligned(Ni)); @@ -582,5 +620,8 @@ void SheathBoundarySimple::transform(Options& state) { // Additional loss of energy through sheath // Note: energy_source already includes previously set values set(species["energy_source"], fromFieldAligned(energy_source)); + + // Add the total sheath power flux to the tracker of y power flows + add(species["energy_flow_ylow"], fromFieldAligned(ion_sheath_power_ylow)); } }