# Understanding Davis artificial viscosity

I'm solving 2D Euler's equations in Cartesian coordinates: $$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0, \qquad U = \left( \begin{array}{c} \rho \\ \rho u \\ \rho v \\ e \end{array} \right), \quad F = \left( \begin{array}{c} \rho u \\ p+\rho u^2 \\ \rho uv \\ (e+p)u \end{array} \right), \quad G = \left( \begin{array}{c} \rho v \\ \rho uv \\ p+\rho v^2 \\ (e+p)v \end{array} \right)$$

using finite difference MacCormack method with Davis artificial viscosity which dampens non-physical oscillations. The MacCormack method is pretty simple and not shown here. Davis viscosity works good in my simulations, but is quite difficult; I'm trying to fully understand it in order to apply it cylindrical coordinates later. I hope that some of you may help me with unclear parts of it (shown in bold font below).

Viscous terms are added after 2nd MacCormack step as follows (formulated in 1D for simplicity): $$U_j^{n+1}=U_j^{(m)} + \left( D_{j+\frac{1}{2}}-D_{j-\frac{1}{2}} \right),$$ where $U_j^{(m)}$ denote field values (at $j$-th coordinate point and $(n+1)$-th time layer) calculated by MacCormack scheme as usual, and $D_{j+\frac{1}{2}}, D_{j-\frac{1}{2}}$, are forward and backward differences of U with non-linear coefficients: $$D_{j+\frac{1}{2}} := \frac{1}{2} C(\nu_j) \left( 2 - \phi(r_j^+) - \phi(r_{j+1}^-) \right) \Delta U_{j+\frac{1}{2}}, \quad \Delta U_{j+\frac{1}{2}} = U_{j+1} - U_j \\ D_{j-\frac{1}{2}} := \frac{1}{2} C(\nu_{j-1}) \left( 2 - \phi(r_{j-1}^+) - \phi(r_{j}^-) \right) \Delta U_{j-\frac{1}{2}}, \quad \Delta U_{j-\frac{1}{2}} = U_{j} - U_{j-1}$$

Where $C(\nu)$ is a limited coefficient based on Courant number: $$C(\nu) := \begin{cases} \nu(1-\nu) & \;\; \text{if} \;\; \nu \le 0.5 \\ 0.25 & \;\; \text{otherwise} \end{cases}, \qquad \nu = \rho(A(U_j)) \frac{\Delta t}{\Delta x},$$

$A(U)$ is the Jacobian matrix of vector $F$, $\rho(A(U))$ is it's spectral radius (for Euler's equations it equals $\text{max}(|u|+a, |u|-a)$, where $a$ is the speed of sound).

$\phi(r)$ is simple limiting function: $\begin{cases} 0, & r < 0, \\ 2r, & 0 \le r \le 0.5, \\ 1, & r >0.5\\ \end{cases}$

And finally $r_j^+$,$r_j^-$ are slope ratios: $$r_j^+ := \frac{\left( \Delta U_{j-\frac{1}{2}}, \Delta U_{j+\frac{1}{2}}\right)}{\left( \Delta U_{j+\frac{1}{2}}, \Delta U_{j+\frac{1}{2}}\right)}, \qquad r_j^- := \frac{\left( \Delta U_{j-\frac{1}{2}}, \Delta U_{j+\frac{1}{2}}\right)}{\left( \Delta U_{j-\frac{1}{2}}, \Delta U_{j-\frac{1}{2}}\right)}$$ Where $\left( \cdot , \cdot \right)$ is scalar product. These ratios are positive in monotone areas and negative in non-monotone areas, so that in general coefficients like $\left( 2 - \phi(r_j^+) - \phi(r_{j+1}^-) \right)$ are equal 2 in non-monotone areas; about 1 in monotone, but "curvy" areas; and equal 0 in monotone (and rectilinear) areas.

If both these ($2-\phi-\phi$) coefficients are equal 2, and both $C(\nu)$ are equal 0.25 we have an effective second derivative of the field value: $$D_{j+\frac{1}{2}}-D_{j-\frac{1}{2}} = 0.25 \left( U_{j+1} - 2U_j + U_{j-1} \right).$$

But where is the division by $\Delta x^2$? If the viscosity would be simply the 2nd derivative, the scheme would look like $$U_j^{n+1}=U_j^{(m)} + \frac{\Delta t}{\Delta x^2} \left( U_{j+1} - 2U_j + U_{j-1} \right).$$ Is $C(\nu)$ effectively equals $\frac{\Delta t}{\Delta x^2}$?

And, what if $D_{j+\frac{1}{2}}$ is non-zero and $D_{j-\frac{1}{2}}$ is zero? Does it still yield an effective 2nd derivative?

• $C$ and $D$ just stand for the left and right coefficients for the variation. You are probably confusing them with convection and diffusion. Refer to (2.6) in this paper epubs.siam.org/doi/pdf/10.1137/0908002 Aug 8 '17 at 4:16
• @gpavanb In that paper (which is original paper of the method) the letter D is used again with the different meaning, see the unnumbered formula right below (4.14). I'm using that second meaning here. Aug 9 '17 at 5:24