# How to apply central difference to viscous fluxes in 2D Navier-Stokes equations?

I'm trying to solve 2D unsteady compressible Navier-Stokes equations with finite-difference or finite-volume method. Here is the system, it's 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 \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}\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}\end{array} \right)$$ where $$\tau$$ is viscous stress tensor: $$\tau_{xx} = \frac{4}{3}\mu \frac{\partial u}{\partial x} - \frac{2}{3}\mu \frac{\partial v}{\partial y} \\ \tau_{xy} = \mu \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right) \\ \tau_{yy} = \frac{4}{3}\mu \frac{\partial v}{\partial y} - \frac{2}{3}\mu \frac{\partial u}{\partial x}$$

Suppose that we discretize non-viscous flux terms (i.e. those parts of F and G vectors that do not contain derivatives) independently with some finite difference or finite-volume scheme, maybe with WENO reconstruction. Then viscous fluxes may be discretized and added as source terms.

With some of derivatives, it's pretty straightforward, for example:

$$\frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial x} \right) = \mu \frac{\partial^2 u}{\partial x^2} \approx \mu \frac{u_{i+1} - 2u_{i} + u_{i+1}}{\Delta x^2}$$ (second order approximation for uniform $$\mu$$).

But there are more complex derivatives in the viscous fluxes, such as: $$\frac{\partial}{\partial x} \left( \mu \frac{\partial v}{\partial y} \right), \; \frac{\partial}{\partial x} \left( \mu u \frac{\partial u}{\partial x} \right), \; \frac{\partial}{\partial x} \left( \mu u \frac{\partial v}{\partial y} \right), \; \frac{\partial}{\partial x} \left( \mu v \frac{\partial u}{\partial y} \right), \; \frac{\partial}{\partial x} \left( \mu v \frac{\partial v}{\partial x} \right)$$ and so on. And here are the questions:

1. What is the most simple way to discretize all these derivatives under finite-difference approach? I assume it would yield 2nd-order approximation.
2. How to use high-order central differences for them under FD approach?
3. What is the simplest way under finite-volume approach? I'm aware that we should use Gauss theorem in some form; there is some info in Blazek's book, but it's not detailed enough to be easily understandable for me. Are there more detailed books/papers on this?
4. How to apply high-order schemes in FV approach? It seems that method given in Blazek's book gives only 2nd-order approximation.
5. How to discretize these derivatives for non-uniform viscosity $$\mu$$ and heat conductivity $$k$$?
• Use a finite volume approach for this. Dont expand the second derivatives, that will lose conservation property. – cfdlab Nov 18 '20 at 12:23
• @cfdlab I understand that finite-volume approach is more "natural" in the sence of conservation. But how to compute those derivatives in finite-volume approach? – omican Nov 18 '20 at 18:50

Just open up the parentheses, e.g., $$\partial_{x} (\alpha \partial_{x} v) = (\partial_{x} \alpha) (\partial_x v) + \alpha \partial^2_x v$$, where $$\alpha=\mu$$ or $$\alpha=\mu u$$ etc., and apply your central differences:

$$(\partial_{x} \alpha) (\partial_x v) = (\alpha_{i+1}-\alpha_{i-1})(v_{i+1}-v_{i-1})/(4h^2)$$;

$$\alpha \partial^2_x v = \alpha_i (v_{i+1}+v_{i-1}-2v_i)/h^2$$

• Thanks! This way is the simplest, probably. But what about conservation? Is it "bad practice" to write and approximate source terms in non-conservative way? What problems may it introduce? – omican Nov 18 '20 at 11:01
• For enforcing conservation (as an algebraic identity) in the numerical scheme it is better to think in terms of the finite-volume paradigm (which you can mimic in a finite-difference scheme, at least on a uniform grid). For some classes of problems a conservative scheme is a great benefit, e.g., for shock front propagation. But for many problems a conservative scheme is not necessary. But as the grid resolution is increased, any numerical solution should converge to the exact solution of the differential equation which has those conservation properties. – Maxim Umansky Nov 18 '20 at 16:04

Use a finite volume method. Define $$\delta_x \phi(x,y) = \frac{\phi(x+\Delta x/2,y) - \phi(x-\Delta x/2,y)}{\Delta x}$$ $$a_x \phi(x,y) = \frac{\phi(x+\Delta x/2,y) + \phi(x-\Delta x/2,y)}{2}$$ etc.

For example, consider $$\tau_{xx}$$ which is required at $$(i+1/2,j)$$. $$(\tau_{xx})_{i+1/2,j} = \mu_{i+1/2,j} \left[ \frac{4}{3} \delta_x u_{i+1/2,j} - \frac{2}{3} a_x a_y \delta_y v_{i+1/2,j}\right]$$ This is just central differencing, and $$\delta_x u_{i+1/2,j} = \frac{u_{i+1,j} - u_{i,j}}{\Delta x}$$ $$a_x a_y \delta_y v_{i+1/2,j} = \frac{1}{2}\left[ \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y} + \frac{v_{i+1,j+1}-v_{i+1,j-1}}{2\Delta y} \right]$$ You can use this idea to write down approximations of all other terms.

• For a finite-volume method it is better to make velocity-like quantities face-centered rather than cell-centered, that will suppress some numerical artifacts – Maxim Umansky Nov 19 '20 at 16:04