General difference approximations to the fluid dynamic equations require an artificial viscosity in order to converge to a steady state. This artificial viscosity serves two purposes. One is to suppress high frequency noise which is not damped by the central differences. The second purpose is to introduce an entropy-like condition so that shocks can be captured. These viscosities need a coefficient to measure the amount of viscosity to be added. In the standard scheme, a scalar coefficient is used based on the spectral radius of the Jacobian of the convective flux. However, this can add too much viscosity to the slower waves. Hence, it is suggested that a matrix viscosity be used. This gives an appropriate viscosity for each wave component. With this matrix valued coefficient, the central difference scheme becomes closer to upwind biased methods

Turkel, Eli

English

NASA Technical Reports Server

NASA Contractor Report 181712
ICASE REPORT NO. 88-53
ICASE
IMPROVING THE ACCURACY OF CENTRAL DIFFERENCE SCHEMES
Eli Turkel
Contract No. NASI-18107
September 1988
|I_ASJ-£Lq-161712| ]MPEOVIIqG I[1 ;ICCUi_ACi CF
I.I:I_'I_Ai I_IEFI_kE_C} 5CE£I_£S kital EepJct
|_/ISA) 15 _ SSCI 12A
I_89-IC577
Uncla_
G3/6_ 0 16 ¢,49 7,
INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND E_GINFJRING
NASA Langley Research Center, ltasptoa, Vir$inLs 23665
Operated by the Universities Space Research Association
Hat_d Aeronaut,sand
Since _S_a_on
Lan_ey Meeeereh Cemw
I.taml:_n,VJrgJnJa 23665
https://ntrs.nasa.gov/search.jsp?R=19890001206 2020-03-20T05:58:42+00:00Z

IMPROVING THE ACCURACY OF CENTRAL DIFFERENCE SCHEMES
Eli Turkel 1
School of Mathematica_ Sciences
Sackler Faculty of Exact Sciences
Tel-Aviv Univer,,ity
ABSTRACT
Central difference approximations to the fluid dynamic equations require an artificial
viscosity in order to converge to a steady state. This artificial viscosity serves two purposes.
One is to suppress high frequency noise which is not ,iamped by the central differences. The
second purpose is to introduce an entropy-like con, tition so that shocks can be captured.
These viscosities need a coefficient to measure the amount of viscosity to be added. In the
standard scheme, a scalar coefficient is used based o a the spectral radius of the Jacobian of
the convective flux. However, this can add too mucll viscosity to the slower waves. Hence,
we suggest using a matrix viscosity. This gives an appropriate viscosity for each wave
component. With this matrix valued coefficient, 1he central difference scheme becomes
closer to upwind biased methods.
IThis research was supported by the National Aeronautics and Space Administration under NASA Con-
tract No. NAS1-18107 while the author was in residence at the Institute for Computer Applications in
Science and Engineering (ICASE}, NASA Langley Research Center, Hampton, VA 23665.

I. Introduci,ion
In recent years, central difference schemes h_ve been used with much success to solve
transonic flow problems about aerodynamic shapes. These schemes are second order ac-
curate for sufficiently smooth meshes and have an added artificial viscosity to stabilize
the scheme and reach a steady state. This artificial viscosity is usually a blend of two
terms. One is a fourth difference that stabiliz,;s the even-odd modes that appear with
central differences and constants coefficients. Without this viscosity, one cannot reduce
the residual beyond some level because of a remaining high frequency mode. The second
viscosity term is a nonlinear second difference tl:at limits oscillations in the neighborhood
of shocks. A nonlinear shock detector preserves the second order accuracy of the scheme
in smooth regions.
An advantage of the artificial viscosity approach is that it allows the user control
over the amount of dissipation. Nevertheless, o_m sometimes finds that there is too much
dissipation in the numerical solution. Changing the global constants that appear in the
formulas is not sufficient to construct an artifici_d viscosity that is appropriate in both the
shocked and smooth regions of the flow. For some problems, we need to severely limit
the viscosity in some smooth regions, e.g., neaJ the trailing edge, while still maintaining
stability near the shocks. Hence, although the standard artificial viscosity works well in
most cases, it is not sufficiently flexible to handle more delicate problems.
In order to improve the existing artificial vi:_cosity, we shall make use of ideas used in
the construction of upwind schemes. In particular, we shall replace the scalar coefficient in
the artificial viscosity by a matrix. To prevent difficulties near stagnation points, a cutoff
is introduced that depends on the spectral radius of the matrix. By varying the cutoff,
one can obtain an appropriate average betweerL the original scalar viscosity and the new
matrix viscosity. Since the matrix viscosity reduces the amount of smoothing on the slower
waves, it will improve the total accuracy of the scheme.
as
where
II. Finite Volume Formulation and Artificial Viscosity
The Euler equations for an inviscid compres.':_ible flow can be written in divergence form
Of Og Ohi) Q -I- -I- -+--- = 0 (1)
a-7
Q = (p, pu, pw, E) t (2a)
and for an ideal gas
f = (p=,p=' + ,, p,v, p,=, (E + ,)_)'
g = (pv.puv,pv' + p,pvw,(E + p)v)t
h = (pw,puw,pvw,pw2+ p,(E + p)w)t
p = CA/-- 1)[E - pCu 2 + v 2 + w2)/2].
We can also write (1) in the form
(2b)
(2_)
(2d)
(2e)
OQ
O--T + div(F) = O. (lb)
We integrate (1) over a three dimensional cell and consider Q_j,k as an approximation
to the average of Q over the cell. Hence,
cOQi,i,_ f f f divFdV
cgt f f f dE
or using the divergence theorem,
=0
+f f  ..as o.
Hence, the time change of the average Q is governed by the fluxes entering and leaving the
cell.
This finite volume approach leads to a pure central difference method for Cartesian
grids. Though this scheme is stable for constant coefficient hyperbolic equations it is
subject to instabilities that will prevent the convergence to a steady state. To enhance
this convergence a fourth difference viscosity is added to the scheme. The fourth difference
causes oscillations in the neighborhood of shocks. Hence, a shock detector is constructed
and near the shocks the fourth difference is turned off and only the nonlinear second
difference is operative. The total artificial viscosity, V, is the sum of such second and
fourth differences in each coordinate direction.
V,o, = , . - , • + W. , - W,.i-_.k
_,$,k -- _,$,k i,j+ _,k
+V[ij,k+ _- - V--_i,j,,- _- (4)
Hence it is sufficient to describe these terms in the _ direction. Since we only take differ-
ences at neighboring points the artificial viscosity is always in conservation form.
The first difference is defined as
Di+, 'k = Qi+ld,k - Qij, k (5a)
and the second _ difference is defined as
Ei,i,k = Di+, ._ -Di_ (5b)
,_, ½,j,k"
We then form the second and fourth differences. In particular the fourth difference is
formed as a second difference of a second differe:lce with positive weights [3,8]. Hence,
V--_+½,i,k -(_) D , • - (e !4.:- . E,+I "k- e_),kE,,i,k).
= "i+ ½,i,k i+ _,_,k x t-r l,i,lt ,_,
(6)
Let,
Pi+l,i,k - 2pi i,k + Pi-l,i,k (7a)
Vi,j, k
Pi+l,i,k + 2pi i,I, + Pi-l,i,k
t_i,j,_ is used to detect the location of shocks. W_Jen vi,i,k is large then the fourth difference
is reduced. Let,
a,+ ½,i,k = K(*) max(v,_ ,,i,k, _,,i,k, v,+ ,,i,k, v,+ 2,i,k). (8)
We also multiply a by a function of the Mach n umber to reduce a near the surface. Let,
aT a_ a_ where F,G,H are the fluxes in the coordinate system (_,77,_).A = _-_,B = _-q, C = _,
Let A be a measure of the fluxes. The original code chose ,k as
_ = _,7 = _, = p(A) + p(B) + p(C) (9a)
where p is the the spectral radius of the matrix. For problems with a highly stretched
mesh it was found [2,3,8] that for increased accuracy one should choose
A¢= p(A), = , = (9b)
K(2), K (4) are constants that determine the hvel of the second and fourth differences.
These constants are given as input to the code. Then
t(2) =),. ,. a. i • (10a)
i+ ½,./,k ,+ _-,:,k t+ _,_,k
e_),k = _',,i,k max(0, K c4) - a_,i,j, ). (10b)
In order to imitate the upwind type [7] algcrithms we now replace the scalars in (9b)
by matrices. Hence,
he = IAI, _, = II_l ' _, = iel (9c)
where IAI = TIA¢I T-_ when A =TA¢T -_ and h,, is a diagonal matrix with the eigenvalues
of A as its entries. This definition of $ can lead to difficulties when an eigenvalue is near
zero. Hence, we modified A_ to be
A_ = diag(max(ai, qp(A)) ai = e.v. of A (11)
3
where q is a specified constant. When q -- 1, then A_ = p(A) • I and so (9b) is recovered.
When q = 0 then no modification of A_ is done. In general, we found that q = 0.2 gives
good results.
We point out that the use of (9c) does not allow for a constant enthalpy solution and
so enthalpy damping cannot be used [5].
III. Results
We consider the central difference code with Runge-Kutta time stepping [3,5]. As de-
scribed above we use a matrix valued artificial viscosity which approximates TVD type
schemes [6,7,9]. The fourth difference viscosity is still needed to allow the multigrid accel-
eration to quickly reach a steady state. We consider inviscid flow about a NACA0012 with
M_ = 0.8, a -- 1.25 °. A 192 x 32 C mesh is used with 128 points on the airfoil. In [1] it is
pointed out that the standard code smears the weak shock on the lower surface. In Figure
1, we plot the result for the standard scheme, but without enthalpy damping. In Figure
2, we show the same case but using the matrix viscosity. The convergence rate is slowed
down since the fourth order viscosity is not as strong but the shock on the lower surface is
sharper. There is an overshoot on the shock on the upper surface. This is due to fact that
the cutoff (10b) is not sufficiently sharp. One way to improve this is to replace (7) by
[Pi+ld,k -- _i,j,k[ - IPi,l,k -- Pi-l,j,k[
v,,j,k -= ]P,+l,i,k - p,,j,k[ + IP¢,¢,k - Pi-l,i,_[ + c (Tb)
so that vid, k is one at discontinuities. One can also use the matrix coefficient only for
the second difference but use a scalar coefficient for the fourth difference viscosity. This
accelerates the convergence to a steady state but smears the shocks. The results presented
used a four step Runge-Kutta algorithm with the artificial viscosity frozen after the first
stage using a matrix viscosity rather than a scalar viscosity adds about 60% to the total
CPU time.
References
[1] Anderson, W., Kyle, Thomas, J. L., van Leer, B., Comparison of Finite Volume Flux
Vector Splittings for the Euler Equations, AIAA J., Vol. 24, pp. 1453-1460 (1986).
[2] Caughey, D. A., Turkel, E., Effects of Numerical Dissipation on Finite-Volume Solu-
tions of Compressible Flow Problems, AIAA-88-0621 (1988).
[3] Chima, R. B., Turkel, E., Schaffer, S., Comparison of Three Explicit Multigrid Meth-
ods for the Euler and Navier-Stokes Equations, AIAA-87-0602 (1987).
[4] Harten, A., The Artificial Compression Method for Computation of Shocks and Con-
tact Discontinuities III. Self-Adjusting Hybrid Schemes, Math. Comput., Vol. 32, pp.
363-389 (1978).
[5] Jameson, A., Schmidt, W., Turkel, E., Numerical Solutions of the Euler Equations by
Finite Volume Methods using Runge-Kutt._ Time-Stepping Schemes, AIAA-81-1259
(1981).
[6] Osher, S., Tadmor, E., On the Convergence of Difference Approximations to Scalar
Conservation Laws, Math. Comput., Vol. 50, pp. 19-51 (1988).
[7] Roe, P. L., Approximate Riemann Solvers Parametric Vectors and Difference Schemes,
J. Comput. Phys., Vol. 43, pp. 357-372 (1981).
[8] Swanson, R. C., Turkel, E., Artificial Dissipation and Central Difference Schemes for
the Euler and Navier-Stokes Equations, A]AA 8th CFD Conference, AIAA-87-1107-
CP (1987).
[9] Tadmor, E., Convenient Total Variation Diminishing Conditions for Nonlinear Differ-
ences Schemes, SIAM J. Numer. Anal., to appear.
[10] Turkel, E., Acceleration to a Steady State for the Euler Equations, Numerical Methods
for the Euler equations of Fluid Dynamics_ pp. 281-311, F. Angrand, et al. (editors),
SIAM, Philadelphia (1985).
[11] Turkel, E., van Leer, B., Flux Vector Splitting and Runge-Kutta Methods for the
Euler Equations, 9th Inter. Conf. Numer. :Methods Fluid Dynamics, Springer-Verlag
Lecture Notes Physics, Vol. 208, pp. 566-570 (1984).
[12] van Leer, B., Towards the Ultimate Conservative Difference Scheme, II. Monotonicity
and Conservation Combined in a Second-Order Scheme, J. Comput. Phys., Vol. 14,
pp. 361-370 (1974).
[13] Yee, H. C., Upwind and Symmetric Shock-Capturing Schemes, NASA TM 89464
(1987).
G-O
(1)
c_
E7_
O
._J
-2
-4
-6
-8
-10
-12 I
0 200
I I
4-00 600 800
Iterations
J I
1000
c-
.__
c)
.m
N---
_D
O
C.)
o__
d
.35 _'--
.30
.25
.20
.15 -
.10 -
.05 -
0 I
0 200
I I I I
400 600 800
Iterations
1 00(
NACA 0012 FLO_:x3
192 X 33 GRID _4 m
-2.0 -
-1.5 -
-1.0
Cp
0
.5
1.0-
1.5
Figure 1: Scalar viscosity
.800 ALPHA -- 1.250 CD. CL - .02.._,_,_,_,_,_,_,_,_,.5
J
_I I I I I I
0 .2 .4 .6 .8 1.0
x/c
. 65
(9
-(3
(1)
rY
ET_
0
._J
2
0
-2
-4
-6
-8_
-10 --
-12
0
I I I I
200 400 600 800 10',30
Iterations
0
©
©
(D
°__
.._J
.40
.35
,30
.25
.20 -
.1,,.%-
,10 -
,05
0
I
200
I t I I
400 600 soo '1o0(
Iterations
NACA 0012. FLOS,.3-VEPS
192 X 5,3 GRID M -- .800 ALPHa,,.- 1.250
-2.0
-1.5
-1.0
Cp
--.5
0
.5
1.0
1.5
Figure 2: Matrix viscosity
m
_1 I
.2
CD. CL m .023.345 .365
I I I
.6 .8 1.0
Appendix
We present the matrix viscosity in explicit form for three dimensions in general curvi-
linear form. Let
be a change of variables from physical space (x, y, z) to a computational space (_, r/, g) such
that the curvilinear mesh is mapped to a cube. Then (lb) can be rewritten as
where
aQ OF OG OH
-_- + _ + _-_+ -_- =0
G = rl,f + yyg+ rl,h
H = gzf+gvg+gzh.
We next express (A1) in quasilinear form
_-+A +B +C =0
Let
OF OG OH
A- B- c-
aQ aQ aQ"
q = _::u H- _vv + _,w
(A1)
(A2)
(A3)
8
be the three contravariant velocities.
V _= _2_-(u2 + v 2 + w 2) then
0
al V 2 - uq
a2V _ - vq
a3V 2 - wq
-q(h- V _)
al
q- (7- 2)alu
aav - (7 - 1)a2u
alw- (7- 1)a3u
alh- (7- 1)qh
rh:U + rlvv + rl,w
Also define h = E+p as the total enthalpy and letP
a2 as 0
a,w - (7- 1)asv q- (7- 2)azw (7- 1)as
,,_h- (7- 1)qv ,,_h- (,_- 1)qw 7q
(A4
8
al=_,, a2=_iv, as=_,.
For B we get a similar matrix with al = r/z, a2 = r/y, as = r/, and q replaced by r while
forGwehaveal=qz, a2=_'v, a3= G and q replaced by s"
Hence, we can find the absolute value of z_ll three matrices in the same way. Let us
assume that A has eigenvalues
)kl
k2 0
),s
0 ks
(A5)
where Xl = q + x/a_ + a_ + a_c,
redefined by (11) to prevent )_y from approaching zero.
In order to find [A I it is easier to use a tw,:> step procedure. Let
k, = q - V/(t_I_+ a] + a_c and k3 = q. In practice ky, are
AI = T:_AT{ I (A6)
be a symmetric matrix. Since we can symmelrize A, B, G simultaneously the same T1 will
work for all three matrices.
V 2
C
--it
--(q- 1)ulc -('7- 1)vie --('7-- 1)w/c ('7- 1)1c "_
1 0 0 0
0 1 0 0
0 0 1 0
-('7-- 1)u -(7- 1)v -('7- 1)w "7- 1
1 0
C
It 1
C
v 0
TI-I = c
w 0
C
h
C
0 0
C 2
It
0 0
C 2
1)
1 0
C 2
I13
0 1
C 2
1) tO
It2 + v 2 + w 2
2c 2
(A7)
9
xs
As expected T1, T[ 1 do not depend on the matrice al, a2, a - 3.
q ale a2c asc 0
ale q 0 0 0
A1 - a2c 0 q 0 0 (A8)
aac 0 0 q 0
0 0 0 0 q
Since A1 (and B1, C1) is symmetric we can diagonalize it with a unitary matrix T2. T2 will
depend explicitly on ax,a2,a3 and so is different for A1,B1, el. Let A = Ca_ + a22+ a23 -# 0,
then
A al a2 as 0
and
1
T2-
-_ al a2 as 0
0 xl x2 Xs 0
Yl Y2 Y, 0
o o o v_._
2-1 ___ T2 $.
The (xj, yj) are numbers that satisfy the following equations
(A9)
Zl2+ z2+ z2=y[ + v_+ vl : 2_2
alxl + a2x2 + asX3 = alYl + a2Y2 + asys = 0
xlyl + x2y2 + x3Ys = 0
XlYl + x2y2+ xsvs = 0
xlx2 + YlY2 = -2ala2
xlxs + YlYs = -2ala3
(A10)
x2xs + Y2Ys = -2a2as.
It is difficult to give explicit formula for the xi, y_- in all cases since some of the a i may be
zero as long as a_ + al + as2 ¢ 0. However, the final formula does not depend on explicitly
knowing the a i. Given T2 we find that A (A5) is given by
A. = T2AIT_-I= T2TIA(T2T1) -1. (All)
10
We now reverse the process and define
IAI =
0
0
1:_31 / (A12)
where ,_j can be modified eigenvalues of A. Tlen
IAI = (T_TI)-_IAICT_T,). (A13)
Let
AI + A2
U1 -- _ir2 --
2 2
and define the row vectors
ll = (7 - 1)( u2 + v2 + w22 ,--u,--v,--w,i)
(A14)
(A15)
12 = (-q, al, as, as, 0).
We then have the matrices
I = 5 x 5 identity member
g, = (e,, ul,, re1, Wll, hll)'
Zs = (_2, u£2 -I- alt.l, V_s -I- a_l, W_ -Jr- asel,h_ + qel) t
Z3 = (o, axg2,a2e2,a3_,qe2)'
and finally,
Because of the simple nature of the matrices _i it is easy to multiply IAI times a vector.
D e fine,
rx = (1, u,o,w,h) t
r2 = (0, ax,-2, as, q)'.
Let ( , ) denote the standard inner product, then
(_l,rl) : Cs
(_l,r2) : 0
11
(12, rl) = 0
(_2,r2)= _2
and if x is any column vector, then
[( ) I [°,tAIz= Asz+ Ol- _,_ (ll,Z) + ,,2 z) rl+ +_2 _Cl2, _(ll,Z)
- _) ]("la 2 (12,z) r2. (A18)
i2

I_l*Oq_l f_e_C_nau_ s anc_
1. Report No.
NASA CR- 181712
ICASE Report No. 88-53
4. Title and Subtitle
IMPROVING THE ACCURACY OF CENTRAL DIFFERENCE
SCHEMES
7. Author(s)
Eli Turke I
Report Documentation Page
2. Government Accession No. 3. Recipient's Catalog No.
5, Report Date
September 1988
6. Performing Organization Code
8. Performing Organization Report No.
88-53
10. Work Unit No.
505-90-21-01
11. Contract or Grant No.
NASI-18107
13. Type of Report and Period Covered
Contractor Report
9. Performing Organization Name and Address
Institute for Computer Applications in Science
and Engineering
Mail Stop 132C, NASA Langley Research Center
Hampton. VA 23665-52P5
12. Sponsoring Agency Name and Address
National Aeronautics and Space Administration
Langley Research Center
Hampton, VA 23665-5225
15. Supplementary Notes
14. Sponsoring ,_,gency Code
Langley Technical Monitor:
Richard W. Barnwell
Final Report
Submitted to the 11th International
Conference on Numerical Methods
and Fluid Dynamics
16. Abstract
Central difference approximations to the fluid dynamic equations require an
artificial viscosity in order to converge to a steady state. This artificial vis-
cosity serves two purposes. One is to suppress high frequency noise which is not
damped by the central differences. The second purpose is to introduce an entropy-
like condition so that shocks can be captured. These viscosities need a coef-
ficient to measure the amount of viscosity to be added. In the standard scheme, a
scalar coefficient is used based on the spectral radius of the Jacobian of the
convective flux. However, this can add too much viscosity to the slower waves.
Hence, we suggest using a matrix viscosity. This gives an appropriate viscosity
for each wave component. With this matrix valued coefficient, the central dif-
ference scheme becomes closer to upwind biased methods.
17. Key Words (Suggested by Author(s))
artificial viscosity, Runge-Kutta,
Euler equations
19. SecuriW Cla_if. (of this report)
Unclassified
NASA FORM i_l_ OCT 86
18. Distribution Statement
64 - Numerical Analysis
Unclassified - unlimited
_. SecuriW Cla_if. (of this pa_) 21. No. of pa_s
Unclassified 14
22. Price
A02
NASA-Langley, 1988


Improving the accuracy of central difference schemes

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server