Fix for FD WENO method for multi-component compressible flows

I'm solving two-dimensional four-component compressible Navier-Stokes equations with finite-difference WENO approach. The equations are pretty standard:

$$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0;$$

$$U = \left( \begin{array} {c} \rho \\ \rho u \\ \rho v \\ E \\ \rho_1 \\ \rho_2 \\ \rho_3 \end{array} \right), \;\; F = \left( \begin{array} {c} \rho u \\ \rho u^2 + p -\tau_{xx} \\ \rho uv -\tau_{xy} \\ (E + p)u -u\tau_{xx} -v\tau_{xy} - k \frac{\partial T}{\partial x} \\ \rho_1 u \\ \rho_2 u \\ \rho_3 u \end{array} \right), \;\; G = \left( \begin{array} {c} \rho v \\ \rho uv -\tau_{xy} \\ \rho v^2 + p -\tau_{yy} \\ (E + p)v -u\tau_{xy} -v\tau_{yy} - k \frac{\partial T}{\partial y} \\ \rho_1 v \\ \rho_2 v \\ \rho_3 v \end{array} \right)$$

where $$\rho = \rho_1+\rho_2+\rho_3+\rho_4$$.

The following equations of state are used:

$$E = \frac{p}{\gamma-1}+\rho\frac{u^2+v^2}{2},\quad p=\sum_{i=1}^{4}\rho_i\frac{R_0}{\mu_i}T$$ Where $$R_0$$ is universal gas constant and $$\mu_i$$ are gas component molar masses. The second of these EOS is used only to calculate $$T$$ from $$P,\rho_i$$. For now, heat capacity ratio is considered constant and same for all gas components: $$\gamma = 1.4$$.

Here are the Jacobian matrix for inviscid part of $$F$$ and its left and right eigenvectors:

$$J_f = \frac{\partial F(U)}{\partial U} = \left( \begin{array}{ccccccc} 0 & 1 & 0 & 0 & 0 & 0 & 0\\ (\gamma-1)q-u^2 & (3-\gamma)u & (1-\gamma)v & \gamma-1 & 0 & 0 & 0\\ -uv & v & u & 0 & 0 & 0 & 0\\ \Big((\gamma-1)q-h\Big)u & h-(\gamma-1)u^2 & (1-\gamma)uv & \gamma u & 0 & 0 & 0\\ -\frac{\rho_1}{\rho}u & \frac{\rho_1}{\rho} & 0& 0 & u& 0& 0\\ -\frac{\rho_2}{\rho}u & \frac{\rho_2}{\rho} & 0& 0 & 0& u& 0\\ -\frac{\rho_3}{\rho}u & \frac{\rho_3}{\rho} & 0& 0 & 0& 0& u \end{array}\right)$$ $$R_f = \left( \begin{array}{ccccccc} 1 & 1 & 1 & 0 & 0 & 0 & 0\\ u-a & u & u+a & 0 & 0 & 0 & 0\\ v & v & v & -1 & 0 & 0 & 0\\ h-au & q & h+au & -v & 0 & 0 & 0\\ \frac{\rho_1}{\rho}& \frac{\rho_1}{\rho}& \frac{\rho_1}{\rho}&0& 1 & 0 & 0\\ \frac{\rho_2}{\rho}& \frac{\rho_2}{\rho}& \frac{\rho_2}{\rho}&0& 0 & 1 & 0\\ \frac{\rho_3}{\rho}& \frac{\rho_3}{\rho}& \frac{\rho_3}{\rho}&0& 0 & 0 & 1 \end{array}\right)$$ $$L_f = R^{-1}_f = \left( \begin{array}{ccccccc} \frac{(\gamma-1)q+au}{2a^2} & \frac{(1-\gamma)u-a}{2a^2} & \frac{(1-\gamma)v}{2a^2} & \frac{(\gamma-1)}{2a^2} & 0 & 0 & 0\\ \frac{a^2-(\gamma-1)q}{a^2} & \frac{(\gamma-1)u}{a^2} & \frac{(\gamma-1)v}{a^2} & \frac{(1-\gamma)}{a^2} & 0 & 0 & 0\\ \frac{(\gamma-1)q-au}{2a^2} & \frac{(1-\gamma)u+a}{2a^2} & \frac{(1-\gamma)v}{2a^2} & \frac{(\gamma-1)}{2a^2} & 0 & 0 & 0\\ v & 0 & -1 & 0 & 0 & 0 & 0\\ -\frac{\rho_1}{\rho} & 0 & 0 & 0 & 1 & 0 & 0\\ -\frac{\rho_2}{\rho} & 0 & 0 & 0 & 0 & 1 & 0\\ -\frac{\rho_3}{\rho} & 0 & 0 & 0 & 0 & 0 & 1\\ \end{array}\right)$$ Where $$q = \frac{u^2+v^2}{2},\; a=\sqrt{\frac{\gamma p}{\rho}}, \; h=\frac{a^2}{\gamma-1} + q = \frac{\gamma}{\gamma-1}\frac{p}{\rho} + q$$.

The values for the matrices are Roe-averaged between neighbor grid points ($$L$$ and $$R$$) in the following way: $$u = \frac{\sqrt{\rho_L}u_L + \sqrt{\rho_R}u_R}{\sqrt{\rho_L}+\sqrt{\rho_R}}, \quad v = \frac{\sqrt{\rho_L}v_L + \sqrt{\rho_R}v_R}{\sqrt{\rho_L}+\sqrt{\rho_R}}, \quad \frac{\rho_i}{\rho} = \frac{\sqrt{\rho_L}\frac{\rho_{iL}}{\rho_L} + \sqrt{\rho_R}\frac{\rho_{iR}}{\rho_R}}{\sqrt{\rho_L}+\sqrt{\rho_R}},$$ $$E_L = \frac{p_L}{\gamma-1}+\rho_L \frac{u_L^2+v_L^2}{2}, \quad E_R = \frac{p_R}{\gamma-1}+\rho_R \frac{u_R^2+v_R^2}{2}, \quad h = \frac{\sqrt{\rho_L}\frac{\left(E_L+p_L\right)}{\rho_L} + \sqrt{\rho_R}\frac{\left(E_R+p_R\right)}{\rho_R}}{\sqrt{\rho_L}+\sqrt{\rho_R}}, \quad$$ $$q = \frac{\sqrt{\rho_L}q_L + \sqrt{\rho_R}q_R}{\sqrt{\rho_L}+\sqrt{\rho_R}}, \quad a^2 = (\gamma-1)(h-q)$$

Left eigenvector matrix $$L_f$$ is used to project conservative values $$U$$ and their fluxes $$F$$ into local characteristic space.

Then resulting characteristic variables $$\tilde{U}$$ and fluxes $$\tilde{F}$$ are split into positive and negative parts using Roe-LLF splitting: $$\tilde{F}^+=\tilde{F}+\lambda^{max}_j\tilde{U}, \quad \tilde{F}^-=\tilde{F}-\lambda^{max}_j\tilde{U}$$ where $$\lambda^{max}_j$$ is the maximum characteristic speed ($$|u-a|, |u|$$, or $$|u+a|$$, depending on equation). The maximum is found over only 2 grid points ($$R$$ and $$L$$) points individually for each equation.

Then these split fluxes are WENO-reconstructed individually over left- and right-biased stencils, summed up ($$\tilde{F} = (\tilde{F}^+ + \tilde{F}^-)/2$$), transmuted back into original space using $$R_f$$ matrix and used as numerical flux in current Runge-Kutta step.

This approach worked well for single-component Euler/NS equations or for heterogeneous 3-component mixture with same molar masses. But if one of the components has different molar mass, then there are non-physical oscillations near the interface between different gases (for example, only $$\rho_1 \ne 0$$ on one side and only $$\rho_4 \ne 0$$ on the other side).

During my research of high-order numerical methods for compressible multi-component flows, I got overwhelmed by the number of different models and approaches: overestimated equations with $$\frac{1}{\gamma-1}$$ advection, overestimated system with $$\rho$$ and all $$\rho_i$$ components being solved and then corrected with diffusion fluxes, primitive and conservative values formulations, subcell resolution etc...

It seems that there should be not very complex "fix" for the present approach, but I couldn't find it yet. What is the critical component to maintain non-oscillating behavior near gas interfaces? I suspect several things:

1. Jacobian and/or eigenvector matrices should be different. It seems unlikely, since all additional terms I found in literature seem to zero out in the case of constant heat capacities.

2. Roe-averaging should be more complex. Probably, the averaged pressure and/or speed of sound should be calculated differently. Although, when I tried to implement more complex formulas (from Liou, M.S., Van Leer, B., Shuen, J.S. (1990). Splitting of inviscid fluxes for real gases. Journal of Computational Physics, 87(1), 1-24), they resulted in exactly the same speed of sound as formulas above, even at the interface.

3. Proper multi-component equation of state should be used. Maybe, my "naive" approach with constant $$\gamma$$ is incorrect, and proper model with individual molar masses, heat capacities and heat of formations should be used, so that $$\gamma$$ would derive from them and vary in each grid point.

4. Finite difference approach is inherently inconsistent with gas interfaces. I've read somewhere few years ago, that only finite-volume approach is feasible for such flows. It seems unlikely that no solution was found for FD approach though...

• Have you tried global splitting, i.e., estimate $\lambda_{max}$ globally ? Dec 13 '20 at 3:45
• @cfdlab This will not work. The problem is inherently in the equation system. Dec 13 '20 at 22:17
• Yes, I can confirm that global splitting doesn't solve the problem. Dec 15 '20 at 6:24
• Can you show some sample results ? Dec 16 '20 at 9:38

Remark:

With constant/same adiabatic coefficients the species equations are only passive scalars and you basically solve only the standard Euler equations. This is the reason why the oscillations do not occur.

You are facing the problem of spurious pressure und velocity oscillations arising with real equations of state and multicomponent flows (see works of Remi Abgrall and Richard Saurel). The problem is not so easy as you think and is the result of the missing mechanical equilibrium in the Euler equations when you have jumps in the material properties, especially at the cell faces. All methods available in the literature are strictly non-conservative. Moreover, if you want to implement one of these methods, e.g. the double flux method, then there is no free lunch. All these methods have their drawbacks.

Note that none of your four points is related to the actual problem and will even produce additional problems.

1. The Jacobian will change and you have to consider additional thermodynamic gradients. Your Jacobian won't work.

2. You have to enforce all Roe properties, e.g. $$F(U_R) − F(U_L) = \hat{A}(U_R − U_L)$$.

3. See my first remark.

4. Not necessarily.