DiffusiveFlux.h File Reference

Pure-arithmetic primitives for a 1D cylindrical (radial) finite-volume discretization of the diffusive operator \(\frac1r\frac{d}{dr}\!\left(r\,\Gamma\,\frac{d\phi}{dr}\right)\): one face's conductive/diffusive flux, and the annular cell volume between two faces. More...

Go to the source code of this file.

Namespaces
namespace	rizer
namespace	rizer::numeric

Functions
double	rizer::numeric::faceFlux (double r_face, double coef_face, double value_inner, double value_outer, double spacing)
	Conductive/diffusive flux across one radial cell face: r_face * coef_face * (value_outer - value_inner) / spacing.
double	rizer::numeric::annularVolume (double r_east, double r_west)
	Annular cell volume per unit axial length between two face radii, 0.5 * (r_east^2 - r_west^2) [m^2].

Detailed Description

Pure-arithmetic primitives for a 1D cylindrical (radial) finite-volume discretization of the diffusive operator \(\frac1r\frac{d}{dr}\!\left(r\,\Gamma\,\frac{d\phi}{dr}\right)\): one face's conductive/diffusive flux, and the annular cell volume between two faces.

Shared by PlasmaChannel1D (2T reacting channel: generic over component and face coefficient, ghost-node wall BC) and ThermalPlasmaColumn1D (steady LTE column: single T component, hardcoded kappa, algebraic Dirichlet wall row). Both assemble the SAME two formulas below; they differ only in how they derive face radius/spacing/coefficient and in their wall treatment, which stays local to each caller.

Only the PARABOLIC (gradient) flux lives here, the divergence term of NUMERICS.md section 6.4. The hyperbolic flux of the hydrodynamic sub-step is a Rusanov flux of two reconstructed states and is assembled in EulerFlux.h / HydroStep.cpp instead (NUMERICS.md section 3.3).

From the operator to the two formulas. The grid is node-centred with faces midway between nodes (NUMERICS.md section 2): cell \(i\) spans \([r_{i-1/2}, r_{i+1/2}]\), written \([r_w, r_e]\) below, west being the smaller radius.

Step 1, integrate over the cell rather than over an interval. In cylindrical symmetry the area element per unit axial length is \(2\pi r\,dr\). Integrating the operator over cell \(i\) and dropping the \(2\pi\) common to both sides of the residual,

\[ \int_{r_w}^{r_e} \frac1r\frac{d}{dr}\!\left(r\,\Gamma\,\frac{d\phi}{dr}\right) r\,dr = \int_{r_w}^{r_e} \frac{d}{dr}\!\left(r\,\Gamma\,\frac{d\phi}{dr}\right) dr = \Big[\, r\,\Gamma\,\frac{d\phi}{dr} \Big]_{r_w}^{r_e} \]

The \(1/r\) of the operator and the \(r\) of the area element cancel, so the volume integral collapses to a difference of two FACE quantities. That is the divergence theorem and it is exact; every approximation enters at step 3.

Step 2, divide by the cell volume. The same area element, again without its \(2\pi\), gives the cell volume per unit axial length

\[ V_i = \int_{r_w}^{r_e} r\,dr = \tfrac12\left(r_e^2 - r_w^2\right) \]

which is annularVolume(). The cell average of the operator is therefore \((F_e - F_w)/V_i\) with \(F = r\,\Gamma\,d\phi/dr\) at a face – the (flux_east - flux_west) / volume both callers write.

Step 3, approximate one face flux. The face sits midway between its two nodes by construction, so a central difference ABOUT THE FACE is second-order accurate for the gradient there even on a stretched grid, and needs only the two values straddling it:

\[ F_e = \left. r\,\Gamma\,\frac{d\phi}{dr}\right|_{r_e} \approx r_e\,\Gamma_e\,\frac{\phi_{i+1} - \phi_i}{r_{i+1} - r_i} \]

with \(\Gamma_e\) the face coefficient (the callers' arithmetic mean of the two neighbouring nodes') and \(r_{i+1} - r_i\) the spacing across the face. That is faceFlux(), argument for argument.

What the form buys. Three properties, none of which a finite-difference stencil for the same operator would have:

• Conservation. One face is shared by two cells and the single number computed for it is added to one and subtracted from the other, so summing the discrete operator over all cells telescopes to the two boundary fluxes alone. No \(\phi\) is created or destroyed by the interior discretization, however coarse the grid.
• The axis costs nothing. At node 0 the west face is at \(r_w = 0\), so \(F_w = 0\) whatever the gradient does: the \(d\phi/dr = 0\) symmetry condition is imposed by the geometry, and the axis needs neither a ghost node nor an equation of its own (NUMERICS.md section 7.1).
• A jumping coefficient is admissible. \(\Gamma\) appears only at faces and only multiplying a flux, never differentiated, so a kappa that moves by orders of magnitude across the channel edge enters as a bounded face average rather than an unbounded derivative.

Units follow from the formulas: \(F\) carries one metre more than the physical flux it represents (W/m for conduction, kg/m/s for species diffusion), \(V_i\) is in m^2, and the quotient is the per-volume source the residual rows want (W/m^3, kg/m^3/s).

Definition in file DiffusiveFlux.h.