energy_advect.hpp File Reference

flow — VoF Part II (WO-P23): CONSISTENT rho c_p T transport, driven by the SAME geometric fluxes as the colour advection of the same step (VOF_PLAN §9 item 6). More...

#include <functional>
#include <Kokkos_Core.hpp>
#include <stdexcept>
#include "policy.hpp"
#include "vof/advect_wy.hpp"
#include "vof/momentum_advect.hpp"
#include "vof/phase_change.hpp"

Go to the source code of this file.

Classes
class	peclet::flow::vof::VofEnergyAdvector

Namespaces
namespace	peclet
namespace	peclet::flow
namespace	peclet::flow::vof

Detailed Description

flow — VoF Part II (WO-P23): CONSISTENT rho c_p T transport, driven by the SAME geometric fluxes as the colour advection of the same step (VOF_PLAN §9 item 6).

Why

The energy scalar of rungs P0/P1 was advected by scalar_transport.hpp's Koren TVD flux, i.e. by a DIFFERENT flux from the one that moved the colour. At a mixed cell the two disagree by O(a), and the temperature update then divides a heat content carried with one volume flux by a heat capacity built from another — an artificial heating/cooling of order d(rho c_p), which at the water/steam rho_l c_pl / (rho_g c_pg) ~ 2000 of the P2/P3 rungs is the whole answer. Malan et al. (2021) and Boyd & Ling (2023) both fix it the same way: advect H = (rho c_p) T with the colour's own fluxes and recover T = H / (rho c_p)(C^{n+1}).

The construction (the exact analogue of <tt>momentum_advect.hpp</tt>, one cell-centred unknown)

Per sweep in direction d, with a_f the face Courant number, F_f the geometric LIQUID volume flux Weymouth-Yue already computed for the colour, and Phi_f = rcp_g (a_f - F_f) + rcp_l F_f the heat-capacity flux of that face:

C_i += (F_-  - F_+ )          + flag_i (a_+ - a_-)
H_i += (Phi_- That_- - Phi_+ That_+) + rcp^_i T_i (a_+ - a_-) ,   rcp^ = rcp(flag_i)

That_f is the plain donor-cell (upwind) temperature of the face and flag_i = H(C^n - 1/2) is the SAME frozen dilation indicator the colour uses; T_i = H_i / rcp(C_i^{new}) recovers the temperature. What is actually EVALUATED is the algebraically identical deviation form T_new = T_i + [Phi_-(That_- - T_i) - Phi_+(That_+ - T_i)] / rcp(C^{new}) — see energyUpdate.

The identity that makes it consistent. With rcp(C) = rcp_g + drcp*C and a uniform T == T0: Phi_- - Phi_+ = rcp_g(a_- - a_+) + drcp(F_- - F_+) and the dilation coefficient rcp^_i = rcp_g + drcp*flag_i, so dH = T0 * drcp * dC = T0 * d(rcp) exactly — H = T0 rcp(C) is preserved and the recovery returns T0 at ANY heat-capacity ratio, in floating point, with no tuning knob. Note the dilation term is NOT gated by the flag (unlike the colour's): in a pure gas cell it is rcp_g T_i div u, which is exactly the term that turns the conservative form into the material derivative rho c_p (dT/dt + u.grad T) the energy balance of an evaporating interface is written in. Dropping it is the classic "artificial heating at a divergence source".

Scope

Cut-cell geometry is NOT composed here (enable_phase_change already refuses an immersed solid); advect throws if the advector carries geometry. Boundedness needs no flux clamp: the flux is the colour's own, which Weymouth's proof already bounds on the pressure cell — the clamp momentum_advect.hpp needs exists only because its control volume is HALF-SHIFTED and so carries a colour the planes were not reconstructed from.

Definition in file energy_advect.hpp.