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...
