NASA/TM-1998-208738

Comments on the Diffusive Behavior of

Two Upwind Schemes

William A. Wood and William L. Kleb

Langley Research Center, Hampton, Virginia

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23681-2199

October 1998

https://ntrs.nasa.gov/search.jsp?R=19990009124 2018-05-11T16:54:48+00:00Z
Available from:

NASA Center for AeroSpace Information (CASI)

7121 Standard Drive

Hanover, MD 21076 1320

(301) 621 0390

National 2echnical Information Service (NTIS)

5285 Port Royal Road

Springfield, VA 22161 2171

(703) 605-3000
Abstract

The diffusive characteristics of two upwind scheines, multi-dimensional fluctuation splitting and locally one-dimensional finite volume, are compared for scalar advection-diffusion

problems. Algorithms for the two schemes are developed for node-based data representation oil median-dual meshes associated with unstructured triangulations in two spatial

dimensions. Four model equations are considered: linear advection, non-linear advection,

diffusion, and advection-diffusion. Modular coding is employed to isolate the effects of the

two approaches for upwind flux evaluation, allowing for head-to-head accuracy and efficiency

comparisons. Both the stability of compressive limiters and the amount of artificial diffusion

generated by the schemes is found to be grid-orientation dependent, with the fluctuation

splitting scheme producing less artificial diffusion than the finite volume scheme. Convergence rates are compared for the combined advection-diffusion problem, with a speedup of

2.5 seen for fluctuation splitting versus finite volume when solved on the same mesh. However, accurate solutions to problems with small diffusion coefficients can be achieved on

coarser meshes using fluctuation splitting rather than finite volume, so that when comparing

convergence rates to reach a given accuracy, fluctuation splitting shows a speedup of 29 over

finite volmne.

Nomenclature

A, B

c

P

h

j-1

M

fi,

n

Q

r_

&

t

X,y

F

7

0

Flux Jacobians

Nodal update coefficients

Convective flux vector, /? = F(x, y, u)

Mesh spacing

Cartesian unit vectors

Inverse Jacobian of the coordinate transformation

Edge length

Averaging function for limiters

Outward unit normal to control cell

Length-scaled inward normal

Fluctuation ratio

Distance vector from node to face

Median-dual area about node i

Time

Dependent variable

Cartesian coordinates

Curvilinear advection speeds

Boundary of control cell

Limiter bound

Small parameter for van Albada limiter

Finite element linear shape function

Advection velocity vector
//

T

m

Q

V

Diffusion coefficient

Curvilinear coordinates

Timestep

Artificial dissipation--fluctuation splitting

Artificial dissipation--finite volume

Element advective fluctuation

Fluctuation components

Limited fluctuations

Viscous fluctuation

Limiter function

Area of control cell

Gradient operator

Introduction

Upwind discretizations for advection equations typically introduce artificial numerical

dissipation into the solution. When combined advection-diffusion problems are considered,

this dissipation introduced in the discretization of the advection terms should be less than

the true physical diffusion. To this end the diffusive characteristics of upwind schemes are

investigated and their performance in resolving solutions to advection-diffusion prot)lems

with small diffusion coefficients is analyzed.

Two node-based, median-dual methods for modeling convective fluxes are considered.

The first is a traditional locally one-dimensional approximate Riemann solver finite volume

(FV) scheme. 1 Locally one-dimensional schemes applied on multidimensional domains are

known to introduce excess dissipation when discontinuities are not aligned with the mesh. 2

The second method is the NNL 3 fluctuation splitting (FS) scheme, also referred to in

the literature as a residual distribution scheme. FS has a more-compact stencil than FV for

second-order formulations and exhibits "zero cross-diSusion"t in a grid-aligned condition.

Both of these attributes should lead to less introduced dissipation as compared with FV.

The sensitivity of FS to grid orientation and res_lting production of cross-diffusion is

investigated in the present report. The use of compressive limiter functions is also tested

with both algorithms. Local timesteps based on posi|ivity arguments are implemented for

both first- and second-order discretizations of the imp!icit matrix.

Formulation of FS schemes for diffusion problems is a recent research area. 4,5 The present

study seeks to quantit_ the relative merits of using _ low-diffusion advection operator to

resolve advection-diffusion problems with small diffu, ion coefficients. Lessons learned on

these problems will guide the development of computer :odes for solving compressible viscous

fluid dynamic problems. A similar approach for central difference schemes with explicit

numerical dissipation has recently been taken by Efraimsson. 6

t "Zero cross diffusion" refers to the practice of adding artificial diffusion terms in the streamwise direction

only, as opposed to adding artificial dissipation in both the stre:_mwise and cross-stream directions.
Governing Equations

Tile non-linear advection-diffusion equation,

u, + V. f = V. (uVu) (1)

is cast as a hypert)olic conservation law, to which steady-state solutions are sought.

Finite Volume

In FV form, using the divergence theorem Eqn. 1 t)ecomes,

where Q is the median dual at)out node i and F is the t)oundary of _. Using mass lumping

to the nodes, similar to an explicit finite element treatment, 8 the temporal evolution is

evaluated on a time-invariant mesh as,

u, dQ = Si-_ -+ -- - 'u_) (3)

, T

The discretization of the convective flux, k_, is performed using Barth's implementation 1

of the upwind, locally one-dimensional, approximate Riemann solver of Roe. 9

f aC_?,g

where the artificial dissipation provides the upwinding,

= _ I.i,*_ + bGl(uo, . - _,.,) (5)

with f_,= i_i + h.vj. Out and in refer to states on the outside and inside of _ at the face. A

and B are the flux ,]acobians,

OF (1) 0F (2)

A--0_ ' B- 0u (6)

and (.4,/3) represent their conservative linearizations at the cell face. 9

Piecewise linear reconstruction from the nodal unknowns to the cell faces as,

u:_c_ = u,: + _',V'. • _" (7)

provides second-order spatial accuracy in smoothly-varying regions of the solution. Mediandual gradients of the dependent variable, X7u, are obtained from the unweighted least squares

procedure outlined by' Barth. Following Bruner and Walters, 14 the limiter is supplied an

argument equal to half the argument Barth uses, namely,

'_), = t/) (2 (8)

?lmin / raax
where u ''i'O'_'_ is tile minimum (resp. maximum) of ui and all distance-one neighbors. The

most restrietiv(, limiting from choosing the minimum or maximum is used.

In casting the limiter argument in this form, Bvuner equates the Barth limiter with

Superl)ee, for a limiter argument less than or equal to one. The present authors incorrectly

identified the Barth linfiting with the non-symmetri( Chakravarthy and Osher la limiter in

Ref. 7. The Barth limiting is non-symmetric, but takes the form,

@ P_<O

_/' = if 0 _. _ <

q

q

(9)

for tile limiter cast as Eqn. 8.

Two methods for evaluating the diffusion term are incorporated into FV. Tile more

compact of the two, the finite eleInent discretization, is discussed in the following section.

The less-compact diffusion formula is obtained by discretizing the last term of Eqn. 2, in a

lnanner similar to Eqn. 4,

_i P (Vui,, + VUo,_t). fi AF (10)

face8

The diffusion coefficient is averaged over the length of the face. The gradients from Eqn. 7

are not limited before averaging at the control-volume faces in Eqn. 10, as suggested by

Anderson and Bonhaus. t°

Fluctuation Splitting

Tile NNL FS scheme is presented as a slight re-interpretation of tile work of Sidilkover and

Roe. a Tile current interpretation is as a volume integral over triangular elements, without

recourse to tile divergence theorem. Tile discretized equations, however, are identical.

Integrating Eqn. 1 over an eleinent, where f_ is nox tile area of tile triangular element,

For linear variation of the dependent variable over the element, tile temporal evolution is,

(711, + _- ) (12)

f_

ut d_ = l_'/it = _ u2, u3,

where u l, 'u2, and Ua correspond to the three nodes defining element fL

Defining local curvilinear coordinates, { and r/, parallel to sides 1-2 and 2--3, respectively

(Fig. 1), the divergence of the convective flux can be written,

V-F=F_ ')+F_ 2)- I @'2"*_,,-h,'Z) (13)

J-I

Defining the scaled inward normal, n = -hh,, where h is a mesh edge length, the divergence

(Eqn. 13) be('omes,

v.P= 2--fi l (-h12n2"r_+h23nl'lW,_) (14)
If F is linear or quadratic in u, then for a linear variation of u over the element,

V • fi df_ = a'A21'u +/3Aa2U

where tile difference operator is defined A2_u = we - Ul and the advection speeds are,

1 - /_), fl=l

a: = -_(n2_A + n2_ _(n,x,4 + nl./))

.4 and/_ are now the conservatiw, linearizations over the triangular element. 11

The advective fluctuation can be defined,

,=-£V-Fd_

The fluctuation can be split,

where,

(15)

(16)

(17)

with,

_=0{+_" (18)

0 _ = --_'A21tt, 0 n = --3/ka2U

Following SidilkoveP 2 the fluctuation is limited to achieve a second-order scheme,

O_" = O_ + c)'{J(Q) = O_ (1 _) )

&,l' = 0,1 __ O"_/,(Q) = O'J (1 - _J_(Q))

19)

(20)

(21)

cO_ (22)

Upwinding is achieved through the introduction of the artificial dissipation terms,

¢_ = sign(a)0 C, $" = sign(3)0 '/" (23)

Combining Eqn. 12 with a distribution scheme for Eqn. 17 and sumlning over all elements,

the contributions to nodal time derivatives can be written in the form,

1

Sl_/I, +--- _(0_* - 6 _) --}-COG

s2_,_,+-- _(0_"+ ) + 5( - 6") + coe

$3_/3, _ 2(Or/" _+_6r/) + COE (24)
oi"in a morecorot)actform,

1 [i(3 - i)(O c + (-1)i_ ¢) + (-4 + 5i -/2)(0 "" - (-1)i0'_)] + COE (25)

where COE stands for contributions from other elements containing these nodes.

A finite element treatment, sinlilar to Tomaich, 4 is employed to obtain the diffusive

fluctuation,

,,, v. (26)

Assuming piecewise-linear data and an element-averaged diffusion coefficient leads to a diffusive fluctuation of zero for the triangular element. Introducing the linear nodal shape

flmctions Ji, such that _-_i_10i = 1, the elemental diffusive fluctuation can be expressed

i=l 0_ = O, where

0v, = ./_ 0_V. (OVui dl_ (27)

Integrating by parts,

Ov_=frOiOVu'fidF-f#Vu'VvOidQ (28)

The boundary integral in Eqn. 28 will cancel on summing contributions for interior nodes.

The remaining volume integral can be evaluated anal) tically,

3

- #

0,,_ - uVu'2 ni+l -- 4--_ Z ujnj+l " ni+l (29)

j=l

Distributing this diffusive fluctuation to the nodes and keeping only tile larger of the I)hysical

or artificial dissipation leads to the update formula,

SI _tlt _ _ -[- max - ,0_ + COE

$2u2, _ 0C+02 _" + max ((4_-_))2 ,0_,2 + COE

_3'_3t_--07'*--2 + Inax (_ , Ov3)+_'OE

(30)

Boundary Condb;ions

Explicit Dirichlet inflow boundary conditions are ,_'mployed. Advective outfow boundaries are treated for free convection through tile boundary nodes, allowing boundary nodes

to be handled in the same manner as interior nodes. For the diffusion terms a Neumann

outflow I)oundary is applied with zero gradient, achieved by setting the boundary integral

in Eqn. 28 to zero.

6
Limiter functions

Minmod, van Albada, 15 Superbee, and 7 j6 symmetric limiters are utilized for FV (Eqn. 7)

and FS (Eqns. 20 and 21) in tile form of symmetric averaging functions related to the limiter

as,

The van Albada averaging flmction is,

M = (pq+ e2)(P+ q)

p2 + q2 + 2_2

where the small parameter c varies like c2 _ Ax 3, and serves to reduce tile limiting in smooth

regions.

The averaging function for the ?' limiter, of which tile Minmod (7 = 1) and Superbee

(_ = 2) are special cases, is,

I 0 pq<_O

7p 7lpl_<lq[

M(p,q) = q if IPl-<lql-< 71Pl

p Iql-<lpl-< 71q[

7q 7]ql<-lpl

(31)

Timestep

Both schemes are formulated either as Gauss-Seidel time-relaxation or forward Euler

time-evolution algorithms.

The nodal ul)dates for the discrete system can be formed as a sum of contributions from

all nodes.

i =

j jCi

For positivity 17 each of the coefficients in Eqn. 32 must be non-negative.

Advective Timestep restriction

In the FV context the nodal update (Eqn. 32) can be rearranged into the form of Eqn. 3,

Si(u[+_ - u[) = (ci- l)ui+ -- E cj'uj (33)

T T j#i

For the upwind, edge-based algorithm considered here, each @cj will be related to a positivedefinite coefficient equal to zero for outflowing faces and related to the wavespeed for inflowing faces, yielding the restriction r _> 0 on the timestep. The remaining term can be
expressed,

S'(c,- 1)= - (34)

T

k abou_ i

where the Ck coefficients are also positive-definite, either zero for inflowing faces or related

to the wavespeed for outflowing faces. Rearranging and imposing the positivity constraint,

Ci :> O, yields tile timestep restriction,

1 - Ck = Ci -_'"0 (35)

k about i

&

r < _k .boo_ick (36)

For FS, tile nodal updates are assembled from Eqi_. 24 as,

sa(_'+* - u_)= _ c,(, j - _,) (37)

7" i

j¢i

In this case the cj coefficients are formed as contributions from the fluctuations in the

triangles to both tile left and the right of mesh edge _. The positivity restriction on r is

found to have a similar form as for finite volume (Eqn. 36),

&

r _<_2j+_cj (as)

Local time-stepping based on positivity is shown t_ yield stable, yet non-converging, solutions in some second-order cases (see Results section). Robust convergence is obtained by

using the first-order o's in Eqns. 36 and 38, even for se(ond-order-accurate spatial discretizations. This is equivalent to the common practice of using a first-order Jacobian discretization

in an time-implicit scheme.

Diffusive Timestep Restriction

Unfortunately, tile finite element formulation for t_e diffusive terms (Eqn. 29) cannot be

guaranteed to preserve local positivity on obtuse triangles (see Barth_). Considering only

the contributions from the current node, the coefficient for the diffusion term can be written,

=u,,(1& _, i-_] (39)

The appropriate edge length is tile side of tile elem, mt opposite the current node. The

resulting timestep restriction is,

&

r < (40)

-- ET _,g2

In a similar manner the timestep restriction from Eqn 10 is,

Si

7-<

3uAF 2

-- ET 4f_

(41)
Results

Linear Advection

The linear advection equation is obtained froin Eqn. 1 by setting v = 0 and /_ = _tt,

yielding,

'u, + V. (Au) = 0 (42)

A divergence-less advection velocity is considered, such that V-_ = 0. Equation 42 can then

be written,

'at + ._- Vu = 0 (43)

Uniform Advection

Uniform advection of the Heavyside function at -45 degrees, ,_ = (1,-1), on a cutcartesian mesh is shown for first-order FS, second-order FS, and second-order FV in Figs. 2 4,

respectively. The mesh is shown as the dashed 1)ackground, and equally-spaced contours vary

on [0,1], the minimum and maximum solution values. Tile spread of the contour lines with

spatial evolution is indicative of the amount of dissipation introduced into the solution by

the discretization of the convective terms.

Second-order FS is seen to be greatly superior to first-order, as expected, reproducing the

exact sohltion in this case with no introduced dissipation. Also, FS represents a significant

reduction in numerical diffusion versus the corresponding FV scheme, with both results

employing the Minmod limiter.

However, the "zero cross-diffusion" results of Fig. 3 with FS are misleading. In Fig. 5

the advection velocity has been rotated counter clockwise by 90 degrees on the same grid.

Clearly, the artificial dissipation introduced by the FS scheme has been increased.

The corresponding FV solution is shown in Fig. 6. While the change in contour spreading

for the FV scheme between Figs. 4 and 6 is less dramatic than the change in spreading for

the FS scheme in Figs. 3 and 5, the FS results still exhibit less diffusion than the FV results,

comparing Figs. 5 and 6.

Employing the compressive Sut)erbee limiter with the FS scheme yields the results of

Fig. 7. In this case the discontinuity is confined to a 2-3 cell stencil, and does not grow in

space. Applying the Superbee limiter to FV cannot eliminate all artificial dissipation for

this case, as is possible with FS. The FV results (not shown) corresponding to Fig. 7 spread

the discontinuity over four (:ells by the outflow boundary, with a continually broadening

trend.

However, while it is possible to use the Superbee limiter with FS for this case, compressive

limiters can be unstable on different grid orientations. For example, no degree of compression

is stable for the case of Fig. 3. This potential for instability is related to global positivity,

as discussed by Sidilkover and Roe. 3

The effect of using a general unstructured grid is investigated in Figs. 8 and 9. The

unstructured grid in this case was generated using VGRID. ls'19 The FS solution exhibits

less dissipation, but is not as smooth as the FV solution. While the FS scheme preserves

contact discontinuities over larger spatial ranges than the FV scheme, FS does not appear

to degenerate gracefully with regard to extreme coarsening of the unstructured mesh for

this test case. This behavior could have negative implications for applications employing

multigrid convergence acceleration.
Circular Advection

Circular advection is achieved by setting _ = (y, - c). A decaying sine-wave input profile

is used,

u(x, 0) = (el sin 7_x) 2

Results for tile two schemes, using the Minmod limiter, are presented on the worse-case

cut-cartesian mesh in Figs. 10 and ll. Again, the FS results are considerably less diffusive

than tile FV solution.

Tile circular-adveetion problem is also applied on an unstructured mesh. The input

profile for this case consists of both a top-hat function and a decaying sine wave, allowing

comparisons between tim schemes for both sharp discontinuities and smooth gradients. The

input profile is,

{ (c 2xsin(27rx)) 2 -0.5_<x<0

u(x, 0) = 0 -0.6 < x < -0.5

0.4 -0.8 _< x < -0.6

0 -1 _< x < -0.8

Results for this case are displayed in Fig. 12 for FS and Fig. 13 for FV, both using the

Minmod limiter. FS performs significantly better at preserving the top-hat distribution.

FS also does a better job of maintaining the minimum and maximum values of the sine

distribution, though both schemes do well on the smooth gradient portion of the sine wave.

Non-linear Advection

The non-linear adveetion equation is obtained fror:l Eqn. 1 by setting F = (_, u) with

u = 0. In non-conservative form the equation is written,

ttt + _t_tx + uy =0

A coalescing shock problem is considered, with an anti-syinmetric input profile,

= = 0

u(x, 0) = -2x- 1 on x= (-1,0)

The exact solution to this problem contains symmetric expansion fans on the sides and a

1 1

compression fan at the inflow that coalesces into a vertical shock at (x, y) = (-i, 5)"

The first mesh is cut-cartesian containing 26 × 2{; nodes. The FS and FV solutions,

both using the Mininod limiter, are presented in Figs. 14 and 15. Both algorithms exhibit

the same grid dependence on the amount of artificial dissipation as seen before, with the

left-half solutions having more diffnsion than the righ; halves, due to the grid orientation.

Both inethods perform the same in the compression-fro region, coalescing into a shock to

within the accuracy of the input-profile diseretization.

The shock is more sharply defined by FS than by FV. Figure 14 has the correct shock

speed, with nearly the entire gradient captured in on,_ cell thickness. In contrast, Fig. 15

10
showsa slightly incorrect shock speedwhen using FV, as tile shockprogressesto tile left

beyondthe coalescencepoint, eventhough the discretization is conservative.Tile incorrect

shockspeedresults from a non-symmetricdistribution of the dependentvariable to the left

and right of the shock,causedt)y tile excessiveartificial diffusion generatedon the gridmisaligned(left-hand) side.

Contoursof the absolutevalueof the error arepresentedin Figs. 16and 17. Errors from

both computed solutions showa lack of symmetry, again reflectingthe grid dependenceof

the artificial diffusion terms. The error levelsfrom FS are lessthan from FV. The shock

curvature in the FV solution at the coalescing point is clearly visible in Fig. 17, resulting in

significant downstream errors in the shock location as compared with the FS errors.

This problem is repeated on a 25 x 25 mesh with symmetric diagonal cuts, favorably

aligned with the advection directions. The FS and FV solutions, Figs. 18 and 19, are in

good agreement. Plots of the absolute error contours, Figs. 20 and 21, show FS to be a

little more accurate than FV for this case.

The final mesh for this case is a truly unstructured triangulation containing 847 nodes

and 1617 cells. The nodes are clustered to the outflow boundary, with a bias towards the

left-hand side. The FS solution is presented in Fig. 22, showing very accurate and crisp

shock resolution and good symmetry in the solution contours despite the mesh-clustering

bias. In contrast, the FV solution in Fig. 23 has a more-diffuse shock and again an incorrect

The FV solution is also shock speed, with the outflow shock offset to the left, of z = -_.

somewhat less symmetric than the FS solution.

Linear Diffusion

Choosing F = 0 the heat-conduction equation is obtained from Eqn. 1,

ut = V. (vVu)

The test problem, a steady-state boundary value problem on a unit square, is taken from

Tomaich. 4 The Dirichlet boundary values are,

u(-1, y)=0, u(O,y) = sin(Try)

u(z,O) =0, u(z, 1) = -sin(rr.r)

The analytical solution on x = [-1, 0], y = [0, 11 is,

1

u(x, y) - sinh 7r [sinh(rr(x + 1)) sin(roy) + sinh(rry)sin(rr(x + 1))]

Both diffusion discretizations, Eqns. 10 and 29, are compared on a 438-node unstructured

mesh. Figures 24 and 25 plot the absolute value of the error in the converged solutions

using Eqns. 10 and 29, respectively. A carpet plot of the solution, using the finite element

formulation, is presented in Fig. 26.

The finite element treatment is clearly more accurate, and is used to discretize the diffusion terms for both FV and FS in the following section. The average-gradient results in

Fig. 24 appear to exhibit a decoupling mode, similar to odd/even deconpling for structured

meshes.

11
Linear Advection-Diffusion

The final test case is a linear advection-diffusion problem of Smith and Hutton. 2° The

flux function is F = Au, with,

A: (2y(1- x2), -2x(1 - y2))

The streamlines for this problem, while not truly circular, are similar in orientation to the

circular advection problem. Tile inflow profile is,

u(x,0)=l+tanh(20x+ 10)

The diffusion coefficient is chosen to be a constant, u = 10 -3. The domain is the unit square

in the second quadrant. No closed-form solution is known to the authors for this problem.

A sequence of five unstructured meshes is considered. The meshes have no preferred

clustering or stretching and have nominal node-spacings of 0.1, 0.05, 0.025, 0.0125, and

0.00625, labeled as Meshes A-E, respectively. The number of nodes for each mesh, along

with tile solution times for both FS and FV on a 195 MHz SGI R10000 CPU are listed in

Table 1.

L2-norms of the artificial and physical viscosities computed using both FS and FV are

presented for each mesh in Table 2. Notice that the norm of the artificial dissipation

for both FV and FS drops lower than the norm of the physical dissipation on Meshes D

and E. Since tile algorithms select only the larger of the physical or artificial dissipation

(Eqn. 30), Table 2 suggests both schemes are grid resolved on Mesh D. However, the norm

of the physical dissipation is smaller for FV than FS on each mesh A D. The physical

viscosity is driven by the solution curvature, suggesting FS maintains the solution profile

sharper than FV on the coarser meshes. A comparison of outflow profiles will soon veri_

this interpretation.

Further evidence of a grid-resolved FS solution is seen in Figs. 27 and 28. The FS

solution on Mesh E at tile outflow boundary is presel ted along with tile inflow profile and

the corresponding pure-advection (v = 0) FS solution in Fig. 27. Tile pure-advection

solution is seen to replicate the inflow profile, with a clear separation from the diffused,

u = 1() -:_, solution. Plotting only the FS results with respect to grid refinement, Fig. 28

shows a convergence of the outflow profile by Mesh C for FS.

The accuracy of FS and FV are compared in Fig. 2.(t, where the outflow solutions from FS

and FV are plotted tbr Meshes C and E. Taking the _xid-converged FS Mesh-E solution to

be the "truth" solution, it is clear that FS reaches th( grid converged solution on a coarser

mesh than FV.

Computational efficiencies of the two algorithms a_e compared in Fig. 30, where the L2-

norm of the residual is plotted versus CPU time for the fine-mesh FS and FV solutions, along

with the FS convergence history on Mesh D. The Me_h-E FS solution converges in 760 sec.

The corresponding FV solution takes 2.5 times longer than FS, due, in part, to the need

to reconstruct gradient information at each node with FV for second-order spatial accuracy.

However, considering the solution time to reach a given accuracy, it is more reasonable to

compare the FS solution time on Mesh D to tile fines ;-mesh FV solution. The FS Mesh-D

solution took only 64 sec, a factor of 29 times less than FV on Mesh E, and still shows better

accuracy than tlle fine-mesh FV solution.

12
An even greater speedupis seenwith FS in conjunction with the van Albada limiter,

wherenow the Mesh-Bsolution over-plotsthe curvefrom tile finest grid, shownin Fig. 31.

The correspondingFV result using the van Albada limiter on Mesh B is included, and

clearlyfalls short of tile FSaccuracy.The FV casewasrepeatedwith the highly-compressive

Superbeelimiter with little improvementin accuracy.The solution time for FS on MeshB

is about one second,yielding a speedupfactor of 2 3 ordersof magnitude over FV.

The final set of resultsaddressesconvergenceissueswhile pushing the positivity limits.

Figure 32 comparestwo convergencehistories for the second-orderFS on MeshB. The

non-converging,though stable,convergencehistory is the result of using strict positivity argumentsto setthe timestep (Eqn. 38). The resulting solutionis boundedand approximately

correctbut oscillator3.'.Limiter "ringing" is consideredto be a contrit)utor to this behavior,

and the higher-orderdiscretization for the implicit matrix could be reducing the diagonal

dominance,and hencestability, of the Gauss-Siedeliteration.

Full convergenceis achievedby using first-order positivity coefficients,which are not

dependenton the limiters. The resulting local timestepswill not be as largeastrue secondorder positivity would allow, but appearto be morerobust.

Summary of Results

Fluctuation splitting and finite volume schemes are compared in detail as applied to scalar

advection, diffusion, and advection-diffusion problems. The fluctuation splitting scheme is

seen to introduce less artificial dissipation while treating advection terms, allowing for more

accurate resolution of weakly dissipative advection-diffusion problems. The ability to resolve

solutions to these problems on coarser meshes makes the fluctuation splitting scheme the

preferred choice over finite volume.

Linear advection test problems are utilized to investigate tile dependence of artificial

diffusion production on grid orientation. Both fluctuation splitting and finite volume are

shown to exhit)it grid dependencies, but with fluctuation splitting producing less artificial

dissipation on all grids considered.

A non-linear coalescing shock problein further explores grid dependencies as cases are

constructed that result in incorrect shock speeds for finite volume. Fluctuation splitting

shows correct shock speeds for all grids and provides tighter shock capturing than finite

volume.

An adveetion-diffusion problem with small physical dissipation (diffusion coefficient of

10 -a) is considered where the reduction in artificial dissipation with fluctuation splitting

results in a significant accuracy improvement over finite volume. Convergence times are

compared, showing a speedup of 2.5 for fluctuation splitting over finite volume on identical

grids, using a point Gauss-Seidel relaxation. However, a grid convergence study shows

fluctuation splitting has better resolution of the solution on a coarser mesh than finite volume

does on finer meshes, resulting is a speedup of 29 for fluctuation splitting over finite volmne.

Based upon these significantly reduced solution times for solving model problems, as compared to the current state-of-the-art finite volume method, fluctuation splitting is considered

a worthwhile scheme to pursue for modeling fluid dynamic problems.

13
References

lBarth, T. J., "Aspects of Unstructured Grids a'ad Finite-Volume Solvers for the Eulet and Navier-Stokes Equations," Computational Fluid Dynamics, No. 1994-04 in Lecture

Series, von Karman Institute for Fluid Dynamics, 1994.

2Mesaros, L. M., Multi-Dimensional Fluctuation Splitting Schemes for the Euler Equations on Unstructured Grids, Ph.D. thesis, University of Michigan, USA, 1995.

3Sidilkover, D. and Roe, P. L., "Unification of Some Advection Schemes in Two Dimensions," Report 95-10, ICASE, Hampton, Feb. 1995.

4Tomaich, G. T., A Genuinely Multi-Dimensional Upwinding Algorithm for the NavierStokes Equations on Unstructured Grids Using a Compact, Highly-Parallelizable Spatial Discretization, Ph.D. thesis, University of Michigan, USA, 1995.

5Carette, J.-C. and Deconinck, H., "Adaptive Hybrid Remeshing and SUPG/MultiDUpwind Solver for Compressible High Reynolds Number Flows," AIAA Paper 97-1857, May

1997.

6Efraimssom G., "A 2D Analysis of the Influence of Artificial Viscosity Terms on Sohltions of tile Euler Equations," Journal of Computational Physics, Vol. 138, No. 1, Nov.

1997, pp. 103.

7Wood, W. A. and Kleb, W. L., "Diffusion Characteristics of Upwind Schemes on Unstructured Triangulations," AIAA Paper 98-2443, Jun. 1998.

SBathe, K.-J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Inc.,

Englewood Cliffs, USA, 1982.

9Roe, P. L., "Approximate Riemann Solvers, Parameter Vectors, and Difference

Schemes," Journal of Computational Physics, Vol. 43, October 1981, pp. 357-372.

l°Anderson, W. K. and Bonhaus, D. L., "An Implicit Upwind Algorithm for Computing

Turbulent Flows on Unstructured Grids," Computers and Fluids, Vol. 23, No. 1, ,]an. 1994,

pp. 1 21.

LIDeconinck, H., Roe, P. L., and Struijs, R., "A Multidimensional Generalization of Roe's

Flux Difference Splitter for the Euler Equations," Co'_nputers and Fluids, Vol. 22, No. 2/3,

1993, pp. 215 222.

12Sidilkover, D., "A Genuinely Multidimensional U:)wind Scheme and Efficient Multigrid

Solver for the Compressible Euler Equations," Report 94-84, ICASE, USA, Nov. 1994.

l_Chakravarthy, S. R. and Osher, S., "High Resolution Applications of the Osher Upwind

Scheme for the Euler Equations," AIAA Paper 83-1943, Jun. 1983.

14Bruner, C. W. S. and Walters, R. W., "Paralleli::ation of the Euler Equations on Unstructured Grids," AIAA Paper 97 1894, Jun. 1997.

lSvan Albada, G. D., van Leer, B., and Roberts W. W., "A Comparative Study of

Computational Methods in Cosmic Gas Dynamics," Report 81-24, ICASE, Hampton, Aug.

1981.

16Hirsch, C., Numerical Computation of Internal and External Flows--Volume 2: Computational Methods for Inviscid and Viscous Flows, John Wiley & Sons Ltd., 1990.

17Spekreijse, S., "Multigrid Solution of Monotone Second-Order Discretization of Hyperbolic Conservation Laws," Mathematics of Computati_,n, Vol. 49, 1987, pp. 135 155.

14
lSPirzadeh,S., "Three-DimensionalUnstructuredViscousGrids bytile Advancing-Layers

Method," AIAA Journal, \7ol. 34, No. 1, Jan. 1996, pp. 43-49.

19Pirzadeh, S., "Progress Toward a User-Oriented Unstructured Viscous Grid Generator,"

AIAA Paper 96 0031, Jan. 1996.

2°Smith, tR.. M. and Hutton, A. G., "The Numerical Treatment of Advection: A Performance Comparison of Current Methods," Numerical Heat Transfer, Vol. 5, 1982.

15
Table 1 Grids and solution times for advection-diffusion problem.

CPU seconds

Mesh Nodes FS FV

A 134 < 1 < 1

B 495 1 1

C 1,928 5 8

D 7,529 64 145

E 28,915 760 1880

16
Table 2 L2-norms (×105) of artificial and physical viscosities for advection-diffusion

problem.

FS FV

I1¢11_ I1¢,,11,,Mesh 11_'112II_ll_

(art.) (phys.) (art,.) (phys.)

1274 215 A 1918 190

597 265 B 640 176

192 161 C 144 119

54 76 D 46 66

13 36 E 18 36

17
ry

]_'23

2

1

X

Figure 1 Fluctuation splitting element nomenclature.

18
0.5

U=0

0.1

U=I

i_ X=(1,-1)

spacing:

0.05 - 0.95

0.1 increments

-0.5 0

X

Figure 2 First-order fluctuation splitting, uniform advection.

19
0.5

U=O

0_1

U=I

-0.5 0

X

z=(i,-i )

Contour

spacing:

0.05 - 0.95

0.1 increments

Figure 3 Second-order fluctuation splitting, uniform advection.

2O
Y

0.5

U=O

°-1 -o.5

U=I

o

X

Z=(1,-1)

Contour

spacing:

0.05 - 0.95

0.1 increments

Figure 4 Second-order finite volume, uniform advection.

21
Y

0.5

U=O

0.1

/ ///_L,

×

_,=(1,1)

Contour

spacing:

0.05 - 0.95

0.1 increments

Figure 5 Second-order fluctuation splitting, uniform advection.

22
y ,I

0.5

U=O

0-1 -().5

U=I

/

X

_.=(1,1)

Contour

spacing:

0.05 - 0.95

0.1 increments

Figure 6 Second-order finite volume, uniform advection.

23
Y

0.5

U=O

0.1 -().5 r _ ; , J

U=l 0

X

_.=(1,1)

Contour

spacing:

0.05 - 0.95

0.1 increments

Figure ? Second-order fluctuation splitting with compressive limiter.

24
Contour

spacing:

0.05 - 0.95

0.1 increments

= i i

0-I -0.5 0

U=I X

Figure 8 Fluctuation splitting on unstructured mesh.

25
1 --

0.5

U=O

0.1

/

4).5 _ ' 0

U=1 X

_.=(1,1)

Contour

spacing:

0.05 - 0.95

0.1 increments

Figure 9 Finite volume on unstructured mesh.

26
Y

U=0

0.1

0.5

f

f

J

j_J ff

-0.5 0

U = (e x sin _x) 2 X

X=(y,-x)

Contour

spacing:

0.025 - 0,375

0,05 increments

Figure 10 Fluctuation splitting, circular advection.

27
Y

0.5

U=O

1- f

yf

f

/_J • /j _

f

_ _,= (y,-x)

Contour

:- _ - ' 0.025 - 0.375

0.05 increments

0-1 -0.5 0

U = (e =sin _x) 2 X

Figure 11 Finite volume, circular advection.

28
0.5

U=0

1

. , _, , - ,

;_=(y,-x)

Contour

spacing:

0.025 - 0.375

0.05 increments

0-1 -0.5

U = (e _ sin(2_x)) 20 U=O U=0.4 U=O X

Figure 12 Fluctuation splitting on unstructured mesh, circular advection.

29
0.5

U=0

0.1

U=0

Z=(y,-x)

_ Contour

spacing:

0,025 - 0,375

_'i_'_ 0.05 increments

i

-0.5

U = 0.4 U = 0 X U = (¢ _' sin(2_x)) 20

Figure 13 Finite volume on unstructured mesh, circular advection.

3O
0.5

Contour

spacing:

-1 -1

0.1 increments

U=O U=O

0-1 -0.5 0

U = -2x-1 X

Figure 14 Fluctuation splitting, Burgers equation.

31
v iii_i i

o.s !!:i! :

U--O 0

U = -2x -I

-0.5

X

i

o

Contour

spacing:

-1 - 1

0.1 increments

U=O

Figure 15 Finite volume, Burgers equation.

32
1

y_ ii ii ! i

0._ I, ; , ,, J . | , J , * I I

i!_! i i o_iii ii'i io !

×

Contour

spacing:

0-1

0.1 increments

Figure 16 Fluctuation splitting, Burgers equation, absolute error.

33
0.5

0.1 -0.5 0

X

Figure 17 Finite volume, Burgers eqtation, absolute error.

34
0.5

Contour

spacing:

-1 -1

0.1 increments

U=O U=O

°-1 -0.5 o

U = -2x - 1 X

Figure 18 Fluctuation splitting, Burgers equation, symmetric mesh.

35
I

Y

Contour

spacing:

-1 -1

0.1 increments

U=O U=O

Figure 19 Finite volume, Burgers equation, symmetric mesh.

36
0.5

1

Contour

spacing:

0-1

0.1 increments

- i 7 .... • i ii ,•" i ,• • ," iii

0-1 -0.5 0

X

Figure 20 Fluctuation splitting, Burgers equation, absolute error.

37
Y

0.5

01 0s

X

Contour

spacing:

0-1

0.1 increments

Figure 21 Finite volume, Burgers equation, absolute error.

38
1-1 , f . , -0.5 ..... , _ . , 01

. : i .'-i 'i:'' :i',i i', Contour

• ,": ' i" _" '::_' _ .' " ,'," :, spacing:

, • =- " : " _ , -I-I

i i 0.1 increments

0.5 _: : • ' _ ' i 0.5

U=O o o U=O

-0.5

U = -2x -1 X

Figure 22 Fluctuation splitting, Burgers equation, unstructured mesh.

39
-1 -0 5 o

1 _ _ , _..:i:,_' i :!: . i..:,, :,,, \_ , 1

::_:, ! i!:_ i:_!_: !!ii

• :, •.-_ . _ .:' , .:, spacing:

, ' " . ' " _,: -1-1

:,_ :: ;- : .' " . , ": 0,1 increments

0.5 : • • ' _ " ' :', 0.5

u=o :. ....... ...... u=o

0 0

-0.5

U = -2x - 1 X

Figure 23 Finite volume, Burgers equation, unstructured mesh.

4O
1-1

0.5

0.1

, • -0.5 , , 01

_ ____-]_ o-o02

"_ "' '___/('_Ilt_ 0.002 incr.

._ 05

×

Figure 24 Pure=diffusion problem error, diffusion terms from Eqn. 10.

41
0,5

I -0.5 01

0. I I ' J '-0.5 ' ' '

K

Contour

spacing:

0 - 0.02

0.002 incr.

0.5

00

Figure 25 Pure-diffusion problem error, diffusion terms from Eqn. 29.

42
_).25

Figure 26 Heat equation solution using finite element formulation.

is 0.1.

Contour increment

43
2 -- -- _

FS, v=lO "3

1.75 " " FS, v=O

Inflow profile

I

0.25

"_ ".-,

' i i i I i i i _ L i L i = I L I J J .J 0.4 0.45 0.5 0.55 0.6 0.65

Y

Figure 27 Fluctuation splitting profiles on finest nmsh, advection-diffusion problem.

44
U

 eSh

_ _ A

i \ - - E

1.25 \\\\_

0.75 - _\ \

, \\

0.5 _" "-,\

8._.... o'._... o._.... o_._ .... o.__,-o!o__--_._

Y

Figure 28 Fluctuation splitting grid convergence, advection-diffusion problem.

45
U

2

FS, Mesh E

1,75 _,_ FS Mesh C

_ .......... FV' Mesh E

15 "_ ----- FV:MeshC

\

0.5 - _

i O.25

O I L _ J J I J I I L I I I _ I L I I i I I , I I , 1 I "l _

0.35 0.4 0.45 0.5 ( 55 0.6 0.65

Y

Figure 29 Fluctuation splitting and finite volume for advection-diffusion problem.

46
100

10-1

I Res 12 FV, mesh E

10-z

1o. \ \

/_m_e_o/ _,_e°_ \

10-5 _ I l

lO'e 0 i, _ L ' EO0 ' I ' IO_)L L , 1500 ' "_'2000

CPU(sec)

Figure 30 Convergence histories for advection-diffusion problem.

47
2__

1.75 "_ "_ _ _ FS, minmod, mesh E

FS, van Albada, mesh B

FV, van Albada, mesh B

1.5

1.25 \\

1

0./5

0.5

0.25 _

o-- .... I .... I .... I ....

035. 0.4 0.45 0.5 C55 0.6 0.65

Y

Figure 31 Advection-diffusion results using van Albada limiter.

48
REPORT DOCUMENTATION PAGE FormApproved OMB No. 0704-0188

Pul_ic reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, seamhing existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Sand comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis

Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DO 20503.

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

October 1998 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Comments on the Diffusive Behavior of Two Upwind Schemes WU 242-80-01-01

6.AUTHOR(S)

William A. Wood and William L. Kleb

7.PERFORMINGORGANIZATION NAME(S)ANDADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

9.SPONSORING/MONITORING AGENCYNAME(S)ANDADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

8. PERFORMING ORGANIZATION

REPORT NUMBER

L-17777

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA/TM-1998-208738

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Category 64 Distribution: Standard

Availability: NASA CASI (301) 621-0390

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

The diffusive characteristics of two upwind schemes, multi-dimensional fluctuation splittingand locally

one-dimensional finite volume, are compared for scalar advection-diffusion problems. Algorithms for the two

schemes are developed for node-based data representation on median-dual meshes associated with unstructured

triangulations in two spatial dimensions. Four model equations are considered: linear advection, non-linear

advection, diffusion, and advection-diffusion. Modular coding is employed to isolate the effects of the two

approaches for upwind flux evaluation, allowing for head-to-head accuracy and efficiency comparisons. Both the

stability of compressive limiters and the amount of artificial diffusiongenerated by the schemes is found to be

grid-orientation dependent, with the fluctuation splittingscheme producing less artificial diffusion than the finite

volume scheme. Convergence rates are compared for the combined advection-diffusion problem, with a speedup

of 2.5 seen for fluctuation splittingversus finite volume when solved on the same mesh. However, accurate

solutions to problems with small diffusion coefficients can be achieved on coarser meshes using fluctuation

splitting rather than finite volume, so that when comparing convergence rates to reach a given accuracy, fluctuation

splitting shows a speedup of 29 over finite volume.

14. SUBJECT TERMS

CFD, Fluctuation Splitting

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

19. SECURITY CLASSIFICATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

54

16. PRICE CODE

A04

20. LIMITATION OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102
100

10-1

IResl,

10,2

10 _

IO_

10"s

10 -6

Second-order positivity

coefficients

25 50 75 100

Iteration

Figure 32 Convergence rates using first- and second-order positivity coefficients.

49