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

Question

I've been working on modelling the NS equations in order to simulate shock waves. The equations are set up on the form: \begin{equation} \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} = S(U) \end{equation} where \begin{equation} 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} \end{equation}

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: \begin{equation} \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})) \end{equation}

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.

Answer 1

Consider a FV method as an approximation of the integral conservation law. Starting from the one-dimensional, scalar conservation equation \begin{equation} u_t + f(u,\nabla u)_x =0, \end{equation} with the definitions of the step sizes $\Delta x$ and $\Delta t$ \begin{equation} t_n = n \Delta t,~~~x_i = i \Delta x,~~~x_{i+1/2}=\frac{1}{2}\left(x_i + x_{i+1} \right), \end{equation} 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

\begin{equation} 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, \end{equation}

you can write the differential equation in discrete form as

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

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)$

\begin{equation} 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} \end{equation}

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

\begin{equation} 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]. \end{equation}

which should recover your approximation.

Edit: You also can use a material law for $\mu$

\begin{equation} 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]. \end{equation}