Discretizing the viscous component in 1 - D Navier stokes compressive flow

I've been working on modelling the NS equations in order to simulate shock waves. The equations are set up on the form: $$$$\frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} = S(U)$$$$ where $$$$U = \begin{bmatrix}\rho\\ \rho v\\E\end{bmatrix}, F(U) = \begin{bmatrix}\rho v\\ \rho v^2 + p\\v(E+p)\end{bmatrix}, S(U) = \begin{bmatrix} 0\\ \frac{4}{3} \frac{\partial}{\partial x}\left(\mu\frac{\partial v}{\partial x}\right)\\\frac{4}{3} \frac{\partial}{\partial x}\left(\mu v\frac{\partial v}{\partial x}\right)+ \frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right)\end{bmatrix}$$$$

Wherein E it the total energy $$E = \rho \left(u+\frac{u^2}{2}\right)$$ and the viscosity ($$\mu$$) and thermal conductivity ($$k$$) are calculated based on the given density and temperature. Pressure is calculated with an EOS. I've utilized the FORCE centered flux by Toro to calculate the discretized fluxes $$F_{i-\frac{1}{2}}$$ and $$F_{i+\frac{1}{2}}$$. What I'm unsure of is wheter I have used the right discretization technique for the viscous and conductive contributions. I've simply split the discretization up in this manner: $$$$\frac{4}{3} \frac{\partial}{\partial x}\left(\mu\frac{\partial v}{\partial x}\right) \approx \frac{4}{3} \frac{1}{\Delta x^2}(\mu_{i+\frac{1}{2}}(v_{i+1}-v_i)-\mu_{i-\frac{1}{2}}(v_{i}-v_{i-1}))$$$$

Is this the right way to do it? The flux discretization is based on the paper "MUSTA schemes for multi-dimensional hyperbolic systems: Analysis and improvements" by Titarev and Toro.

• Seems to be correct. The parabolic flux in normal direction can be approximated with the central scheme, as done here. However, the tangential contributions, arising in 2D or 3D, would result in an inconsistency. Sep 22 at 19:16
• Remark: I would write $\frac{\partial F(U,\nabla U)}{\partial x}$ instead of $S(U)$. Sep 22 at 19:32
• To me it seems to be correct. Can I ask you why you have used this particular scheme? Sep 23 at 12:09
• I'm modelling a shock wave where the results will be compared with Non-equillibrium molecular dynamic simulations (NEMD). The initial data in NEMD is discontinuous, so I'm using the FORCE flux in order to avoid oscillations. Sep 23 at 12:37

Consider a FV method as an approximation of the integral conservation law. Starting from the one-dimensional, scalar conservation equation $$$$u_t + f(u,\nabla u)_x =0,$$$$ with the definitions of the step sizes $$\Delta x$$ and $$\Delta t$$ $$$$t_n = n \Delta t,~~~x_i = i \Delta x,~~~x_{i+1/2}=\frac{1}{2}\left(x_i + x_{i+1} \right),$$$$ and $$f$$ representing the hyperbolic and parabolic flux contributions. After integration over $$[x_{i-1/2},~x_{i+1/2}]\times[t_{n},~t_{n+1}]$$ with the abbreviations

$$$$u_i^n :=\frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} u(x,t^n)~dx,~~~~f_{i+1/2}^n :=\frac{1}{\Delta t} \int_{t_{n}}^{t_{n+1}} f(u(x_{i+1/2},t))~dt,$$$$

you can write the differential equation in discrete form as

$$$$u_i^{n+1}=u_i^{n}-\frac{\Delta t}{\Delta x}\left( f_{i+1/2}^n - f_{i-1/2}^n \right).$$$$

Due to the ambiguity of the physical flux at the cell edges, the physical flux $$f$$ is replaced by an arbitrary numerical flux $$g=g(u_L,u_R)$$

$$$$u_i^{n+1}=u_i^{n}-\frac{\Delta t}{\Delta x}\left( g_{i+1/2}^{n} - g_{i-1/2}^{n} \right). \label{eq:diskret1}$$$$

Now for the parabolic flux $$f = \frac{4}{3}\mu\frac{\partial v}{\partial x }$$ you can use central approximations, e.g., the arithmetic mean and a central finite differences

$$$$g(u_L,u_R)=\frac{4}{3} \left[\frac{1}{2}\left(\mu_L + \mu_R\right)\cdot \frac{\displaystyle \left(v_R - v_L\right) }{\displaystyle \Delta x}\right].$$$$

Edit: You also can use a material law for $$\mu$$
$$$$g(u_L,u_R)=\frac{4}{3} \left[ \mu\left( \frac{1}{2}\left(\rho_L + \rho_R\right), \frac{1}{2}\left(T_L + T_R\right)\right) \cdot \frac{\displaystyle \left(v_R - v_L\right) }{\displaystyle \Delta x}\right].$$$$
• Thanks! Follow-up question. Would it be possible to replace the arithmetric mean of the viscosities with $\mu_{1/2} = \mu(rho_{1/2},T_{1/2})$ by using the calculated $rho_{1/2}$ and $T_{1/2}$ from the equation: $$Q^{1/2} = \frac{1}{2}\left((Q_L + Q_R) + \frac{\Delta t}{\Delta x}(F_R -F_L) \right)$$ from the paper I've used? Sep 24 at 7:00
• @Twm1995 Yes, this is even more consistent. Remark: A simple arithmetic mean of $T$ or $\rho$ would be sufficient. Sep 24 at 7:01