# 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. 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? 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? 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. 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 Nov 19 '20 at 16:04