Structure-Preserving Finite Volume Schemes

• The paper advances a finite volume scheme that exactly preserves discrete conservation laws and curl-free properties in continuum systems.
• It employs mimetic discretization with discrete gradient and curl operators on unstructured Voronoi/Delaunay meshes to enforce thermodynamic and geometric constraints.
• Numerical tests confirm robust asymptotic-preserving behavior, accurate Fourier heat conduction limits, and strict entropy consistency.

A structure-preserving finite volume scheme is a numerical framework designed to exactly maintain discrete analogs of key invariant properties—such as conservation laws, stationary involutions, thermodynamic compatibility, and asymptotic limits—of continuum partial differential equations. These schemes are distinguished by their mimetic discretization of differential operators and fluxes, designed to enforce, at the mesh and time-discretization level, the geometric, thermodynamic, and physical structure of the underlying PDE system. Recent developments focus on hyperbolic thermodynamically compatible systems, such as compressible heat-conducting flow models formulated in the Godunov–Romenski SHTC framework. Emphasis is placed on (i) exact preservation of discrete curl-free involutions, (ii) construction of compatible discrete grad–curl operators on unstructured Voronoi grids, (iii) design of semi-discrete cell solvers for auxiliary variables, and (iv) discrete enforcement of energy conservation and entropy production, including rigorous asymptotic preserving (AP) limits toward classical physical laws such as Fourier heat conduction.

1. Hyperbolic System and Thermodynamical Constraints

The governing system augments compressible Euler equations with a hyperbolic heat conduction law using the thermal impulse vector $\mathbf{j}(\mathbf{x},t)$. The state vector is

$$\tilde{w} = (\rho,\;\rho\mathbf{u},\;E,\;\mathbf{j})^T,$$

with specific energy $E$ including thermal and kinetic contributions: $$E = \frac{1}{2}\rho\|\mathbf{u}\|^2 + \rho\,\varepsilon(\rho,\eta) + \frac{1}{2}\|\mathbf{j}\|^2,$$ and the entropy density $s = \rho\,\eta$.

The hyperbolic evolution equations are: \begin{align*} &\partial_t\rho + \nabla\cdot(\rho\mathbf{u}) = 0,\ &\partial_t(\rho\mathbf{u}) + \nabla\cdot\Bigl(\rho\mathbf{u}\otimes\mathbf{u} + \Bigl(p+\frac{1}{2}(\rho\,\alpha'(\rho)-1)|\mathbf{j}|^2\Bigr)\mathbf{I}+\mathbf{j}\otimes\mathbf{j}\Bigr)=\mathbf{0},\ &\partial_t\mathbf{j} + \nabla\cdot(\mathbf{j}\otimes\mathbf{u}) + \nabla\bigl(\theta(\rho,\eta)\bigr) + \nabla\times(\mathbf{j}\times\mathbf{u}) = -\frac{1}{\tau}\mathbf{j},\ &\partial_t E + \nabla\cdot\Bigl((E+p)\mathbf{u}+(\mathbf{j}\cdot\mathbf{u})\mathbf{j}+\mathbf{q}\Bigr)=0, \end{align*} supplemented with the entropy law: $$\partial_t(\rho\eta)+\nabla\cdot(\rho\eta\,\mathbf{u}+\mathbf{j}/\theta) = \frac{1}{\theta\,\tau}\|\mathbf{j}\|^2\,\ge0.$$ A crucial involution is the curl-free propagation for $\mathbf{j}$: $$\nabla\times\mathbf{j}(t=0)=\mathbf{0}\;\Longrightarrow\;\nabla\times\mathbf{j}(t)=\mathbf{0}\;\forall\,t>0.$$ In the stiff relaxation limit ($\tau\to0$), the system is consistent with Fourier’s law: $$\mathbf{q} = -\tau\,\theta\,\nabla\eta \approx -K\,\nabla T.$$

2. Mimetic Discrete Curl–Gradient Operators on Voronoi Grids

Let $\{\omega_c\}$ denote primal Voronoi cells and $\{\omega_p\}$ dual Delaunay triangles. Discrete operators are constructed to enforce compatibility:

• Discrete gradient (cell $\to$ node):

$$(\nabla_h^c\phi)_p = \frac{1}{|\omega_p|} \sum_{c\in\mathcal{C}(p)}l_{cp}\,\mathbf{n}_{cp}\,\phi_c$$

• Discrete curl (node $\to$ cell):

$$(\nabla_h^p\times\mathbf{j})_c = -\frac{1}{|\omega_c|} \sum_{p\in\mathcal{P}(c)}l_{cp}\,\mathbf{t}_{cp}\cdot\mathbf{j}_p$$

where $\mathbf{t}_{cp}$ is the local tangent vector and $l_{cp}\mathbf{n}_{cp}$ is a dual-edge normal.

Compatibility: Direct algebra yields the discrete mimetic identity: $$\nabla_h^p\times(\nabla_h^c\phi) = 0\quad\forall\;\{\phi_c\}$$ ensuring that curl-free initial data for $\mathbf{j}_p$ remains curl-free under the update. This discrete preservation is nontrivial on general Voronoi/Delaunay meshes and is a distinguishing feature of the method.

3. Cell Solver and Entropy Production for the Thermal Impulse

The update for the dual-cell $\mathbf{j}_p$ is

$$\frac{d\mathbf{j}_p}{dt} = -\frac1{|\omega_p|}\sum_{c\in\mathcal{C}(p)}l_{cp}\,\mathbf{n}_{cp}\,\varphi_c -\frac1{|\omega_p|}\sum_{c\in\mathcal{C}(p)}|\omega_{cp}| (\nabla_h^p\times\mathbf{j})_c\times\mathbf{u}_c -\frac1{\tau}\mathbf{j}_p,$$

where $\varphi_c = \mathbf{j}_c\cdot\mathbf{u}_c+\theta_c$ is computed with a nodal-Lagrangian-type jump relation: $$\varphi_p-\varphi_c = z_p(\mathbf{j}_c-\mathbf{j}_p)\cdot\mathbf{n}_{pc},\qquad z_p = \max|\lambda(\mathbf{w}_p)|.$$ This generates a symmetric positive-definite ($2\times2$) system for the primal-cell $\mathbf{j}_c$, mimicking subcell force balances from classical staggered hydrodynamics.

Entropy production and energy conservation are enforced via compatible flux modifications, drawing on Abgrall-type thermodynamic corrections: $$\tilde h_f\cdot n_{fc} = \hat h_f\cdot n_{fc} - \alpha_f(p_d-p_c),$$ where $\alpha_f$ is chosen so that a discrete entropy balance is satisfied exactly for smooth solutions. Local non-negativity of entropy increment in each dual cell is obtained by pairing the update for $\mathbf{j}_p$ with its entropy-conjugate $\beta_p$: the resulting terms yield negative-definite Hessian dissipation and path-integral entropy production, ensuring strict thermodynamic consistency.

4. Asymptotic-Preserving and Subflux Transfer Mechanisms

IMEX time discretization of the relaxation term delivers the asymptotic-preserving (AP) property: as $\tau\to0$, the discrete scheme recovers

$$\mathbf{j}_p^n = -\tau\,\nabla_h^c\theta^n + \mathcal{O}(\tau^2)$$

yielding the correct Fourier limit for heat conduction.

To integrate the dual variables into primal cell updates for entropy correction, each Delaunay triangle is subdivided into subcells, building a subcell finite-volume balance: $$\frac{d\mathbf{j}_{mp}}{dt} = -\frac{1}{|\omega_{mp}|}\left(\sum \hat F^j_{fm} + \sum \hat B^j_{fm}\right),$$ with equivalence constraints solved using the graph Laplacian and its pseudo-inverse. This enables consistent application of entropy-compatible flux corrections across the entire scheme.

5. Numerical Validation and Structure-Preservation Results

A suite of tests confirms the effectivity and structure-preserving attributes:

• Manufactured solutions: First-order convergence in all variables on Voronoi grids, with machine-zero discrete curl of $\mathbf{j}$.
• Riemann problems: Preservation of $\nabla_h\times\mathbf{j}=0$ and sharp wave capturing on quasi-one-dimensional grids.
• Vortex and dissipation: The method maintains curl-freeness and semi-discrete entropy law to machine tolerance; non-compatible variants lose structure.
• Cylindrical explosion: Correct asymptotic Fourier relaxation in the heat-conducting limit, matching reference solutions.

These results confirm that discrete curl-involution, AP behaviour for Fourier conduction, and thermodynamic compatibility are enforced in practice. The scheme is robust and implementable with provable structure-preservation properties (Boscheri et al., 28 Jul 2025).

6. Significance and Research Trajectory

Contemporary structure-preserving finite volume frameworks on unstructured and Voronoi meshes align with the broader movement toward mimetic discretization in hyperbolic thermomechanics. They provide critical robustness for multi-physics applications—compressible flow, heat conduction, elasticity—where geometric involutions and thermodynamic laws must be satisfied even as the mesh and temporal discretization evolve. The distinguishing contribution lies in the exact discrete preservation of continuous PDE constraints, notably curl-free involutions, energy conservation, entropy production, and Fourier limits.

Extensions to higher-order accuracy, alternative physical regimes (e.g., magnetohydrodynamics, nonlinear acoustics), and general grid topologies are ongoing themes. Structure-preserving FV schemes are central to advancing reliable and physically consistent simulation of multiscale continuum phenomena.