PCNSFVHLLC

The derivation of the porous HLLC discretization that follows is based extensively on the material in (Toro, 2009), drawing mostly from chapters 2, 3, and 10. Details pertinent to the MOOSE implementation of the free-flow HLLC discretization can be found in CNSFVHLLCBase.

Solution Properties Across the Contact Wave

HLLC restores the middle contact wave to the HLL formulation. Generalized Riemann Invariants reveal what quantities change or are constant across the wave. We perform the Generalized Riemann Invariants analysis on the porous Euler equations in the following way: we convert the ϵp\epsilon \nabla p term into (Iϵp)pϵ\nabla \cdot\left(\bm{I}\epsilon p\right) - p \nabla \epsilon and ignore the latter term when composing the flux vector F\bm{F} (the term is instead treated as part of a source vector S\bm{S}). Then we define our conserved variable set (for one-dimension for simplicity here)

U=[ρϵρϵuρϵet] \bm{U} = \begin{bmatrix} \rho \epsilon\\ \rho \epsilon u\\ \rho \epsilon e_t \end{bmatrix}(1)

and our flux vector

F=[ρϵuρϵu2+ϵpρϵu(et+pρ)]\bm{F} = \begin{bmatrix} \rho \epsilon u\\ \rho \epsilon u^2 + \epsilon p\\ \rho \epsilon u \left(e_t + \frac{p}{\rho}\right) \end{bmatrix}

and our source vector

S=[0pϵ0]\bm{S} = \begin{bmatrix} 0\\ p \nabla \epsilon\\ 0 \end{bmatrix}

Armed with these definitions we can write the Euler equations succinctly as

Ut+Fx=S. \bm{U}_t + \bm{F}_x = \bm{S}.(2)

We can also write Eq. (2) in a quasi-linear form

Ut+AUx=S\bm{U}_t + \bm{A}\bm{U}_x = \bm{S}

where A\bm{A} is the Jacobian matrix of partial derivatives of F\bm{F} with respect to U\bm{U}. It can be shown for an ideal gas that the eigenvalues of A\bm{A} are

λ1=uc , λ2=u , λ3=u+c\lambda_1 = u - c\ ,\ \lambda_2 = u\ ,\ \lambda_3 = u + c

where cc is the speed of sound in the medium. The corresponding eigenvectors are

K(1)=[1uchtuc] , K(2)=[1u12u2] , K(3)=[1u+cht+uc] \bm{K}^{(1)} = \begin{bmatrix} 1\\ u - c\\ h_t - uc \end{bmatrix} \ ,\ \bm{K}^{(2)} = \begin{bmatrix} 1\\ u\\ \frac{1}{2}u^2 \end{bmatrix} \ ,\ \bm{K}^{(3)} = \begin{bmatrix} 1\\ u + c\\ h_t + uc \end{bmatrix}(3)

where ht=et+p/ρh_t = e_t + p/\rho is the total specific enthalpy. The second eigenvector K(2)K^{(2)} corresponds to the middle contact wave. For a general m×mm \times m system, with variable set:

W=[w1,w2, ... ,wm]\bm{W} = \left[w_1,w_2,\ ...\ ,w_m\right]

and eigenvectors

K(i)=[k1(i),k2(i), ... ,km(i)]\bm{K}^{(i)} = \left[k_1^{(i)},k_2^{(i)},\ ...\ ,k_m^{(i)}\right]

the i-thi\text{-th} Generalized Riemann Invariants are given by the (m1)\left(m - 1\right) ODEs:

dw1k1(i)=dw2k2(i)= ... =dwmkm(i) \frac{dw_1}{k_1^{(i)}} = \frac{dw_2}{k_2^{(i)}} =\ ...\ = \frac{dw_m}{k_m^{(i)}}(4)

Taking our conserved variable set (Eq. (1)) and the contact wave eigenvector from Eq. (3) and substituting into Eq. (4) yields the relations

d(ρϵ)1=d(ρϵu)u=d(ρϵet)12u2\frac{d\left(\rho\epsilon\right)}{1} = \frac{d\left(\rho\epsilon u\right)}{u} = \frac{d\left(\rho\epsilon e_t\right)}{\frac{1}{2}u^2}

These equalities can be algebraically manipulated to yield the following relationships across the contact wave:

u=constant , ϵp=constant u = constant\ ,\ \epsilon p = constant(5)

We will use Eq. (5) when constructing the porous HLLC fluxes below.

Porous HLLC Fluxes

For discontinuous wave solutions over a wave-speed SiS_i associated with the λi\lambda_i characteristic, the Rankine-Hugoniot conditions state that the flux changes according to

ΔF=SiΔU\Delta \bm{F} = S_i\Delta \bm{U}

We can apply the Rankine-Hugoniot conditions to help us establish a discretization for the porous HLLC fluxes. Applying Rankine-Hugoniot conditions over the left wave results in

FL=FL+SL(ULUL) \bm{F}_{*L} = \bm{F}_L + S_L\left(\bm{U}_{*L} - \bm{U}_L\right)(6)

Analogously over the center contact wave:

FR=FL+S(URUL) \bm{F}_{*R} = \bm{F}_{*L} + S_*\left(\bm{U}_{*R} - \bm{U}_{*L}\right)(7)

and the right wave

FR=FR+SR(URUR) \bm{F}_R = \bm{F}_{*R} + S_R\left(\bm{U}_R - \bm{U}_{*R}\right)(8)

The star fluxes can be written FK\bm{F}_{*K} = F(UK)\bm{F}\left(\bm{U}_{*K}\right), where KK denotes LL or RR and preliminarily we will write

UK=[ρKρKuKρKet,K] \bm{U}_{*K} = \begin{bmatrix} \rho_{*K}\\ \rho_{*K} u_{*K}\\ \rho_{*K} e_{t,*K} \end{bmatrix}(9)

Our goal is to eventually express UK\bm{U}_{*K} and consequently FK\bm{F}_{*K} in terms of known left and right quantities, We can leverage the information from Eq. (5) to help us construct star region solutions

{ϵLpL=ϵRpRuL=uR=u \begin{cases} \epsilon_L p_{*L} &= \epsilon_R p_{*R}\\ u_{*L} &= u_{*R} = u_* \end{cases}(10)

Per (Toro, 2009) it is justifiable to select that the middle wave speed be equal to the star region velocity

S=u S_* = u_*(11)

Manipulating the mass, momentum, and energy components of Eq. (6) and Eq. (8), we can construct equations for the star region density, pressure, and total specific energy as functions of the known left and right states and the middle wave speed SS_* (where we have substituted SS_* anyplace we encountered uLu_{*L} or uRu_{*R}). The star region density relationships are given by

ρK=ρKSKuKSKS; \rho_{*K} = \rho_K \frac{S_K - u_K}{S_K - S_*};(12)

the star region pressure relationships are given by

pK=pK+ρK(SKuK)(SuK); p_{*K} = p_K + \rho_K\left(S_K - u_K\right)\left(S_* - u_K\right);(13)

and the total specific energy relationships are given by

et,K=et,K+(SuK)[S+pKρK(SKuK)] e_{t,*K} = e_{t,K} + \left(S_* - u_K\right)\left[S_* + \frac{p_K}{\rho_K\left(S_K - u_K\right)}\right](14)

Substituting Eq. (12), Eq. (13), and Eq. (14) into Eq. (9) and using u=Su_* = S_* from Eq. (11), we arrive at the vector expression for UK\bm{U}_{*K}:

UK=ϵKρK(SKuKSKS)[1Set,K+(SuK)[S+pKρK(SKuK)]]\bm{U}_{*K} = \epsilon_K \rho_K \left(\frac{S_K - u_K}{S_K - S_*}\right) \begin{bmatrix} 1\\ S_*\\ e_{t,K} + \left(S_* - u_K\right)\left[S_* + \frac{p_K}{\rho_K\left(S_K - u_K\right)}\right] \end{bmatrix}

We must now establish an equation for SS_*. Combining Eq. (13) with the pressure information from Eq. (10), we arrive at:

S=ϵRpRϵLpL+ρLϵLuL(SLuL)ρRϵRuR(SRuR)ϵLρL(SLuL)ϵRρR(SRuR) S_* = \frac{\epsilon_R p_R - \epsilon_L p_L + \rho_L \epsilon_L u_L \left(S_L - u_L\right) - \rho_R \epsilon_R u_R \left(S_R - u_R\right)}{\epsilon_L \rho_L \left(S_L - u_L\right) - \epsilon_R \rho_R \left(S_R - u_R\right)}(15)

Left and right wave speeds SLS_L and SRS_R are computed in the same way as for free flow as outlined in CNSFVHLLCBase. With SLS_L and SRS_R the final HLLC flux can be constructed:

FHLLC={FLif 0SLFLif SL0SFRif S0SRFRif 0SR \bm{F}_{HLLC} = \begin{cases} \bm{F}_L & \textrm{if } 0 \leq S_L \\ \bm{F}_{*L} & \textrm{if } S_L \leq 0 \leq S_* \\ \bm{F}_{*R} & \textrm{if } S_* \leq 0 \leq S_R \\ \bm{F}_R & \textrm{if } 0 \geq S_R \end{cases}(16)

Note that although in the derivation above we assumed a one-dimensional setup, the intermediate solution states UK\bm{U}_{*K} can be generalized to three-dimensions in a way analogously to the multidimensional free-flow intermediate solutions expressed in CNSFVHLLCBase. Indeed, the porous intermediate states can be simply expressed as

UK,porous=ϵKUK,free\bm{U}_{*K,porous} = \epsilon_K \bm{U}_{*K,free}

and with SS_* expressed according to Eq. (15).

The mass component of Eq. (16) is implemented in PCNSFVMassHLLC, momentum in PCNSFVMomentumHLLC, and fluid energy in PCNSFVFluidEnergyHLLC.

References

  1. Eleuterio F Toro. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media, 3rd edition, 2009.[BibTeX]