Compressible viscous flow and entropy-stable viscous fluxes

On the convex half-space $\rho >0$, set $g(\rho ,m)=\mathrm{\mid }m{\mathrm{\mid }}^2\mathrm{/}(2\rho )$. For any direction $(a,b)\in \mathbb{R}\times {\mathbb{R}}^d$, differentiating twice along $(\rho +ta,m+tb)$ gives $\frac{d^2}{d{t}^2}g(\rho +ta,m+tb){\mathrm{\mid }}_{t=0}=\mathrm{\mid }b-(m\mathrm{/}\rho )a{\mathrm{\mid }}^2\mathrm{/}\rho \ge 0$. Thus $g$ is convex on that half-space.

Let $U_0=({\rho }_0,m_0,E_0)$ and $U_1=({\rho }_1,m_1,E_1)$ lie in $\mathcal{G}$, and let $0\le t\le 1$. Then ${\rho }_t=(1-t){\rho }_0+t{\rho }_1>0$ and convexity from step 1.1 gives $g({\rho }_t,m_t)\le (1-t)g({\rho }_0,m_0)+tg({\rho }_1,m_1)$. The strict endpoint inequalities $E_i>g({\rho }_i,m_i)$ imply $E_t>g({\rho }_t,m_t)$, so $U_t\in \mathcal{G}$.

The function $g$ is continuous where $\rho >0$, so the strict inequalities $\rho >0$ and $E-g(\rho ,m)>0$ define an open set there. Step 2.1 proves convexity, and choosing $t=1\mathrm{/}2$ proves the arithmetic-mean assertion.

Set $\epsilon =E-\mathrm{\mid }m{\mathrm{\mid }}^2\mathrm{/}(2\rho )$ and $u=m\mathrm{/}\rho$. Rewriting the entropy gives $\eta =-c_v\rho \log ((\gamma -1)\epsilon )+c_v\gamma \rho \log \rho -s_0\rho$. For an increment $\delta U=(a,b,c)$, its first internal-energy variation is $\delta \epsilon =c-u\cdot b+\mathrm{\mid }u{\mathrm{\mid }}^{2}a\mathrm{/}2$, and its second directional variation is $D^2\epsilon [\delta U,\delta U]=-\mathrm{\mid }b-ua{\mathrm{\mid }}^2\mathrm{/}\rho$. The chain rule therefore yields $D^2\eta [\delta U,\delta U]=c_v(\frac{\gamma }{\rho }{a}^2-\frac{2}{\epsilon }a\delta \epsilon +\frac{\rho }{{\epsilon }^2}(\delta \epsilon {)}^2+\frac{1}{\epsilon }\mathrm{\mid }b-ua{\mathrm{\mid }}^2)$. The first three terms have matrix $\left(\begin{array}{cc}\gamma \mathrm{/}\rho & -1\mathrm{/}\epsilon \\ -1\mathrm{/}\epsilon & \rho \mathrm{/}{\epsilon }^2\end{array}\right)$, whose leading diagonal entry is positive and determinant is $(\gamma -1)\mathrm{/}{\epsilon }^2>0$. The last term is nonnegative. If the full quadratic form is zero, then $a=\delta \epsilon =0$ and $b=0$; the formula for $\delta \epsilon$ then gives $c=0$. Thus the Hessian is positive definite.

Differentiating the same expression for $\eta$ gives $\beta =c_v\rho \mathrm{/}\epsilon =1\mathrm{/}T$ and $v=(c_v\gamma -s-\beta \mathrm{\mid }u{\mathrm{\mid }}^2\mathrm{/}2,\beta u,-\beta )$. Positive definiteness from step 1.1 makes $D_{U}v=D_U^2\eta$ invertible, so the inverse function theorem gives a smooth local inverse near each state. Since $\mathcal{G}$ is convex by [F2], integrating the positive Hessian along the segment between any two distinct states also gives strict convexity of $\eta$ on $\mathcal{G}$.

Fix $j$ and an arbitrary variation $dU=(d\rho ,dm,dE)$, and write $d{u}_j=(d{m}_j-u_{j}d\rho )\mathrm{/}\rho$. Differentiating the flux in [F1] gives $d{F}_j=(u_{j}d\rho +\rho d{u}_j,u_{j}dm+md{u}_j+e_{j}dp,u_j(dE+dp)+(E+p)d{u}_j)$. With $v_m=\beta u$ and $v_E=-\beta$ from step 2.1, the $dp$ terms in $v^{T}d{F}_j-u_j{v}^{T}dU$ cancel, and the coefficient of $d{u}_j$ is $\rho v_{\rho }+\beta \rho \mathrm{\mid }u{\mathrm{\mid }}^2-\beta (E+p)=-\rho s=\eta$, using $E=\epsilon +\rho \mathrm{\mid }u{\mathrm{\mid }}^2\mathrm{/}2$, $p=(\gamma -1)\epsilon$, and $\beta \epsilon =c_v\rho$. Thus $v^{T}d{F}_j=u_j{v}^{T}dU+\eta d{u}_j=d{q}_j$ for every variation, proving ${\mathrm{\nabla }}_U{q}_j=v^{T}D{F}_j$ directly.

Direct substitution of $v$, $F_j$, and $q_j$ gives $v^T{F}_j-q_j=(\rho (c_v\gamma -s-\beta \mathrm{\mid }u{\mathrm{\mid }}^2\mathrm{/}2)+\beta \rho \mathrm{\mid }u{\mathrm{\mid }}^2-\beta E+\rho s)u_j$. Using $E=\epsilon +\rho \mathrm{\mid }u{\mathrm{\mid }}^2\mathrm{/}2$ and $\beta \epsilon =c_v\rho$, the coefficient simplifies to $c_v(\gamma -1)\rho =R\rho$. This proves ${\psi }_j=\rho R{u}_j$ and completes the stated entropy-pair and potential identities.

For any entropy-variable gradient array $Z$, set $G_{k\mathrm{\ell }}=Z_{k,m_{\mathrm{\ell }}}+u_{\mathrm{\ell }}{Z}_{k,E}$ and $c_k=Z_{k,E}$. The formula for $v$ in [F2] gives, when $Z_k={\mathrm{\partial }}_{k}v$, $G_{k\mathrm{\ell }}=T^{-1}{\mathrm{\partial }}_k{u}_{\mathrm{\ell }}$ and $c_k=T^{-2}{\mathrm{\partial }}_{k}T$, since ${\mathrm{\partial }}_k(\beta u_{\mathrm{\ell }})+u_{\mathrm{\ell }}{\mathrm{\partial }}_k(-\beta )=\beta {\mathrm{\partial }}_k{u}_{\mathrm{\ell }}$ and ${\mathrm{\partial }}_k(-1\mathrm{/}T)=T^{-2}{\mathrm{\partial }}_{k}T$.

For a general array $Z$, define ${\tau }_{k\mathrm{\ell }}(Z)=2\mu T{G}^0(Z{)}_{k\mathrm{\ell }}+\zeta T\mathrm{tr}G(Z){\delta }_{k\mathrm{\ell }}$. The displayed bilinear form expands to ${\mathcal{B}}_U(Z,W)={\sum }_{k,\mathrm{\ell }}(W_{k,m_{\mathrm{\ell }}}+u_{\mathrm{\ell }}{W}_{k,E}){\tau }_{k\mathrm{\ell }}(Z)+\kappa T^2{\sum }_k{W}_{k,E}{Z}_{k,E}$. Thus the linear functional in $W$ is pairing against the flux with zero mass row, momentum entries ${\tau }_{k\mathrm{\ell }}(Z)$, and energy entry $u\cdot {\tau }_{k\cdot }(Z)+\kappa T^2{Z}_{k,E}$. By step 1.1 this is precisely the physical viscous flux when $Z=\mathrm{\nabla }v$.

The linear map ${\mathcal{R}}_U$ in the Statement gives ${\mathcal{B}}_U(Z,W)=({\mathcal{R}}_{U}Z)\cdot ({\mathcal{R}}_{U}W)$ term by term. Its adjoint therefore defines the block matrix $K(U)={\mathcal{R}}_U^{\ast }{\mathcal{R}}_U$, which is symmetric and positive semidefinite because $Z^{T}K(U)Z=\mathrm{\mid }{\mathcal{R}}_{U}Z{\mathrm{\mid }}^2\ge 0$. In particular, for $W=Z$ the quadratic form equals $2\mu T\mathrm{\mid }{G}^0(Z){\mathrm{\mid }}^2+\zeta T(\mathrm{tr}G(Z){)}^2+\kappa T^2{\sum }_k{Z}_{k,E}^2\ge 0$, using $T>0$, $\mu >0$, $\zeta \ge 0$, and $\kappa >0$.

Take $Z=\mathrm{\nabla }v$. Then $\mathrm{\nabla }u=TG(Z)$ and $\mathrm{\nabla }T=T^{2}c(Z)$ by step 1.1. Since $\tau$ is symmetric, $\tau :\mathrm{\nabla }u=2\mu T^2\mathrm{\mid }{G}^0(Z){\mathrm{\mid }}^2+\zeta T^2(\mathrm{tr}G(Z){)}^2$. Dividing by $T$ and adding $\kappa \mathrm{\mid }\mathrm{\nabla }T{\mathrm{\mid }}^2\mathrm{/}{T}^2=\kappa T^2\mathrm{\mid }c(Z){\mathrm{\mid }}^2$ gives exactly ${\mathcal{B}}_U(Z,Z)$, proving the production identity and its nonnegativity.

If $U_L=U_R$, then $\mathrm{\Delta }v=0$ and [F1] gives $F_j^{\mathrm{e}\mathrm{c}}(U_L,U_R)=F_j(U_L)$. If $U_L\mathrm{\ne }{U}_R$, convexity [F3] keeps the segment $U(t)=U_L+t(U_R-U_L)$ in $\mathcal{G}$, and $\mathrm{\Delta }{v}^T(U_R-U_L)={\int }_0^1(U_R-U_L{)}^T{D}_U^2\eta (U(t))(U_R-U_L)dt>0$ by [F2]. Hence $\mathrm{\Delta }v\mathrm{\ne }0$, so the denominator in [F1] is nonzero; also $\bar{U}\in \mathcal{G}$. This proves well-definedness and consistency.

Exchanging $U_L$ and $U_R$ leaves $\bar{U}$ fixed and changes the signs of both $\mathrm{\Delta }v$ and $\mathrm{\Delta }{\psi }_j$. Thus the scalar correction coefficient in [F1] changes sign while its vector factor changes sign, leaving the flux unchanged. The diagonal branch is also symmetric.

Fix a state $X\in \mathcal{G}$ and write $U_R=X+h\mathrm{/}2$, $U_L=X-h\mathrm{/}2$ for small $h$. The entropy compatibility in [F2] implies $D_U{\psi }_j(X)=F_j(X{)}^T{D}_{U}v(X)$. Symmetric Taylor expansions therefore give $\mathrm{\Delta }v=D_{U}v(X)h+O(\mathrm{\mid }h{\mathrm{\mid }}^3)$ and $\mathrm{\Delta }{\psi }_j=F_j(X{)}^T{D}_{U}v(X)h+O(\mathrm{\mid }h{\mathrm{\mid }}^3)$, so $N:=\mathrm{\Delta }{\psi }_j-\mathrm{\Delta }{v}^T{F}_j(X)=O(\mathrm{\mid }h{\mathrm{\mid }}^3)$. Positive definiteness of $D_{U}v(X)$ implies $\mathrm{\mid }\mathrm{\Delta }v\mathrm{\mid }\ge c\mathrm{\mid }h\mathrm{\mid }$ in a sufficiently small neighborhood of $h=0$.

When $\mathrm{\Delta }v\mathrm{\ne }0$, dotting [F1] with $\mathrm{\Delta }v$ cancels its denominator and gives $\mathrm{\Delta }{v}^T{F}_j(\bar{U})+\mathrm{\Delta }{\psi }_j-\mathrm{\Delta }{v}^T{F}_j(\bar{U})=\mathrm{\Delta }{\psi }_j$. When $\mathrm{\Delta }v=0$, step 1.1 gives $U_L=U_R$, so both sides are zero.

On that neighborhood, the correction $N\mathrm{\Delta }v\mathrm{/}\mathrm{\mid }\mathrm{\Delta }v{\mathrm{\mid }}^2$ is $O(\mathrm{\mid }h{\mathrm{\mid }}^2)$ and tends to zero on the diagonal. Differentiating the symmetric Taylor remainders gives ${\mathrm{\partial }}_{h}N=O(\mathrm{\mid }h{\mathrm{\mid }}^2)$, ${\mathrm{\partial }}_{X}N=O(\mathrm{\mid }h{\mathrm{\mid }}^3)$, ${\mathrm{\partial }}_h\mathrm{\Delta }v=O(1)$, and ${\mathrm{\partial }}_X\mathrm{\Delta }v=O(\mathrm{\mid }h\mathrm{\mid })$. The quotient rule then shows that the correction derivatives are bounded and tend to zero as $h\to 0$; its diagonal extension is $C^1$. Away from the diagonal the denominator is nonzero by step 1.1 and the formula is smooth. Thus the flux is locally Lipschitz everywhere on $\mathcal{G}\times \mathcal{G}$.

By [F2] and the balance law [F1], ${\mathrm{\partial }}_t\eta +{\sum }_j{\mathrm{\partial }}_j{q}_j=v^T{\sum }_j{\mathrm{\partial }}_j{F}_j^{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}={\sum }_j{\mathrm{\partial }}_j(v^T{F}_j^{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}})-{\sum }_j({\mathrm{\partial }}_{j}v{)}^T{F}_j^{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}$. The second equality is the product rule applied componentwise to each finite-dimensional flux vector.

From the explicit entropy variables in [F2] and viscous flux in [F1], the momentum and stress-work terms cancel pointwise: $v^T{F}_j^{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}=(T^{-1}u)\cdot {\tau }_{j\cdot }-T^{-1}(u\cdot {\tau }_{j\cdot }+\kappa {\mathrm{\partial }}_{j}T)=-\kappa \frac{{\mathrm{\partial }}_{j}T}{T}$. Substitution into step 1.1 and the factorization identity [F3] gives the displayed local balance and its nonpositive right-hand side.

For each direction set $Q_j=H{D}_j$. Then [F1] gives $Q_j^T=-Q_j$ and $Q_j\mathbf{1}=H{D}_j\mathbf{1}=0$. Skew-symmetry also gives ${\mathbf{1}}^T{Q}_j=-(Q_j\mathbf{1}{)}^T=0$.

Multiply the semidiscrete equation by $H_{ii}$ and sum in $i$. For the Euler term, each component of ${\sum }_{i,\mathrm{\ell }}(Q_j{)}_{i\mathrm{\ell }}{F}_j^{\mathrm{e}\mathrm{c}}(U_i,U_{\mathrm{\ell }})$ is zero: swapping $i$ and $\mathrm{\ell }$ changes the sign of $Q_j$ by step 1.1 and leaves the flux unchanged by [F2]. The viscous term sums to ${\mathbf{1}}^{T}H{D}_{j}W={\mathbf{1}}^T{Q}_{j}W=0$ for each nodal flux array $W$. Thus the weighted sum of each conservative component has derivative zero.

By the chain rule, the Euler contribution to the entropy derivative is $-2{\sum }_{j,i,\mathrm{\ell }}(Q_j{)}_{i\mathrm{\ell }}{v}_i^T{F}_j^{\mathrm{e}\mathrm{c}}(U_i,U_{\mathrm{\ell }})$. Skew-symmetry and flux symmetry [F2], after swapping the two indices in one copy of the sum, rewrite this as $-{\sum }_{j,i,\mathrm{\ell }}(Q_j{)}_{i\mathrm{\ell }}(v_i-v_{\mathrm{\ell }}{)}^T{F}_j^{\mathrm{e}\mathrm{c}}(U_i,U_{\mathrm{\ell }})$. The entropy-potential identity [F2] changes the summand to $(Q_j{)}_{i\mathrm{\ell }}({\psi }_j(U_i)-{\psi }_j(U_{\mathrm{\ell }}))$. Both its row sum and column sum vanish by step 1.1, so the total Euler contribution is zero.

For each viscous direction, $v^{T}H{D}_{j}W=-(D_{j}v{)}^{T}HW$ by $D_j^{T}H=-H{D}_j$. Applying this to $W_i={\sum }_k{K}_{jk}(U_i)(D_{k}v{)}_i$ shows that the remaining entropy derivative is exactly the negative sum in the Statement. Each nodal quadratic form is nonnegative by [F3], and each weight $H_{ii}$ is positive; hence the derivative is nonpositive.

Fix an interior face $x_f=x_{i+1\mathrm{/}2}$ and put $U_f=U(x_f)$. Taylor's theorem at $x_f$ gives $U(x_f\pm h\mathrm{/}2)=U_f\pm (h\mathrm{/}2)U^{\mathrm{\prime }}(x_f)+(h^2\mathrm{/}8)U^{\mathrm{\prime }\mathrm{\prime }}(x_f)+O(h^3)$, uniformly in the face index because $U^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}$ is bounded. Averaging the two expansions and using [F1] yields $U_{i+1\mathrm{/}2}=U_f+O(h^2)$.

In [F3], $u=m\mathrm{/}\rho$, $\epsilon =E-\mathrm{\mid }m{\mathrm{\mid }}^2\mathrm{/}(2\rho )$, and $s=c_v\log ((\gamma -1)\epsilon \mathrm{/}{\rho }^{\gamma })+s_0$. These functions and $\beta$ are smooth on $\mathcal{G}$, where $\rho ,\epsilon >0$, so $v\circ U$ is $C^3$ on $[0,L]$ with bounded third derivative. Its centered Taylor expansion gives $(v(U(x_f+h\mathrm{/}2))-v(U(x_f-h\mathrm{/}2)))\mathrm{/}h={\mathrm{\partial }}_{x}v(U_f)+O(h^2)$ uniformly in $i$.

The block $K_{11}(U)$ is smooth on $\mathcal{G}$ by its explicit finite-dimensional factorization [F2], hence has bounded first derivatives on a compact neighborhood of the image of $U$. Step 1.1 gives $K_{11}(U_{i+1\mathrm{/}2})=K_{11}(U_f)+O(h^2)$. Multiplying this by the bounded difference quotient in step 1.2 and using [F1] proves the stated $O(h^2)$ face-flux estimate.