RadialDiffusion.cpp

Go to the documentation of this file.

#include "models/nrp/RadialDiffusion.h"
 
#include "models/nrp/PlasmaChannel1D.h"
#include "numeric/DiffusiveFlux.h"
#include "numeric/TableLookup.h"
#include "thermo/PlasmaState.h"
 
#include <cstddef>
 
namespace rizer
{
 
    namespace
    {
 
        // -------------------------------------------------------------------------- //
        //   File-local helpers                                                       //
        // -------------------------------------------------------------------------- //

        template <typename Coefficient>
        double eastFaceFlux(const PlasmaChannel1D &domain,
                            Cantera::span<const double> x, std::size_t component,
                            std::size_t node, double ghost, Coefficient coef)
        {
            const std::size_t n_points = domain.nPoints();
            const double value_here = x[domain.index(component, node)];
 
            const double value_east = (node < n_points - 1)
                                          ? x[domain.index(component, node + 1)]
                                          : ghost;
            const double radius_east = domain.faceRadius(node + 1);
            const double spacing_east = domain.faceSpacing(node + 1);
            // One past the last node names the wall ghost half-cell.
            const double coef_east =
                0.5 * (coef(value_here, node) + coef(value_east, node + 1));
            return numeric::faceFlux(radius_east, coef_east, value_here, value_east,
                                     spacing_east);
        }

        template <typename Coefficient>
        double faceDivergence(const PlasmaChannel1D &domain,
                              Cantera::span<const double> x, std::size_t component,
                              std::size_t node, double ghost, Coefficient coef)
        {
            const double value_here = x[domain.index(component, node)];
            const double flux_east =
                eastFaceFlux(domain, x, component, node, ghost, coef);
            const double radius_east = domain.faceRadius(node + 1);
 
            double flux_west = 0.0;
            double radius_west = 0.0; // [m]
            if (node > 0)
            {
                const double value_west = x[domain.index(component, node - 1)];
                radius_west = domain.faceRadius(node);
                const double coef_west =
                    0.5 * (coef(value_west, node - 1) + coef(value_here, node));
                flux_west =
                    numeric::faceFlux(radius_west, coef_west, value_west,
                                      value_here, domain.faceSpacing(node));
            }
            return (flux_east - flux_west) /
                   numeric::annularVolume(radius_east, radius_west);
        }

        double centredGradient(const PlasmaChannel1D &domain,
                               Cantera::span<const double> x, std::size_t component,
                               std::size_t node, double ghost)
        {
            if (node == 0)
            {
                return 0.0; // axis symmetry: d(phi)/dr(r=0) = 0
            }
            const std::size_t n_points = domain.nPoints();
            const double value_here = x[domain.index(component, node)];
            const double value_west = x[domain.index(component, node - 1)];
            const double gradient_west =
                (value_here - value_west) / domain.faceSpacing(node);
 
            const double value_east = (node < n_points - 1)
                                          ? x[domain.index(component, node + 1)]
                                          : ghost;
            const double gradient_east =
                (value_east - value_here) / domain.faceSpacing(node + 1);
            return 0.5 * (gradient_west + gradient_east);
        }
 
    } // namespace
 
    // ------------------------------------------------------------------------------ //
    //   Free functions                                                               //
    // ------------------------------------------------------------------------------ //
 
    void addRadialDiffusion(const PlasmaChannel1D &domain,
                            Cantera::span<const double> x,
                            const RadialDiffusionCoefficients &c, std::size_t j,
                            double lhsTg, double lhsTe, double &dTg, double &dTe,
                            double *dYdt, double *scratch)
    {
        // Four steps, each adding onto the caller's rates:
        //   conduction in each pool, then the three terms species diffusion
        //   drives (A: the species equation itself, B: the compressibility it
        //   forces through DP/Dt, C: the enthalpy it carries with it).
        const std::vector<double> &rho = domain.density();
        const std::vector<double> &Y_amb = *c.Y_amb;
        const std::vector<double> &molecular_weights = *c.molecular_weights;
 
        // ---- conduction, each pool on its own table -------------------------
        if (!c.kappa->empty())
        {
            dTg += faceDivergence(domain, x, 0, j, c.ambient_temperature,
                                  [&](double T, std::size_t)
                                  { return evalAtState(*c.kappa, T, "radial conduction", "kappa_H(Tg)", "K",
                                                       "Extend the kappa_H_T/kappa_H_v table."); }) /
                   lhsTg;
        }
        if (!c.kappa_e->empty())
        {
            dTe += faceDivergence(domain, x, 1, j, c.ambient_temperature,
                                  [&](double T, std::size_t)
                                  { return evalAtState(*c.kappa_e, T, "radial conduction", "kappa_e(Te)", "K",
                                                       "Extend the kappa_e_T/kappa_e_v table."); }) /
                   lhsTe;
        }
        if (!(c.species_diffusivity > 0.0))
        {
            return;
        }
 
        const double Tg = x[domain.index(0, j)];
        const double Te = x[domain.index(1, j)];
        const double We = molecular_weights[c.electron_index];
        const std::size_t n_points = domain.nPoints();
 
        // Term A (the species equation) and Term B (the compressibility source
        // it drives in both energy equations) both need the divergence of the
        // Fickian flux, div(rho*D*grad(Y_k)). The density sits INSIDE the
        // divergence, where no constant can be pulled out of it unless it
        // happens to be uniform: div(rho*D*grad(Y_k))/rho_j and
        // D*div(grad(Y_k)) are different quantities in general, and only the
        // first conserves species mass -- only its face fluxes telescope
        // against the neighbouring cells'. The face coefficient is averaged
        // from the two adjacent nodes' rho*D, and the wall ghost resolves to
        // the ambient density the constructor closed from p_amb.
        for (std::size_t k = 0; k < c.n_species; k++)
        {
            scratch[k] = faceDivergence(
                domain, x, 2 + k, j, Y_amb[k],
                [&](double, std::size_t n)
                {
                    const double density =
                        (n < n_points) ? rho[n] : c.wall_ghost_density;
                    return density * c.species_diffusivity;
                });
        }
 
        // ---- Term A: Fickian species diffusion ------------------------------
        // Converted to dY_k/dt by dividing through by the node's density.
        for (std::size_t k = 0; k < c.n_species; k++)
        {
            dYdt[k] += scratch[k] / rho[j];
        }
 
        // ---- Term B: species-diffusion compressibility ----------------------
        // From expanding DP/Dt
        // for the two-temperature ideal-gas EOS under Fick's law
        // (PHYSICS.md section 5, the species-diffusion compressibility term).
        // Built on the
        // same flux divergence Term A uses.
        dTe += (GasConstant * Te / We) * scratch[c.electron_index] / lhsTe;
 
        double gas_compressibility_sum = 0.0;
        for (std::size_t k = 0; k < c.n_species; k++)
        {
            if (k == c.electron_index)
                continue;
            gas_compressibility_sum += scratch[k] / molecular_weights[k];
        }
        dTg += GasConstant * Tg * gas_compressibility_sum / lhsTg;
 
        // ---- Term C: enthalpy diffusion -------------------------------------
        // Species heat transport. The full
        // multicomponent energy equation's heat flux includes +sum_k h_k*j_k
        // (species carry their own sensible enthalpy as they diffuse); with
        // Fick's law j_k = -rho*D*grad(Y_k), its divergence -div(sum_k h_k*j_k)
        // reduces to a mixed-gradient PRODUCT rather than a divergence (see
        // PHYSICS.md section 5, the enthalpy-diffusion term):
        //   electron: +rho*D*cp_e * (dYe/dr)(dTe/dr)
        //   gas:      +rho*D*(dTg/dr) * sum_{k!=e} cp_k * (dYk/dr)
        // Being a product of local gradients rather than a divergence, this one
        // takes the node's own density whether or not the density varies with
        // radius.
        const std::vector<double> &species_heat_capacity = *c.species_cp;
        const double Tg_gradient =
            centredGradient(domain, x, 0, j, c.ambient_temperature);
        const double Te_gradient =
            centredGradient(domain, x, 1, j, c.ambient_temperature);
        const double Ye_gradient = centredGradient(
            domain, x, 2 + c.electron_index, j, Y_amb[c.electron_index]);
        dTe += rho[j] * c.species_diffusivity *
               species_heat_capacity[c.electron_index] * Ye_gradient *
               Te_gradient / lhsTe;
 
        double gas_enthalpy_diffusion_sum = 0.0;
        for (std::size_t k = 0; k < c.n_species; k++)
        {
            if (k == c.electron_index)
                continue;
            const double Yk_gradient =
                centredGradient(domain, x, 2 + k, j, Y_amb[k]);
            gas_enthalpy_diffusion_sum +=
                species_heat_capacity[k] * Yk_gradient;
        }
        dTg += rho[j] * c.species_diffusivity * Tg_gradient *
               gas_enthalpy_diffusion_sum / lhsTg;
    }
 
    WallEnergyOutflow wallEnergyOutflow(const PlasmaChannel1D &domain,
                                        Cantera::span<const double> x,
                                        const RadialDiffusionCoefficients &c,
                                        const double *species_enthalpy)
    {
        // Each term is the east-face flux of the outermost node, formed by the
        // same eastFaceFlux() call its divergence makes, negated: that flux is
        // positive INTO the node, and these report what leaves.
        const std::size_t wall_node = domain.nPoints() - 1;
        WallEnergyOutflow outflow;
        if (!c.kappa->empty())
        {
            outflow.conduction -= eastFaceFlux(domain, x, 0, wall_node,
                                               c.ambient_temperature,
                                               [&](double T, std::size_t)
                                               { return evalAtState(*c.kappa, T, "radial conduction", "kappa_H(Tg)", "K",
                                                       "Extend the kappa_H_T/kappa_H_v table."); });
        }
        if (!c.kappa_e->empty())
        {
            outflow.conduction -= eastFaceFlux(domain, x, 1, wall_node,
                                               c.ambient_temperature,
                                               [&](double T, std::size_t)
                                               { return evalAtState(*c.kappa_e, T, "radial conduction", "kappa_e(Te)", "K",
                                                       "Extend the kappa_e_T/kappa_e_v table."); });
        }
        if (!(c.species_diffusivity > 0.0))
        {
            return outflow;
        }
 
        const std::vector<double> &rho = domain.density();
        const std::vector<double> &Y_amb = *c.Y_amb;
        const std::size_t n_points = domain.nPoints();
        outflow.species_mass.assign(c.n_species, 0.0);
        for (std::size_t k = 0; k < c.n_species; k++)
        {
            const double mass_inflow = eastFaceFlux(
                domain, x, 2 + k, wall_node, Y_amb[k],
                [&](double, std::size_t n)
                {
                    const double density =
                        (n < n_points) ? rho[n] : c.wall_ghost_density;
                    return density * c.species_diffusivity;
                });
            outflow.diffusion -= species_enthalpy[k] * mass_inflow;
            outflow.species_mass[k] = -mass_inflow;
        }
        return outflow;
    }
 
} // namespace rizer