Understanding Davis artificial viscosity

OlegSutyrin · August 6, 2017, 04:24

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

-th coordinate point and

-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},$

is the Jacobian matrix of vector

, $\rho(A(U))$ is it's spectral radius (for Euler's equations it equals $\text{max}(|u|+a, |u|-a)$ , where 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

,

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?

FMDenaro · August 6, 2017, 05:10

I think you should see the artificial viscosity as a term that disappears as dx->0. Therefore, multiply and divide by dx^2 and you get a second order derivative multiplied by dx^2.
You wrote dt/dx^2 but the dimensions are no longer congruent

OlegSutyrin · August 6, 2017, 06:00

Thanks for such a fast answer, FMDenaro!

Quote:

Originally Posted by FMDenaro

I think you should see the artificial viscosity as a term that disappears as dx->0. Therefore, multiply and divide by dx^2 and you get a second order derivative multiplied by dx^2.

So, you mean that artificial viscosity is not like physical viscosity (which effectively is a second derivative with coefficient that doesn't depend on grid), but more like physical viscosity multiplied by

?
But spurious oscillations do not disappear as $dx \to 0$ . If artificial viscosity disappears as $dx \to 0$ , how would it dampen the oscillations then?

Quote:

Originally Posted by FMDenaro

You wrote dt/dx^2 but the dimensions are no longer congruent

I can't say that I understand that completely. Could you elaborate?

FMDenaro · August 6, 2017, 06:10

The artificial viscosity has nothing to do with the physical viscosity. It is just an added terms that has the aim of dumping oscillations. But it has to be consistent with the original PDE equation as h->0.
I suggest a reading in the fundemental textbooks about the numerical solution of Euler equations.

OlegSutyrin · August 6, 2017, 06:20

Yeah, it seems that deep diving into the theory is inevitable :-)
Thanks again!

Martin Hegedus · August 7, 2017, 14:54

As a general statement, artificial viscosity is a consequence of numerical differentiation. It is not a physical term. It can be determined, in a hand waving way, from the wave equations. In a nutshell, central difference is a result of taking the limit when coming from the right and left.

In a sense, dF/dx is a function of lim (+ and -)->0 [dF/dx(+) + dF/dx(-)] + lim (+ and -)->0 [dF/dx(+) - dF/dx(-)]. If everything was perfect, and assuming no discontinuities, [dF/dx(+) - dF/dx(-)] would be zero. But numerically it is not. Thus the C*(U(i+1)-2*U(i)+U(i-1)) where C is a knob.

OlegSutyrin · August 12, 2017, 05:34

Quote:

Originally Posted by Martin Hegedus

If everything was perfect, and assuming no discontinuities, [dF/dx(+) - dF/dx(-)] would be zero. But numerically it is not.

Thanks for the explanation! As far as I know, ENO/WENO methods try to alleviate this very problem by choosing most smooth stencil over some region.

Quote:

Originally Posted by Martin Hegedus

Thus the C*(U(i+1)-2*U(i)+U(i-1)) where C is a knob.

"Knob" has too many meanings in English language, so I can't understand the meaning of this phrase :-)

FMDenaro · August 12, 2017, 07:21

I suggest having a reading to the texbook of Leveque (for example pag.71-72) to start with the concept of artificial viscosity explicitly added to the discretization in order to dump oscillations and numerical viscosity present in the local truncation error that is implicitly induced by the discretization.

OlegSutyrin · August 12, 2017, 07:33

Quote:

Originally Posted by FMDenaro

I suggest having a reading to the texbook of Leveque (for example pag.71-72) to start with the concept of artificial viscosity explicitly added to the discretization in order to dump oscillations and numerical viscosity present in the local truncation error that is implicitly induced by the discretization.

Could you provide the name of the book? There are several books by Randall J. LeVeque in the web, but I'm not sure which one do you mean (pages 71-72 of them contain something different).

FMDenaro · August 12, 2017, 07:37

http://www.cambridge.org/catalogue/c...=9780521009249

OlegSutyrin · August 12, 2017, 07:57

Quote:

Originally Posted by FMDenaro

http://www.cambridge.org/catalogue/c...=9780521009249

Ah, I missed this one because I was looking only on finite difference methods :-) Thanks!

August 6, 2017, 04:24	Understanding Davis artificial viscosity	#1
OlegSutyrin Member Oleg Sutyrin Join Date: Feb 2016 Location: Russia Posts: 41 Rep Power: 10	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 -th coordinate point and -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},$ is the Jacobian matrix of vector , $\rho(A(U))$ is it's spectral radius (for Euler's equations it equals $\text{max}(\|u\|+a, \|u\|-a)$ , where 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 , 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?

August 6, 2017, 05:10		#2
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	I think you should see the artificial viscosity as a term that disappears as dx->0. Therefore, multiply and divide by dx^2 and you get a second order derivative multiplied by dx^2. You wrote dt/dx^2 but the dimensions are no longer congruent OlegSutyrin likes this.

August 6, 2017, 06:10		#4
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	The artificial viscosity has nothing to do with the physical viscosity. It is just an added terms that has the aim of dumping oscillations. But it has to be consistent with the original PDE equation as h->0. I suggest a reading in the fundemental textbooks about the numerical solution of Euler equations. OlegSutyrin likes this.

August 7, 2017, 14:54		#6
Martin Hegedus Senior Member Martin Hegedus Join Date: Feb 2011 Posts: 500 Rep Power: 19	As a general statement, artificial viscosity is a consequence of numerical differentiation. It is not a physical term. It can be determined, in a hand waving way, from the wave equations. In a nutshell, central difference is a result of taking the limit when coming from the right and left. In a sense, dF/dx is a function of lim (+ and -)->0 [dF/dx(+) + dF/dx(-)] + lim (+ and -)->0 [dF/dx(+) - dF/dx(-)]. If everything was perfect, and assuming no discontinuities, [dF/dx(+) - dF/dx(-)] would be zero. But numerically it is not. Thus the C(U(i+1)-2U(i)+U(i-1)) where C is a knob. OlegSutyrin likes this.

August 12, 2017, 07:21		#8
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	I suggest having a reading to the texbook of Leveque (for example pag.71-72) to start with the concept of artificial viscosity explicitly added to the discretization in order to dump oscillations and numerical viscosity present in the local truncation error that is implicitly induced by the discretization. OlegSutyrin likes this.

August 6, 2017, 06:20		#5
OlegSutyrin Member Oleg Sutyrin Join Date: Feb 2016 Location: Russia Posts: 41 Rep Power: 10	Yeah, it seems that deep diving into the theory is inevitable :-) Thanks again!

August 12, 2017, 07:37		#10
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	http://www.cambridge.org/catalogue/c...=9780521009249 OlegSutyrin likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Problem with divergence	TDK	FLUENT	13	December 14, 2018 06:00
How can use MacCormack & Baldwin artificial viscosity in term of generalized nonortho	Nam Tran Ha	CFD Freelancers	0	December 21, 2016 03:31
Artificial viscosity	rads	CFX	2	July 31, 2006 07:02
artificial viscosity	rvndr	Main CFD Forum	1	March 2, 2004 10:24
Artificial Viscosity	Girish Bhandari	Main CFD Forum	5	January 13, 2001 13:22