A Central Differential Flux with High-Order Dissipation for Robust Simulations of Transcritical Flows

TL;DR

This paper introduces the CDHD numerical flux solver, combining central flux and high-order dissipation, to improve the robustness and accuracy of transcritical flow simulations while reducing energy errors.

Contribution

The paper presents a novel primitive-variable flux solver with high-order dissipation that enhances stability and energy conservation in transcritical flow simulations.

Findings

Reduces energy conservation error by two orders of magnitude.

Effectively suppresses spurious pressure oscillations.

Robustly handles smooth transcritical flows and shock waves.

Abstract

The simulation of transcritical flows remains challenging due to strong thermodynamic nonlinearities that induce spurious pressure oscillations in conventional schemes.While primitive-variable formulations offer improved robustness under such conditions, they are always limited by energy conservation errors and the absence of systematic high-order treatments for numerical fluxes. In this paper, we introduce the Central Differential flux with High-Order Dissipation (CDHD), a novel numerical flux solver designed for primitive-variable discretization. This method combines a central flux for advection with a minimal, upwind-biased dissipation term to stabilize the simulation while maintaining formal accuracy. The dissipation term effectively suppresses oscillations and improves stability in transcritical flows. Compared to traditional primitive-variable approaches, CDHD reduces the energy conservation error in two order of magnitude. When incorporated into a hybrid framework with a conservative shock-capturing scheme, the method robustly handles both smooth transcritical phenomena and shock waves. Numerical tests validate the accuracy, stability, and energy-preserving capabilities of CDHD, demonstrating its potential as a reliable tool for complex real-gas flow simulations.

Tables

Table 1. Table 1: Convergence study results for L1-Error.

$N_x$	$dx$	$L_1$-Error($\rho$)	Order	$L_1$-Error($\rho u$)	Order	$L_1$-Error($E$)	Order
32	$3.1250\times {10}^{-2}$	$1.7647\times {10}^1$	N/A	$1.7647\times {10}^3$	N/A	$6.9015\times {10}^6$	N/A
64	$1.5625\times {10}^{-2}$	$8.3227\times {10}^0$	1.0843	$8.3227\times {10}^2$	1.0843	$3.1120\times {10}^6$	1.1490
128	$7.8125\times {10}^{-3}$	$2.6924\times {10}^0$	1.6281	$2.6924\times {10}^2$	1.6281	$9.9529\times {10}^5$	1.6447
256	$3.9062\times {10}^{-3}$	$9.5399\times {10}^{-1}$	1.4969	$9.5399\times {10}^1$	1.4969	$3.4484\times {10}^5$	1.5292
512	$1.9531\times {10}^{-3}$	$2.6181\times {10}^{-1}$	1.8654	$2.6181\times {10}^1$	1.8654	$9.2446\times {10}^4$	1.8993
1024	$9.7656\times {10}^{-4}$	$5.8015\times {10}^{-2}$	2.1740	$5.8015\times {10}^0$	2.1740	$2.0015\times {10}^4$	2.2075

Equations

$$\frac{\partial \mathbf{W}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{W})}{\partial x}=0,$$

$$\frac{\partial \mathbf{V}}{\partial t}+B(\mathbf{V})\frac{\partial \mathbf{V}}{\partial x}=0,$$

$$B(\mathbf{V})=\left[\begin{array}{ccc}u& \rho & 0\\ 0& u& 1/\rho \\ 0& \rho a^2& u\end{array}\right],$$

$$P=\frac{\rho R_{u}T}{M_w-b\rho }+\frac{a\alpha (T){\rho }^2}{M_w^2+2M_{w}b\rho -b^2{\rho }^2},$$

$$a=0.45724\left(\frac{R_u^2{T}_c^2}{P_c}\right),b=0.07780\left(\frac{R_u{T}_c}{P_c}\right).$$

$$e(T,\rho )$$

$$a^2$$

$$\frac{d}{dt}$$

$$+\mathbf{F}\left(\mathbf{W}\left(x_{i+1/2},t\right)\right)-\mathbf{F}\left(\mathbf{W}\left(x_{i-1/2},t\right)\right)=0,$$

$$\frac{d{\mathbf{W}}_i}{dt}=-\frac{1}{\Delta x}\left({\mathbf{F}}_{i+1/2}-{\mathbf{F}}_{i-1/2}\right),$$

$$\mathbf{W}{i}^{n+1}=\mathbf{W}{i}^n-\frac{\Delta t}{\Delta x}\left[{\mathbf{F}}_{i+\frac{1}{2}}-{\mathbf{F}}_{i-\frac{1}{2}}\right].$$

$$\Delta t=C_{\mathrm{cfl}}\frac{\Delta x}{S_{\max }^{(n)}},$$

$$\begin{array}[]{l}S_{0}=\left\{v_{i-2},v_{i-1},v_{i}\right\},\\ S_{1}=\left\{v_{i-1},v_{i},v_{i+1}\right\},\\ S_{2}=\left\{v_{i},v_{i+1},v_{i+2}\right\}.\end{array}$$

$${\beta }_0$$

$${\beta }_1$$

$${\beta }_2$$

$${\alpha }_0=\frac{d_0^{-}}{{\left(ϵ+{\beta }_0\right)}^2},{\alpha }_1=\frac{d_1^{-}}{{\left(ϵ+{\beta }_1\right)}^2},{\alpha }_2=\frac{d_2^{-}}{{\left(ϵ+{\beta }_2\right)}^2},$$

$${\omega }_k=\frac{{\alpha }_k}{{\sum }_{j=0}^2{\alpha }_j}.$$

$$P^{\ast }$$

$$u^{\ast }$$

$${\rho }_{i+\frac{1}{2}}^L$$

$${\rho }_{i+\frac{1}{2}}^R$$

$$s_{i+\frac{1}{2}}^L$$

$$s_{i+\frac{1}{2}}^R$$

$$\frac{P_{i+\frac{1}{2}}}{P_i}>1+\epsilon \text{ and }{s}_{i+\frac{1}{2}}^L<0,$$

$$\frac{P_{i-\frac{1}{2}}}{P_i}>1+\epsilon \text{ and }{s}_{i-\frac{1}{2}}^R<0,$$

$$\frac{\partial \mathbf{V}}{\partial t}+\mathbf{B}(\mathbf{V})\frac{\partial \mathbf{V}}{\partial x}=0,$$

$${\mathbf{L}}_v=\underset{\text{accuracy }}{\underbrace{{\mathbf{L}}_v^{\text{central }}}}+\underset{\text{stability }}{\underbrace{{\mathbf{L}}_v^{\text{diss }}}}.$$

$${{\partial }_x\mathbf{V}|}_i$$

$${\mathbf{L}}_v^{\text{central }}(i)$$

$${\mathbf{H}}_{i+1/2}^{\pm }$$

$$\Psi (s)$$

$${\mathbf{L}}_v^{\mathrm{diss}}(i)=-\frac{1}{\Delta x}\left({\mathbf{H}}_{i+1/2}^{-}+{\mathbf{H}}_{i-1/2}^{+}\right).$$

$${\bar{v}}_i=\frac{1}{\Delta x}{\int }_{x_{i-1/2}}^{x_{i+1/2}}v(\xi )d\xi =v(x_i)+O(\Delta x^2),$$

$$v_{i+1/2}^{\pm }=v(x_{i+1/2})+O(\Delta x^5).$$

$$V_L$$

Full text

1 Introduction

The accurate simulation of fluids in the transcritical regime is essential for the design and analysis of advanced engineering systems[1], such as liquid rocket engines[2, 3, 4, 5], supercritical power cycles[6, 7, 8], and next-generation chip cooling systems[9, 10]. As shown in Figure 1, the regime located near the pseudo-boiling or Widom line[11, 12], is characterized by drastic, liquid-like to gas-like variations in thermodynamic and transport properties. These strong thermodynamic nonlinearities introduce new fluid phenomena [13, 14, 15, 16] but present significant challenges for numerical simulations[17]. The application of widely used conservative schemes often generates severe unphysical pressure oscillations[18] that can render the simulation unstable. Even in ideal gas flows, similar oscillations can be observed in regions with rapid property change, such as material interfaces[19] or in multicomponent flows[20, 21, 22]. Thus, accurately capturing these transcritical phenomena requires not only a precise real-gas Equation of State (EoS), such as the Peng-Robinson[23, 24] model utilized in this paper, but also a numerical scheme that is both accurate and robust.

To address these numerical challenges, numerous attempts have been made to suppress these unphysical oscillations. For low-velocity transcritical flows, researchers have widely adopted the low-Mach-number approximation[25, 26, 27] of the governing equations. For instance, Pasquale[28] use a low-Mach approximation in the direct numerical simulation to fully resolve the interaction between turbulent and transcritical layers in a nitrogen jet. While effective and robust for low-velocity flows, this approximation can introduce non-negligible errors when the flow Mach number exceeds $0.1$. Thus, many new numerical methods have been developed based on non-conservative formulations. The double flux method[29] that is originally designed to reduce the pressure oscillation in multicomponent flows is used for the simulation of transcritical flows[18, 30, 31, 32]. Instead of using the total energy equation, the pressure evolution equation[33, 34, 35, 36] is extensively employed for transcritical flow simulation. Kawai et al. successfully employ a high-order compact differencing scheme to solve the pressure evolution equation for a robust and accurate simulation of transcritical flows[36]. However, these non-conservative methods, including the double-flux and pressure evolution formulations, come at the cost of energy conservation error, which can lead to an incorrect prediction of shock wave locations and strengths. To reduce this conservation error, hybrid methods[37, 38, 22, 39] are also developed by researchers. Bradley Boyd and Dorrin Jarrahbashi [39] propose a hybrid scheme that switches between a double-flux in transcritical region and a traditional fully conservative method for other simulation domains. An emerging research topic in recent years is the development of fully conservative schemes to effectively suppress pressure oscillations during real-fluid simulations[40]. Eric Ching et at. develop a fully conservative discontinuous Galerkin method for supercritical, real-fluid flows [41] that uses an $L^{2}$ projection of primitive variables[42, 43, 44] to reduce pressure oscillations. While this approach maintains stability, it can still result in small deviations from pressure equilibrium. More recently, Fujiwara et al. derive a pressure-equilibrium condition[45] from the governing equations and constructed numerical fluxes to implicitly maintain this condition. Terashima et al. extend this method to real fluids[46], proposing an approximately pressure-equilibrium-preserving scheme with full energy conservation. However, these methods[45, 46, 47] are not easily extended to general transcritical flows.

Despite these advances, the pressure evolution equation remains a promising approach for the robust simulation of transcritical flows. An open source code called CUBENS[48] is recently developed based on the pressure evolution equation for the simulation of real fluid including transcritical flows. However, several key challenges persist. Firstly, discretizing the pressure evolution equation can lead to an incorrect prediction of shock wave locations and strength. Secondly, the energy conservation error should be minimized for industrial applications where heat transfer is critical. Finally, a robust, high-order, and easily extendable numerical flux solver for the primitive variable form of the equation has yet to be developed.

This paper presents a hybrid numerical method designed to address these challenges. We utilize a hybrid framework that switches between a non-conservative, primitive-variable update in smooth regions and a fully conservative update in shock regions. The Primitive Variable Riemann Solver (PVRS)[49] is utilized as a sensor to dynamically switch to a fully conservative scheme (fifth-order WENO[50] with Roe flux[51]) to ensure accurate and robust shock capturing. For the primitive-variable update, we propose a new interface flux solver called Central Differential flux with High-order Dissipation (CDHD) for the pressure evolution equation. This method uses a high-order central flux to provide accuracy for the convection term and augments it with a small, high-order, upwind-biased dissipation term to provide the necessary stabilization for transcritical gradients without degrading the formal order of accuracy. A rigorous error analysis shows that the overall order of accuracy is determined by the central flux, provided that the dissipation term is computed using high-order reconstructed values. Due to its low numerical diffusion, the proposed scheme demonstrates superior energy conservation in smooth regions compared to conventional primitive-variable methods. The hybrid framework ensures robustness for a wide range of flows containing both smooth transcritical phenomena and shocks.

The remainder of this paper is organized as follows: In Section 2, the mathematical model, including the governing equations and thermodynamic closure, is described. In Section 3, we detail the components of the hybrid numerical method. Section 4 provides a formal error analysis justifying the accuracy of the proposed scheme. In Section 5, we present validation results from several one-dimensional (1D) and two-dimensional (2D) test cases to demonstrate the accuracy, stability, and robustness of the proposed scheme. Finally, Section 6 summarizes of the findings.

2 Mathematical Model

2.1 Governing Equations

The governing equations for inviscid transcritical fluid flows are the Euler equations. For clarity, the numerical method is derived and presented for the 1D case. The extension of the proposed method to multiple dimensions is straightforward, as will be demonstrated later with a two-dimensional test case.

The one-dimensional Euler equations can be expressed in a standard conservative form:

[TABLE]

where $\mathbf{W}=[\rho,\rho u,E]^{T}$ represents the vector of conservative variables ($\rho$ for density, $u$ for velocity, $E$ for total energy), and $\mathbf{F}(\mathbf{W})=[\rho u,\rho u^{2}+p,u(E+p)]^{T}$ is the flux. The pressure $p$ is obtained from a real-gas EoS. Alternatively, in quasi-linear primitive form using $\mathbf{V}=[\rho,u,p]^{T}$, the equations become

[TABLE]

where

[TABLE]

and $a$ is the speed of sound, computed via the corresponding EoS of the working fluid.

2.2 Thermodynamic Closure and Equation of State

An accurate description of the thermodynamic nonlinearity for fluids in the transcritical regime is of primary significance. The extreme variations in fluid properties, particularly near the pseudo-boiling or Widom line, require that a thermodynamic model can accurately capture these strong nonlinearities. While the numerical method in this work is developed to be general and compatible with any sufficiently complex EoS, all computations herein will utilize the Peng-Robinson (PR) EoS[23, 24]. The PR-EoS is a cubic EoS that has been shown to reliably predict the thermodynamic properties[52] for hydrocarbons and other fluids in transcritical and supercritical regions. Thus, the PR-EoS should be a suitable choice for the simulation of transcritical fluids.

To be specific, the PR EoS is expressed as

[TABLE]

where $R_{u}$ is the universal gas constant, $M_{w}$ is the molar mass, and $T$ is the absolute temperature. The term $\alpha(T)$ is a temperature-dependent correction factor. The model parameters $a$ and $b$ are determined by the critical pressure $P_{c}$ and critical temperature $T_{c}$ of fluids111This $a$ is the parameter computed by the critical temperature and critical pressure of corresponding fluids for the cubic EoS, rather than speed of sound.:

[TABLE]

The system of Euler equations is closed by defining the thermodynamic properties, such as the specific internal energy $e$ and the speed of sound $c$. These are derived from the chosen EoS using fundamental thermodynamic relations. For detailed derivations for real fluids, refer to the work of Kim et al[52].

The final expression for internal energy $e$ and the speed of sound squared $a^{2}$ are given by[52]

[TABLE]

Here, $e_{0}(T)$ refers to the internal energy of the reference ideal gas state, and the subscript [math] denotes a reference ideal gas state, whose properties can be determined using standard equations, such as the NASA polynomials[18]. The terms $c_{p}$ and $c_{v}$ represent the specific heat capacities at constant pressure and constant volume, while the subscripts $s$ and $T$ denote isentropic and isothermal processes respectively.

3 Numerical Method

This work targets the robust simulation of real-fluid or transcritical compressible flows in which thermodynamic properties vary sharply across a small temperature range near the pseudo-boiling or Wisdom line. During the design of the algorithm for transcritical flows, there are two requirements: (i) strict conservation across shocks and contact discontinuities, and (ii) low spurious pressure-equilibrium (PE) errors when large thermodynamic gradients exist without true discontinuities. To meet both, we employ a hybrid scheme. Cells near the shock waves are updated in a conservative finite volume form using WENO-5 reconstruction[50] and a Roe-type flux solver adapted to the real-gas EoS[51]. Cells in the smooth region, transcritical region are updated in a primitive-variable form with a centered differential flux augmented by a small, upwind-biased dissipation. The primitive variables are reconstructed by WENO-5 to avoid characteristic projections through potentially ill-conditioned real-gas eigenstructures near the critical point. A conservative numerical scheme, such as the finite-volume method, is essential to ensure that the Rankine-Hugoniot jump conditions are satisfied[53]. This guarantees the correct prediction of the shock speed and strength, which is critical for physical accuracy.

Notation: Cell averages are denoted with the subscript $i$. Interface locations are at $x_{i\pm 1/2}$. Superscripts $+$ and $-$ always mean right-biased and left-biased one-sided limits at an interface (e.g., $V^{-}_{1+1/2}$ is the value reconstructed in the left side of cell $i+1/2$). Bold capital $\mathbf{W}$ is the vector of conservative variables, $\mathbf{V}$ is vector of primitive variables, and $B(\mathbf{V})=\partial\mathbf{F}(\mathbf{V})/\partial\mathbf{V}$ is the Jacobian of the inviscid flux written in primitive variables.

3.1 Finite Volume Method

The governing Euler equations in the integral conservative form for a one-dimensional control volume $[x_{i-1/2},x_{i+1/2}]$ are given by[53]

[TABLE]

where $\mathbf{W}$ is the vector of conserved variables and $\mathbf{F}$ is the the corresponding flux vector. On a uniform grid of width $\Delta x$ and defining the cell-averaged value as $\mathbf{W}_{i}(t)=\frac{1}{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}}\mathbf{W}(x,t)dx$, the equation becomes:

[TABLE]

where $\mathbf{F}_{i\pm 1/2}$ represents the numerical flux at the cell interfaces.

Using a forward Euler time integration (or a higher-order method like Runge-Kutta), the semi-discrete form with an explicit update can be written as

[TABLE]

The time step $\Delta t$ is controlled by the Courant-Friedrichs-Lewy (CFL) condition for stability:

[TABLE]

where $S_{\max}^{(n)}=\max_{i}\left\{\left|u_{i}^{n}\right|+a_{i}^{n}\right\}$ is the maximum characteristic wave speed in the domain at time level $n$. The CFL number $C_{\mathrm{cfl}}$ depends on the specific spatial and temporal schemes used.

3.2 WENO Reconstruction

In this work, we employ the WENO-5[50] scheme to reconstruct the primitive variables directly. At each interface $x_{i+1/2}$, we compute right-biased ($v^{+}_{i+1/2}$) and left-biased ($v^{-}_{i+1/2}$) one side limits. In this subsection, superscripts $+$ and $-$ always mean right-biased and left-biased reconstruction at an interface.

To reconstruct a scalar component $v$ for the left biased value $v^{-}_{i+1/2}$, a 5-point stencil $\left\{v_{i-2},v_{i-1},v_{i},v_{i+1},v_{i+2}\right\}$ is used in WENO-5. This can be split into three 3-point sub-stencils:

[TABLE]

Each stencil is associated with a smoothness indicator ($\beta_{0},\beta_{1},\beta_{2}$), which quantifies the variation of variables within each stencil. These indicators are computed as follows:

[TABLE]

The raw weights $\alpha_{k}$ are determined by the smoothness indicators:

[TABLE]

where $\epsilon=1\times 10^{-6}$ is a small value added to avoid division by zero. The constants $d^{-}_{k}$ are the ideal weights that yield fifth-order accuracy in smooth regions. For the left-biased reconstruction $v^{-}_{i+1/2}$, these are $\left\{d_{0}^{-},d_{1}^{-},d_{2}^{-}\right\}=\{0.1,0.6,0.3\}$. The final normalized weights $\omega_{k}$ are then computed:

[TABLE]

The reconstructed value is a weighted sum of polynomial approximations from each sub-stencil.

A similar procedure is followed for the right-biased value at interface $v^{+}_{i+1/2}$, yet using a stencil shifted to the right $\left\{v_{i-1},v_{i},v_{i+1},v_{i+2},v_{i+3}\right\}$ and different ideal weights $\left\{d_{0}^{+},d_{1}^{+},d_{2}^{+}\right\}=\{0.3,0.6,0.1\}$.

3.3 Shock Detection

In this study, shock waves are detected using the PVRS[49]. In PVRS, the Riemann problem is approximately solved by:

[TABLE]

where $\bar{\rho}=\frac{1}{2}\left(\rho_{L}+\rho_{R}\right)$ and $\bar{a}=\frac{1}{2}\left(a_{L}+a_{R}\right)$ are the arithmetically averaged density and speed of sound at the interface, respectively.The propagation speeds of left and right waves are:

[TABLE]

Similar to previous research, the conservative methods are applied at cell $i$ when:

[TABLE]

where the parameter $\varepsilon$ can be selected within the range $(0,0.1)$. In this paper, $\varepsilon=0.05$ is used.

3.4 Conservative Update

For the cells near shock waves, a standard finite-volume WENO-5 scheme with Roe-type real-gas flux is used to update the conservative variable. Primitive variables $\mathbf{V}$ are reconstructed to obtain left ($\mathbf{V}_{L}=\mathbf{V}^{+}$) and right ($\mathbf{V}_{R}=\mathbf{V}^{-}$) states at interfaces, which are then converted to conservative variables $\mathbf{W}_{L}$ and $\mathbf{W}_{R}$. The Roe approximate Riemann solver designed for the real gas is used for the Roe flux. For more details, readers can refer to previous work[51]. It is proved in previous research[37, 38, 39] that, compared with primitive form, the conservative form can maintain an accurate prediction of the speed and strength of shock waves.

3.5 Primitive Update

For cells in the smooth region, we discretize the non-conservative, quasi-linear form of the Euler equations:

[TABLE]

where $\mathbf{V}=[\rho,u,P]^{T}$ is the vector of primitive variables, and $\mathbf{B}(\mathbf{V})$ is the Jacobian matrix of primitive equations. The spatial term $\mathbf{L}_{v}=-\mathbf{B}(\mathbf{V})\partial_{x}\mathbf{V}$ is handled by our proposed CDHD approach. This method split the flux into a central differential term and an upwind dissipative term:

[TABLE]

The central term is approximated using the WENO-5 reconstructed interface state:

[TABLE]

Given that this is a two-point centered difference on staggered locations, its truncation error is $O(\Delta x^{2})$ even though the individual interface states are fifth-order accurate. If desired, a higher order central derivative can be utilized without changing the rest of the formulation.

A small upwind-biased dissipation using path integrals in Dumbser–Osher–Toro Riemann solver (DOTRS)[54, 55] is added for stabilization:

[TABLE]

with $B^{\pm}=\frac{1}{2}(B\pm|B|)$ and a 3-point Gauss-Legendre rule for path integration. The dissipative contribution is then

[TABLE]

In smooth regions, the jump $\mathbf{V}_{R}-\mathbf{V}_{L}=O\left(\Delta x^{5}\right)$ for WENO-5, so that $\mathbf{L}_{v}^{\mathrm{diss}}=O\left(\Delta x^{4}\right)$222For more details of the derivative process about the order of accuracy, please refer to the next section.. Thus the dissipation is subdominant compared to the central term. In fact, the dissipation acts as a gentle upwind stabilization proportional to the spectral radius $|u|+a$, without degrading the formal order of the chosen central derivative.

4 Error Analysis

In this section, we provide a detailed mathematical justification for the behavior of the primitive update in smooth regions. We analyze the truncation error in the central term $\mathbf{L}_{v}^{\text{central}}$, which reveals a second-order accuracy, and the scaling of the state jump $d\mathbf{V}=\mathbf{V}_{R}-\mathbf{V}_{L}$, which leads to DOTRS fluctuations $\mathbf{H}^{\pm}=O(\Delta x^{5})$ and a dissipative term $\mathbf{L}_{v}^{\text{diss}}=O(\Delta x^{4})$. The dissipation magnitude is small enough not to influence the overall order of accuracy, namely the error from the central term dominates. We start by revisiting the key elements of WENO reconstruction and DOTRS, and then proceed with Taylor expansions and error analysis to reveal the behavior. By providing a detail error analysis, one can derive a higher order central flux easily.

4.1 WENO-5 Reconstruction and State Jump Scaling

The WENO-5 scheme, introduced by Jiang and Shu[50], reconstructs interface values from cell-averaged primitives $\mathbf{V}_{i}$ to achieve fifth-order accuracy in smooth regions. For a smooth scalar field $v(x)$, the cell average in cell $i$ (centered at $x_{i}$, with width $\Delta x$) is:

[TABLE]

where the Taylor expansion around $x_{i}$ shows that the approximation error is $O(\Delta x^{2})$.

WENO5 computes left- ($v^{-}_{i+1/2}$) and right-biased ($v^{+}_{i+1/2}$) approximations at interface $x_{i+1/2}$ using a nonlinear weighted average of three candidate polynomials over stencils such that, in smooth regions, the combined reconstruction satisfies

[TABLE]

At interface $i+1/2$ we set:

• •

$V_{L}=V^{+}_{i}$: Right-biased reconstruction from the left cell $i$, approximating $V(x_{i+1/2})+O(\Delta x^{5})$.

• •

$V_{R}=V^{-}_{i+1}$: Left-biased from the right cell $i+1$, also approximating $V(x_{i+1/2})+O(\Delta x^{5})$.

In smooth regions, writing

[TABLE]

with (generally different) constants $c_{L}$ and $c_{R}$, the subtraction of equation (26a) and equation (26b) yields $dV=(c_{R}-c_{L})\Delta x^{5}+O\left(\Delta x^{6}\right)=O\left(\Delta x^{5}\right)$. Lower-order terms are canceled because the scheme is provably fifth-order accurate at the interface. Thus, the jump is

[TABLE]

Near discontinuities the WENO nonlinearity lowers the order and $dV$ can be $O(1)$. Our hybrid method switches computations in such cells to conservative update, and thus the primitive update is used only where the $O\left(\Delta x^{5}\right)$ scaling holds.

4.2 Central Flux Truncation Error

We analyze the truncation error of the central part of the primitive update

[TABLE]

in smooth regions. Two independent $O(\Delta x^{2})$ mechanisms govern the accuracy: (1) the derivative approximation from two midpoint values, and (2) the use of $B\left(\mathbf{V}_{i}\right)$ (a cell averaged value) in place of $B\left(\mathbf{V}(x_{i})\right)$ (an accurate point value).

4.2.1 Central difference from midpoints is second order

Let $V$ be a smooth scalar component. By using exact interface values, one can arrive:

[TABLE]

Thus, the central-difference operator built from the two adjacent midpoints is inherently second order, regardless of how accurately those midpoint values are obtained.

When employing WENO-5 reconstruction in the smooth region, one observes

[TABLE]

so the reconstruction error contributes only

[TABLE]

which is subdominant to the $O(\Delta x^{2})$ from the central difference itself.

4.2.2 Coefficient mismatch of $B\left(\mathbf{V}_{i}\right)$ and $B\left(\mathbf{V}(x_{i})\right)$ is also $O(\Delta x^{2})$

$\mathbf{V}_{i}$ is a cell average $\mathbf{V}_{i}=\mathbf{V}\left(x_{i}\right)+O\left(\Delta x^{2}\right)$ by smoothing $B$:

[TABLE]

By adding and subtracting the exact midpoint values in equation (28), one can obtain

[TABLE]

In the first line of the above equation (33), the second-order central difference produces an $O(\Delta x^{2})$ error. In the second line, $B\left(\mathbf{V}_{i}\right)-B\left(\mathbf{V}\left(x_{i}\right)\right)=O\left(\Delta x^{2}\right)$ times an $O(1)$ gradient also generates an $O(\Delta x^{2})$ error. In the third line, reconstruction defects contribute only $O\left(\Delta x^{4}\right)$ error.Therefore,

[TABLE]

so the overall truncation error of the central term is $O\left(\Delta x^{2}\right)$. Importantly, this second-order behavior is not due to WENO reconstruction (which is 5th order at the interfaces), but to (a) the two-midpoint central differencing and (b) the use of $B(\mathbf{V}_{i})$ evaluated at a cell average.

4.3 Dissipative Term Scaling

Let $d\mathbf{V}=\mathbf{V}_{R}-\mathbf{V}_{L}$ and the straight-line path $\Psi(s)=\mathbf{V}_{L}+sd\mathbf{V},s\in[0,1]$. The DOTRS[54, 55] fluctuations are:

[TABLE]

with eigenvalues $\Lambda=\operatorname{diag}(u-c,u,u+c)$ and $\Psi(s)=\mathbf{V}_{L}+sd\mathbf{V},s\in[0,1]$.

Throughout this section we assume $B$ is $C^{2}$ (twice continuously differentiable) and that the eigen structure of $B(\mathbf{V})$ is smooth along the path $\Psi(s)$ (no eigenvalue sign change), so that $|B|$ and $B^{\pm}=1/2(B\pm|B|)$ are differentiable.

Define $D$ as the Fréchet derivative (directional derivative of a matrix-valued function):

[TABLE]

The small-jump expansion can be derived as follows. Define $f(s)=B^{\pm}(\Psi(s))$. By Taylor expansion in $s$ and the chain rule, one can arrive:

[TABLE]

By integrating the above equation, we obtain

[TABLE]

Multiplying equation (38) by $d\mathbf{V}$ gives

[TABLE]

Hence, $\left\|\mathbf{H}^{\pm}\right\|\leq C_{1}\|d\mathbf{V}\|+C_{2}\|d\mathbf{V}\|^{2}$ , i.e. linear in $d\mathbf{V}$ to the leading order.

By using WENO-5 reconstruction, the reconstructed interface states in smooth region satisfy:

[TABLE]

Thus, $d\mathbf{V}=O\left(\Delta x^{5}\right)$. With the operator norm being controlled by $\|\mathbf{B}\left(\mathbf{V}_{\text{avg }}\right)\|\lesssim C(|u|+c)$, one can arrive

[TABLE]

Thus the dissipation is high-order small and does not limit accuracy when combined with the central term. The total $\mathbf{L}_{v}=\mathbf{L}_{v}^{\text{central}}+\mathbf{L}_{v}^{\text{diss}}$ thus has an overall truncation error $O(\Delta x^{2})$ in smooth regions, remaining spatial accuracy to the second order.

4.4 The Role of DOTRS in First-Order vs. High-Order Schemes

The function of the DOTRS dissipation changes dramatically with the order of the underlying scheme.

4.4.1 Role in First-Order Schemes

In a first-order Godunov-type scheme, states are assumed to be piecewise constant, so $\mathbf{V}_{L}=\mathbf{V}_{i}$ and $\mathbf{V}_{R}=\mathbf{V}_{i+1}$. Jumps across interfaces are large ($d\mathbf{V}=O(1)$). Here, the DOTRS integral acts as the primary Riemann solver, providing the full upwind flux necessary for stability and physical resolution. The resulting spatial discretization $L_{v}(i)=-\frac{1}{\Delta x}(H^{-}_{i+1/2}+H^{+}_{i-1/2})$ is first-order accurate.

4.4.2 Role in This High-Order Scheme

In our high-order method, the primary advection is handled by the accurate central flux term. The role of DOTRS is reduced to that of a targeted stabilizer. Given that high-order reconstruction ensures the interface jump $d\mathbf{V}$ is very small ($O(\Delta x^{5})$) in smooth regions, the DOTRS fluctuation $\mathbf{H}^{\pm}$ also becomes minor. For a small $d\mathbf{V}$, the dissipation behaves similarly to a local Lax-Friedrichs scheme, whereas with a crucial difference that its magnitude is scaled by the high-order jump. The resulting dissipative term in the final update $\mathbf{L}_{v}^{\text{diss}}$ is $O(\Delta x^{4})$. This dissipation is numerically small enough such that it vanishes faster than the error of the second-order central term as the grid is refined. Therefore, it successfully provides stabilization against oscillations without affecting the formal second-order accuracy of the overall scheme.

5 Result and Discussion

In this section, a series of $1\text{D}$ and $2\text{D}$ numerical tests are undertaken to illustrate the accuracy and stability of the proposed scheme. The transcritical nitrogen is adopted as the working fluid in all of the test cases with the thermodynamic properties of fluids obtained by PR EoS. The numerical results are compared with available standard solutions. For a low magnitude dissipation, the WENO-5 reconstruction is used in all cases. The following cases show that the numerical scheme proposed in this paper can achieve an accurate, oscillation free simulation in transcritical regions.

5.1 $1\text{-D}$ Transcritical Advection Cases

5.1.1 Result of CDHD Scheme

A one-dimensional transcritical advection case is considered in a periodic domain $\Omega\in[0,1]$ m. The initial conditions consist of a sharp but continuous temperature gradient centered within the domain, while the velocity and pressure fields are spatially uniform. The setup mimics a droplet in a transcritical environment, which is defined as

[TABLE]

where, $T_{L}=200$ K and $T_{R}=100$ K are the left and right reference temperatures, respectively; $\eta=0.1$ is a smoothing parameter controlling the width of the temperature transition. The density field $\rho(x)$ is computed using the PR EoS.

The proposed CDHD method with a WENO-5 reconstruction is validated against a transcritical convection problem with a known analytical solution. The numerical results at the simulation time $t=0.1\mathrm{\ s}$, corresponding to $10$ periods of convection, are presented in Figure 2. All results are simulated with $512$ grid points.

For the velocity and pressure fields, both the CDHD WENO-5 and the first-order DOTRS schemes demonstrate an excellent agreement with the analytical solution. Both methods achieve a robust simulation, accurately maintaining the constant velocity of $100\mathrm{\ m/s}$ and pressure of $4.0\mathrm{\ MPa}$ across the entire spatial domain. However, significant differences emerge in the temperature and density profiles. The density field exhibits more rapid variations due to the nonlinear thermodynamic effect near the critical point. The CDHD WENO-5 method accurately reproduces both the density and temperature profiles, preserving the sharpness of transitions without introducing spurious oscillations. In contrast, the first-order DOTRS scheme results in noticeable smearing of the interface, limiting its usage in resolving sharp gradients.

Overall, the CDHD method with WENO-5 reconstruction demonstrates superior accuracy and robustness in simulating transcritical flows, particularly in capturing the nonlinearities arising from thermodynamic properties. The method delivers robust performance across all evaluated fields including density, velocity, pressure, and temperature. Our work provides a new strategy for the further development of the high-order scheme for transcritical flows.

5.1.2 Energy Conservation Error

The primitive method has long been criticized for violating the law of energy conservation. However, the results from the proposed CDHD WENO-5 scheme demonstrate substantial improvements in energy conservation. We further examine the temporal evolution of the relative energy conservation error, defined as $\epsilon_{E}=\left|\frac{E(t)-E(0)}{E(0)}\right|$, where $E(t)$ is the integrated total energy at time $t$. Figure 3 presents this error metric over the simulation duration from $t=0$ to $t=0.1\mathrm{s}$ (equivalent to $10$ convection periods) on a semi-logarithmic scale.

The results clearly demonstrate the superior performance of the proposed method in conserving energy. The CDHD WENO-5 scheme maintains a low energy error throughout the simulation, reaching a final error value of approximately $3.95\times 10^{-4}$ at $t=0.1\mathrm{\ s}$. This corresponds to an error of less than $0.04\ \%$, indicating a high degree of energy conservation suitable for long-duration simulations. In contrast, the first-order DOTRS scheme exhibits a significantly larger error that grows to about $10\ \%$ by the final time point. This direct comparison reveals that the energy conservation error of the CDHD WENO-5 scheme is over two orders of magnitude lower than that of the DOTRS scheme at the end of the simulation.

5.1.3 Convergence Test

Table 1 and Figure 4 present the $L_{1}$ errors for $(\rho,\rho u,E)$ under uniform mesh refinement. While the two coarsest grids exhibit pre-asymptotic behavior, the method achieves near-second-order convergence for $N_{x}=256\to 1024$. Specifically, least-squares fits over the three finest meshes yield convergence rates of $r=(2.17,2.17,2.21)$ for $(\rho,\rho u,E)$, which aligns well with our theoretical error analysis.

This second-order accuracy results from the careful balance of several algorithmic components. The WENO-5 reconstruction delivers fifth-order interface states, while the primitive-variable central difference scheme limits the smooth-region accuracy to $O(\Delta x^{2})$. Crucially, the DOTRS term contributes only $O(\Delta x^{4})$ dissipation, effectively stabilizing real-gas gradients without degrading the leading-order accuracy. For applications requiring higher-order accuracy in smooth regions, the overall scheme order can be increased by replacing the central gradient with higher-order central or split-flux differencing.

5.1.4 Role of High-order Dissipation

To isolate the effect of the dissipative term, we repeat the transcritical advection test case using only the central differential flux for $\partial_{x}\mathbf{V}$, omitting the high-order dissipative stabilization. The results in Figures 5 after $600$ time steps with $\text{CFL}=0.8$ clearly demonstrate the necessity of dissipation in the transcritical regime. The small magnitude dissipation proposed in this paper can greatly stabilize the simulation process.

As shown in Figure 5, the central difference scheme without stabilization develops severe numerical instabilities. The density field exhibits spurious oscillations and overshoots near the steep gradient region, with spurious wiggles corrupting the solution throughout the domain. Similarly, the temperature field displays significant nonphysical oscillations that grow rapidly in time. In contrast, the velocity and pressure fields remain relatively stable, staying close to their initial uniform states with only minor perturbations.This selective instability pattern is characteristic of transcritical flow physics, where small perturbations in thermodynamic properties can trigger rapid growth of numerical modes.

These results demonstrate that a small amount of high-order, characteristic-based dissipation is essential for stabilizing the primitive variable update in transcritical flows. The dissipative term proposed in this work provides this stabilization without compromising the leading-order accuracy of the scheme, as evidenced by the convergence analysis presented earlier.

5.2 $1\text{-D}$ Shock Tube Problems

This section demonstrates the capability of the hybrid scheme in this work for accurately capturing the strength and propagation speed of shock waves through a $1\text{-D}$ shock tube test case. We consider a computational domain of $x\in[0,1]$ m with transmissive boundary conditions applied at both ends. All simulations employ a CFL number $0.8$ and are run until $t=5\times 10^{-4}$ s. The initial conditions for the transcritical shock tube problem are specified as a Riemann problem with left and right states:

[TABLE]

This configuration creates a strong pressure and density ratio of the left and right sides in the simulation domain, establishing conditions that span the transcritical regime and challenge numerical methods with steep gradients and rapid thermodynamic variations.

Figure 6 presents the computed density, velocity, pressure, and temperature profiles without using any limiter or filter. The solutions demonstrate the adaptive scheme has the ability to resolve the complex wave structure typical of transcritical shock tubes, including shock waves, contact discontinuities, and expansion regions, while maintaining sharp interfaces without spurious oscillations.

5.3 $2\text{-D}$ Transcritical Advection Cases

To demonstrate the capability of the proposed scheme in multi-dimensional transcritical flows, we consider a 2D droplet convection problem that extends the 1D analysis to more realistic geometries. This test case examines the transport of a subcritical droplet through a supercritical environment under uniform advection.

The computational domain is defined as $\Omega\in[0,1]\times[0,1]$ m with a uniform Cartesian grid of $151\times 151$ cells. The initial state is a circular droplet centered at $(x_{c},y_{c})=(0.5,0.5)$ m with radius $r_{d}=0.15$ m. The temperature field is initialized using a smooth hyperbolic tangent profile to represent the droplet interface:

[TABLE]

where $r=\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$ is the radial distance from the droplet center, $T_{\text{droplet}}=100$ K represents the subcritical liquid phase, $T_{\text{ambient}}=200$ K corresponds to the supercritical gas phase, and $\eta=0.05$ controls the interface thickness. The velocity and pressure fields are initialized as spatially uniform:

[TABLE]

The density field $\rho(x,y)$ is computed from the temperature and pressure using the PR EoS. The simulation runs until $t=0.01\mathrm{s}$ with $\text{CFL}=0.5$, corresponding to one complete convection time. This set up allows the initial state to serve as a analytical solution.

Figure 7 presents the evolution of the transcritical droplet convection, showing the initial and final density distributions along with the final pressure and velocity magnitude fields. The density field evolution reveals that the droplet maintains its circular shape and sharp interface throughout the convection process. The initial density contrast between the subcritical droplet core and the supercritical ambient medium is preserved without significant numerical diffusion. The absence of spurious oscillations near the steep density gradients confirms the efficacy of the WENO-5 reconstruction in handling transcritical interfaces. The pressure and velocity fields exhibit uniformity across the entire domain. The proposed scheme successfully maintains pressure equilibrium while resolving the dramatic thermodynamic variations at the droplet interface. These results confirm that the high-order primitive update with small characteristic-based dissipation transports sharp real-gas gradients cleanly while preserving the constant velocity and pressure fields.

6 Conclusion

We proposed the Central Differential flux with High-Order Dissipation (CDHD) to address spurious oscillations and energy conservation errors in transcritical flow simulations. The method combines a central flux with minimal upwind-biased dissipation, achieving second-order accuracy in smooth regions while substantially improving energy conservation. Embedded in a hybrid framework with a conservative shock-capturing scheme, CDHD demonstrates robust performance across both smooth and shock-containing flows. These results establish CDHD as a reliable and efficient approach for real-gas simulations in engineering and scientific applications.

\printcredits