NASA Technical Memorandum 102800

Extension of a Streamwise

Upwind Algorithm to a

Moving Grid System

Shigeru Obayashi

MCAT Institute, San Jose, California

Peter M. Goorjian and Guru P. Guruswamy

Ames Research Center, Moffett Field, California

April 1990

NationalAeronautics and

Space Administration

Ames Research Center

Moffett Field, California 94()35-1000

EXTENSION OF A STREAMWISE UPWIND ALGORITHM

TO A MOVING GRID SYSTEM

Shigeru Obayashi,* Peter M. Goorjian, and Guru E Guruswamy

Ames Research Center

SUMMARY

A new streamwise up.rind algorithm has been derived to compute unsteady flow fields with the use

of a moving-grid system. T!ae temporally nonconservative LU-ADI (lower-upper-factored, alternatingdirection-implicit) method has been applied for time-marching computations. A comparison of the temporally nonconservative metho(_ with a time-conservative implicit upwind method indicates that the solutions

are insensitive to the conservative properties of the implicit solvers when practical time-steps are used. Using this new method, compwations have been made for an oscillating wing at a transonic Mach number.

The computed results confinn that the present upwind scheme captures the shock motion better than the

central-difference scheme based on the Beam-Warming algorithm. The new upwind option of the code

allows larger time-steps and thus is more efficient, even though it requires slightly more computational

time per time-step than the c_;ntral-difference option.

INTRODUCTION

A code, ENSAERO, is being developed at Ames using the Euler/Navier-Stokes equations for computing the unsteady aerodynamics and aeroelasticity of aircraft. The capability of the code has been demonstrated by computing vortical and transonic flows over flexible swept wings (refs. 1 and 2). The flow fields

were calculated by a time-ac,zurate, finite-difference scheme based on central differencing.

The purpose of this ,:tudy is to enhance the algorithm capability of the present code. In this respect, the use of a new upwind scheme in comparison to the current central-difference (CD) scheme is

investigated. The CD scheme requires an artificial dissipation to stabilize computations. Such artificialdissipation models lead to more dissipative solutions than upwind schemes. In addition, the CD scheme

is sensitive to the amount of dissipation and needs a specific dissipation coefficient for each case. On the

other hand, upwind schemes do not require that any coefficient be specified.

Recently, a streamwi.,e upwind algorithm has been developed and applied to treat steady-state problems of transonic flows over wings (ref. 3) and vortical flows over a delta wing (ref. 4) on fixed grids. The

main feature of the streamwise method is the use of the local stream direction, flow velocity, and pressure

gradient. The switching of flux evaluations always takes place at sonic values, where shock waves may exist. Therefore, this method fcllows the flow physics more closely than conventional upwind methods based

* MCAT Institute, San Jose, California.
ondimensionalsplitting. Thecomputedresultsconfirmthehigherresolutionof thepresentalgorithmover

theCD scheme, as well as over other upwind schemes.

In this paper, the streamwise upwind algorithm has been extended from fixed coordinates to moving

coordinates for computing flows over moving components. The streamwise upwind scheme was applied to

steady-state problems by using the lower-upper-factored, alternating-direction-implicit (LU-ADI) method

in order to accelerate convergence (ref. 4). The same LU-ADI method, which is first-order accurate but

nonconservative in time, is used in the present unsteady computations. In order to check the validity of the

LU-ADI method for computations over moving grids, a conservative implicit version of the streamwise

upwind scheme is also considered. The resulting algorithm has been implemented in the code as a finitevolume, upwind option, in addition to the previous CD option. The updated code has been successfully

applied for computing unsteady transonic flows over an oscillating wing. The computed unsteady pressures

are compared with the experimental data.

GOVERNING EQUATIONS

The nondimensionalized thin-layer Navier-Stokes equations used in this study can be written in

conservation-law form in a generalized body-conforming curvilinear coordinate system ( _, _7, _ ) as

follows:

1

aTQ + o':3,_#+ a,7/_ + a_;G = _..-_-oq(G-----3t' (l)

/-._e

The Euler equations are obtained by setting the viscous flux vector G_ equal to zero. The vector of conserved quantities Q and the inviscid flux vectors/_, F, and G are

r 2 0 l

I p,,u+_,pI

' I + 1

# = 7 1p,,,u+_,pI

LpHU - (tpJ

A

2v

[puy+n_pl

1/,vv /

?= 7/p_ +,7,r'/

LpHV - '0tpJ

A

._w

[ puW + Gp ]

I p,,,w+<,pI

LpHW -- (tp J

where H is the total enthalpy, H = (e + p)/p, and the contravariant velocity components U, V, and

are defined as

= r/t + 'r/=u + _/yv + _Tzw

W = _t + Gu + _v + Gw

2
TheCartesianvelocitycomtxmentsu, v, and w are nondimensionalized by the free-stream speed of sound

coo; the density p is nondimensionalized by the free-stream density Poo; and the total energy per unit volume

e is nondimensionalized by P,ocoo. 2 The viscous flux vector G v is given by

with

#mira3

_mlW(

+

/_mlV¢

_m2((_.u

0+

+

_ra2(tl

_ra2(z

+ (_v + (zw)

1

= + +

1 U2 V2

_3 = "_ ( + + "11)2 )( +

1 ( C2 )(

Pr (_1 - 1)

where Re is the Reynolds number, Pr is the Prandtl number, c is the speed of sound, and J is the transformation Jacobian. Pressure is related to the conservative flow variables Q, through the equation of state

for a perfect gas:

P=(7-1) [ e-_-( P U2 + V2 +w 2) ] (2)

where p is the fluid density aad e is total energy per unit of volume of the fluid.

For the inviscid case. the viscous flux G_ is replaced by 0. For the viscous case, the viscosity

coefficient _ in _v is computed as the sum of #t + #t where the laminar viscosity/_t is taken from the

free-stream laminar viscosity, assumed to be constant for transonic flows, and the turbulent viscosity/_t is

evaluated by the Baldwin-Lomax model. 5

NUMERICAL ALGORITHM

The space-discretizec form of equation (1)can be written as

(3)

where a second-order central-difference evaluation is applied to the viscous term.

The evaluation of the inviscid fluxes is based on the finite-volume cell-centered scheme. To be

consistent with the finite-difference scheme in ENSAERO, the metrics are defined at each grid point where

the flow variables are stored. The surface vector of cell-interface, which is necessary for the finite-volume

formulation, can be obtained by averaging the metrics at the adjoining points. The free-stream preservation

of this metric evaluation was shown in reference 6.
Streamwise Upwind Algorithm

The formula of the streamwise upwind algorithm on the fixed grid system (_t = r/t = (t = 0) can

be written in the cell-interface flux for the r/-direction with a surface vector S normal to the (-( plane as

follows:

1

(Ft [1 + sign(Vi)cos201] + st A*(pqh eslcos2 0t

+ F, [1 - sign(V,)cosZ0_] - st A*(pq)_ e,,cos20_ (4)

Ap ) - [ IvlaQ + (c-IVl)( -p- e, + pav ed) ]sin20

where A *(pq) = p'q* - pq with local sonic values,

2

(q.)2 _ "7_1 (C2 + .L__q2)

p* = p((q*)2 '_ __--'r -7-/

and where q is the velocity magnitude, e8 = ( 1, u, v, w, H) T, ed = (0, k=, k_, kz, k=u + kvv + kzw) 7",

A • =., - .t, V = k=u + kvv + k,w, kz = "O=/IVrll, and so on. The basic scheme is first-order accurate with

I = j and r = j + 1. The variables are averaged unless defined as a difference between left and right states.

The switches Sl,r and the rotation angle 0 will be defined later. The cos 20 terms represent the projection on

the coordinate axis from the streamwise flux vector splitting, and the sin20 terms represent the projection

on the crossflow plane for those terms that use the flux difference splitting. More details for the formulas

are given in reference 4.

To extend equation (4) from a fixed-grid system to a moving-grid system, the flow velocity relative

to the fixed grid, q = (u, v, w), is redefined as the flow velocity measured relative to the moving grid,

q = ( u - xt, v - Yt, w - zt). This is consistent with the modification of the definition of the contravariant

velocity from a fixed-grid system,

Vfixed = _x _ + 711I) + _z _1)

= Vr/-q

to a moving-grid system,

_rmovi _ = tit + rlxU + rl_ V + rlzW

= rl=(u - zt) + r_(v - _t) + _z(w - zt)

= Vr/.q

where r/t = -r/zxt - r/_yt - rhzt. Thus, equation (4) can be rewritten for a moving-grid system by using

the following redefined variables: _2 = (u - xt) 2 + (v - yt) z + (w - zt) 2 and V = k=(u - xt) + k_(v -

!h) + kz( w Z zt): (Note that e, and ea remain unchanged because they originate in the eigenvectors of the

Jacobian OF / OQ.)

4
Equation(4)is writtenin avectorform. If theformulais rewrittenin componentform, thepresent

algorithmcanbesummarizec_with a surfacevectorS andamotionof its centroid,xt = (zt, Yt, zt), as

1

tV l( Fl+ + F( } (5) _(Qt, Q_, Sj+ _., xt]+_.) = 2 J

where

fi,_ ul,r + kxp" l,r

f+ "- l,r Vl,r + ICy Pl,r

f+l,r Wl,r + kz pt:t:r

f,+ Hi,,- kt p+ + (l_r, nlAp VmA2)sin20 |,r -- l,r

fi,_ = cpV)t,_ [1 ± sign(_,,) cosZ0l,,]

:t- 8t,_ A*(p_)z,_ cos20z,_ ± pt,_[Vmlsin20 - 1

_-A1 sin 20

p+ ,, [1 ± sign(_,_) cos20/,,] I

t,_ = l z,_ - _-A2 sin 20

where _ and _" are defined above, and where A1 = (c,, - IY l) mE= ( - I l)p A V, and

kt = -kxxt - kvyt - kzzt. The averaged state (m) is defined for p, u, v, w, and H by the arithmetic

average of the left (1) and right (r) states.

The switches sl and _,_ are defined in the manner of Godunov's method as follows. For _r > 0,

8z = 1 - Cl Cm

sr = (1 - era)(1 - er)

(6)

where

1 [l+sign(M 2

el,m,_ = _- l,_,_- 1)]

and M = _/c.

A simple way to evaluate the rotation angle is to use cos0 = V'/_. In supersonic flow fields,

however, it is important to detect whether the velocity projected to the grid line is beyond the Mach cone.

Thus, _r/_ is replaced by M • V/_ = V/c. If V/c becomes larger than one, cos0 is set to one. To avoid

expansion shocks, the rotation angle is determined by a mixture of averaged (m) and pointwise (l, r) values:

-2

COS2Ol,r = mini ( 1 - ¢) c-_'_+ w _2 , 1 ]

_l,r

(7)
The following relation is used here for evaluating _b6 [0,1] because of the smoothness:

[ 27 (1 - 2-_[q- 1 + (q + 1)P2]}, 0 ] = max _ pl

(8)

where pl and p2 denote upstream and downstream pressures, respectively. The sine is determined by an

arithmetic average of the cosines: sin20 = 1 - 1⁄2(cos20t + c0s20_).

Higher-order schemes are

primitive variables, p, u, v, w, and p. For example,

Pt = {1 + --_[(1 - _)V+ (1 + _)Al}pj

v_ = {1 CJ+'[( 1 + _)_7 + (1 - _)AI}pj+I

4

constructed from a one-parameter family, _, of interpolations of the

(9)

where V and A are backward and forward difference operators, respectively (ref. 7). For the third-order

scheme, _ = _-, Koren's differentiable limiter (ref. 8) is used here. The limiter _b is calculated as

3VpjAp] + e (10)

_bj = 2(Apj - Vpj) 2 + 3_TpyApj + e

where a small constant e, e = 10 -6 typically, is added to prevent the division by zero. The same formulas

are used for the other primitive variables.

LU-ADI Method

One of the time-marching methods for the present upwind scheme is the LU-ADI factorization

method proposed by one of the present authors (ref. 6). The LU-ADI method is a compromise between

ADI and LU factorization. Applied to equation (3), this method is written as,

(T_LADAUAT(1) x (T_LBDBUsTff 1) x (T(LcDcUcT_ 1) x A_) '_= AtR" (11)

where R represents the right-hand side of equation (3). The change of volume in time is neglected, assuming translational movements of the grids. On the left-hand side, the original ADI operator of the BeamWarming method is rewritten using the diagonal form (ref. 9) and the first-order accurate Steger-Warming

flux-vector splitting (ref. 10). For example, in the rl-direction,

A

I + zXtSnB

= T.( r_+ zxtV,3_ + zxt_9_)T_ _

= T,( I- At/5_lj + AtV,3_)(I + Atl3BIj)-_(I + At3+BIi + Atz_9[_)T; 1

= TnLnDnUBT_ I

(12)

This factorization is the approximate LDU (lower-diagonal-upper) factorization. This is more stable than

simple LU factorization since the diagonal element always has the absolute value of the eigenvalues.

6
Approximate Block ADI Method

The LU-ADI method described in the previous section is nonconservative in time owing to the

diagonalization. To investigate the significance of this temporal nonconservativeness, an alternative approach is considered here. A time-conservative method can be constructed using a block-tridiagonal solver

similar to the Beam-Warming method. Since true Jacobians of the numerical fluxes of the present upwind

algorithm are expensive to compute, approximate Jacobians are used here.

To construct an implicit method for the present upwind algorithm, it is easier to start from the vector

form (eq. (4)) rather than from the component form (eq. (5)). From equation (4), the first-order-accurate

flux can be rewritten as

1 IWllj±_

_# =2 7<"_1 { Fj [1 q- sign(r/j)cos20j] -t- sjA'(p_)jcos26jesj

+ i_±_lsin20 QJ + cj+_ -IVj±_I

(13)

_±_ sin20 e,j±_pj

4- pj::_. ( cj+_. - [_.+_: [) sin2 0 eaji _.Vj )

To simplify the formula further, the e_j+_, term in the second line of the above formula is replaced by e,y.

Then, the Jacobian 0F4-/0Q ,.:an be approximated as

/_ Jj IV'olj I - {

2 JJ±r

q-[SjA.(p_)jCOS20j + Cj±_ --[Vj±_lpjsin20] Mj (14)

-_ [_.±_lsin20 1 + pj±_ (cj±_ - I_.±_[)sin20 pJ )

where B = OF/OQ, I is the identity matrix,

M m

10pes 1

p oQ p

I 1 0 0 0 0 I

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

°_(uz+v2+w2) -("r-1)u -("t-1)v -(7-1)w "i

(ooooo) 0 k_ kxk_ kxk_ 0

N = ecl . ( O , k=, kv, k_, O) = 0 k_k= k_ k vkz 0

0 kzkx kzk v k 2 0

0 Vfk= Vfk v vfkz 0

7
and where V y = k:_u + kuv + kzw. The switches and the angles are evaluated identically to the explicit

formulas. Similar to the results discussed in reference 11, this approach yields a method that is less stable

than the LU-ADI method.

Neglecting the viscous terms, the approximate block ADI method can be written as

.(s + At + 4i- - )

×

\ An )

× A tR" \ )

..,..

(15)

In contrast to the LU-ADI method, the resulting method is conservative. However, both methods are firstorder accurate in time. This method requires twice as much CPU time as the LU-ADI method.

RESULTS

Time-accuracy is an essential requirement for aeroelastic computations. Numerical schemes used

for flow calculations in aeroelasticity must guarantee the correct calculation of the amplitude and phase

of unsteady pressures. In order to verify the time-accuracy of the present code, unsteady flows over a

rectangular wing undergoing a prescribed oscillatory motion are computed.

The test cases consider inviscid and viscous unsteady flows over a rectangular wing with a NACA

64A010 airfoil section and an aspect ratio of 4. Steady and unsteady measured data from wind-tunnel tests

of this wing are given in reference 12. The unsteady data are given for the case in which the rigid wing

is oscillating in the pitching mode, a(t) = s,_ - s0 sin(wt), about an axis at x/c,. = 0.5, where cr is

the reference chord length, and w is the pitching frequency in radians per second. The inviscid flow is

computed at Moo = 0.8 with a mean angle of attack sra = 0% pitch amplitude s0 = 1 °, and reduced

frequency k= 0.27 (k= wc,./Uoo). In addition, for the viscous flow, the Reynolds number based on the

wing chord is Re = 2 × 10 6 . The Baldwin-Lomax model (ref. 5) is used to compute the turbulent eddyviscosity coefficient for this case.

The C-H-type grid system is used for the present computations. The (, r/, and _ coordinates represent the chordwise, spanwise, and normal (to the wing surface) directions, respectively. For the inviscid

computations, coarse and fine grids are used. The coarse grid contains 91 x 25 × 25 points in the _, r/, and

if-coordinate directions, respectively, and the fine grid contains 151 × 25 x 34 points. These grids have

a minimum spacing normal to the wing of 0.01 c_. The viscous grid also contains 151 × 25 × 34 poipts

and has the minimum spacing normal to the wing of 0 00005 c_. This yields a value of y+ < 6 at th_ first

grid point above the wing surface. The grids are generated algebraically.

Unsteady computations are started from the corresponding steady-state solution. The stability of

unsteady computations is verified by increasing the numbers of time-steps per cycle until no change is
noticedin theunsteadyresults Theconvergenceof theunsteadycomputationsto aperiodicflow isverified

by comparingtheresultsbetweencycles.Forall casesconsideredhere,theresultsfor thethird cyclegive

pressureprofilesthat areidenticalto thoseof the second-cycleresults. Thus,the numericaltransientis

confirmedto disappearwithin two cycles.

The present upwind results are also compared with the results obtained by using the existing CD

method. Reference 13 discusses the CD results obtained from the same test case. The CD method uses

the diagonal inversions of the Beam-Warming method (ref. 9), and thus it is first-order accurate but nonconservative in time. Its artificial dissipation model consists of second- and fourth-order dissipation terms

controlled by the amount of the second-difference of pressure.

Inviscid Solutions

First, the inviscid coarse-grid (91 x 25 × 25 points) solutions are presented. Figure 1 compares the

computed and measured unsteady, upper-surface pressure coefficients for the real and imaginary parts of

the first Fourier component. Pressure coefficients are shown for various span locations. The computations

were done by using the upwind algorithm with the LU-ADI method (UP-LU). The solid, dashed, and dotted

lines in the figure indicate the results obtained using 1440, 720, and 360 time-steps per cycle of oscillation

(steps/cycle), respectively. (The time-steps per cycle rate of 1440 corresponds to A t _ 0.02 .) The solution

profiles obtained using 1800 ,,teps/cycle coincided with those using 1440 steps/cycle. Thus, the unsteady

pressure profiles converged with respect to time-step sizes at 1440 steps/cycle. In addition, because the

differences between the results with 720 and 1440 steps/cycle are not critical when making comparisons

with experimental data, the computation with 720 steps/cycle is acceptable for numerical efficiency.

Figure 2 shows the analogous CD results. In this case, the second- and fourth-order dissipation

coefficients are fixed at 0.25 and 0.01, respectively (denoted as CD(0.01)). These coefficients are the

values that were used in references 1, 2, 9 and 13. The CD results show less dependence on time-step

size than do the UP-LU resuh s, although the profiles are smeared out at the large-gradient region near the

shock wave and leading-edge regions. In order to check the dependence of the CD method on the amount

of numerical dissipation, co_,putations were tried with half the amount of dissipation coefficients (0.125,

0.005). At 360 steps/cycle, the computation became unstable with these coefficients. At 720 and 1440

steps/cycle, the computations were stable, but the resulting unsteady pressure profiles were still smeared

out.

Figure 3 shows the re,;ults using the approximate block ADI method applied to the present upwind

algorithm (UP-BL). The UP-BL method was not stable for computations using 360 and 720 steps/cycle, and

thus the computations were done using 1080, 2160, and 3600 steps/cycle. The solution profiles obtained

with 3600 steps/cycle finally showed good agreement with the UP-LU result converged at 1440 steps/cycle.

Figure 4 compares tht- three methods, UP-LU, UP-BL, and CD(0.01), fixing the number of timesteps per cycle at 1440. Both UP-LU and UP-BL results give similar profiles at the shock motion, although

the UP-LU method is nonconservative in time. Both upwind methods give better agreement with the

experimental data in the region of the shock motion than the CD method.
Theseresultsindicatethat thetemporallynonconservativeLU-ADI methodis valid for unsteady

computations,evenwith movingshockwaves,whenthenumberof time-stepspercycleis large(i.e.,small

A t). In the present case, any number of time-steps per cycle greater than 720 will give a reasonable result.

The CD method also uses the temporally nonconservative diagonal form. However, the CD result differs

from both the upwind results. This indicates that the solution depends on the numerical dissipation more

than on the time-conservative properties of the methods.

The coarse-grid computations took 18.0, 19.5 and 42.3 _sec per grid point per time-step for the CD,

UP-LU, and UP-BL methods, respectively, on a CRAY-YMP computer using a single processor. Although

the UP-BL method is the most expensive, its accuracy does not appear to compensate for the computational

inefficiency. Thus, the UP-BL method will be dropped in the following computations.

Next, the inviscid computations were done on a finer (151 × 25 × 34 points) grid. Since detailed

comparisons with respect to time-step sizes are done on the coarse grid, the results are highlighted in the

comparison of the UP-LU and CD solutions, as shown in figure 5. Both computations use 720 steps/cycle,

as suggested in the coarse-grid case, because the time-step requirement for the fine-grid case is expected

to be the same as that for the coarse-grid case (note that both grids have the same minimum spacing size

normal to the wing). The fine-grid computations took 15.5 and 16.8/zsec per grid point per time-step for

the CD and UP-LU methods, respectively.

For the CD method shown in figure 5, the two sets of the second- and fourth-order dissipation

coefficients (0.25, 0.01) and (0.125, 0.005), indicated as CD(0.01 ) and CD(0.005), respectively, are used to

illustrate the effect of the artificial dissipation. The CD solutions approach the UP-LU solution at 50% and

77% semispan locations as the dissipation coefficients are reduced. However, the shock profile becomes

too steep and the post-shock oscillation in the plot becomes deep at the root section. Further reduction of

the dissipation coefficients (to 0.0625, 0.0025) made the computation unstable. Compared with the coarsegrid solution shown in figure 4, the CD solutions also tend to converge to the UP-LU solutions because of

the grid refinement. Therefore, the UP-LU method is confirmed to give a less dissipative and thus more

accurate solution than the CD method.

Viscous Solutions

To demonstrate the capability of the UP-LU method for viscous computations, the same case has

been computed with the viscous terms. The viscous grid has the same number of grid points as the inviscid

fine grid, but with a smaller spacing in the _-direction normal to the wing. Because of the stiffness caused

by the clustered viscous grids, time-step sizes used in the CD computations are usually determined by

stability considerations rather than by time-accuracy (ref. 13). Two time-steps/cycle rates were chosen to

demonstrate the robustness and accuracy of the UP-LU method. The results are again highlighted in the

comparison of the UP-LU and CD solutions at specified time-step sizes.

Figure 6 compares the viscous solutions using the UP-LU, CD(0.01), and CD(0.02) methods _,,ith

720 steps/cycle. This number of time-steps per cycle is too small to obtain accurate and stable results for

the CD method, a finding also reported in rcference 13. The numerical instabilities in the CD results are

10
shownby oscillationsin thepressureprofilesof CD(0.01).However,theUP-LUresultshowsthattheUPLU methodrelaxesthisstabilityrequirement.Thus,time-stepsizesfor theUP-LUmethodarenotlimited

by stabilityconsiderationsevenfor theviscouscase.

Figure7 showsthecorrespondingresultswith 1440steps/cycle.This is the time-steprate suggestedin reference13. Here,theCD(0.01)solutionbecomesstableandagreesreasonablywell with the

experimentaldata.To checkwhethertheamountof dissipationis reasonable,thedissipationcoefficients

werereducedby half, aswasdonefor theinviscidcase,buttheCD(0.005)computationdivergedfor the

viscouscase.In figure7,andanalogousto theinviscidresultsshownin figure5,theCDsolutionsapproach

theUP-LU solutionasthedis,,;ipationcoefficientsarereduced.Again theUP-LU methodis confirmedto

give aless-dissipativeandthu_moreaccuratesolutionthantheCD method.Theviscoussolutionsin figure7 aresimilarto theinviscidsolutionsin shownin figure5 becausetheexperimentalflow field doesnot

containstrongviscouseffect(_'ef.12).

TheseresultsalsoindicatethattheCD methodissensitiveto bothtime-stepsizeandto theamount

of artificialdissipationadded.'Fherefore,usersof theCD methodarerequiredto findanadequatecombination of time-stepsizeanddissipationcoefficientsonacase-by-casebasis.In contrast,theUP-LU method

is robust,andusersdo nothaveto find anydissipationcoefficient.

For theviscouscase,the UP-LU computationrequires18.9_zsecper grid point per time stepat

a speedof 141MFLOPSon aCRAY-YMPcomputerusinga singleprocessor,andthe CD computation

requires17.1#sec at a speed or 132 MFLOPS. There is 11% increase of CPU time when the upwind option

of the code is used.

CONCLUSIONS

A new streamwise upwind algorithm has been derived to compute unsteady flow fields with a

moving grid system. The temporally nonconservative lower-upper-factored, alternating-direction-implicit

(LU-ADI) method has been applied to compute flow over an oscillating wing at a transonic Mach number.

A comparison of the temporally nonconservative method with a time-conservative version of the upwind

scheme indicates that the solutions are insensitive to the time-conservativeness of the implicit solvers when

practical time-step sizes are used. The temporally nonconservative upwind method was found to be more

stable and twice as efficient as lhe time-conservative block ADI scheme, even though the block ADI scheme

employed an approximate form of the Jacobians. Comparisons have been made between results obtained

with both upwind schemes, with experimental measurements, and with computed results obtained using

the existing central-difference method.

Comparisons with experimental data show that the upwind algorithm predicts the shock motion

better than the central-difference (CD) method. The CD solutions are also found to be sensitive to the

amount of numerical dissipation, and as a result, the dissipation coefficients must be specified case by case.

In comparison, the present upwind method does not require a dissipation coefficient. Thus, the method that

combines the streamwise upwind and LU-ADI methods is proposed for practical computations. This new

upwind method is incorporated in the aeroelastic code ENSAERO. The new upwind option of the code

allows larger time-steps and thus is more efficient, even though it requires slightly more computational

time per time-step than the CI_ option.

11
REFERENCES

.

.

.

*

°

.

.

.

o

10.

11.

12.

13.

Guruswamy, G. P.: Time-Accurate Unsteady Aerodynamic and Aeroelastic Calculations of Wings

Using Euler Equations. AIAA Paper 88-2281, 1988.

Guruswamy, G. E: Vortical Flow Computations on Swept Flexible Wings Using Navier-Stokes Equations. AIAA Paper 89-1183, 1989.

Goorjian, P. M.: A New Algorithm for the Navier-Stokes Equations Applied to Transonic Flows over

Wings. AIAA Paper 87-1121-CP, 1987.

Obayashi, S.; and Goorjian, P. M.: Improvements and Applications of a Streamwise Upwind Algorithm. AIAA Paper 89-1957-CP, 1989.

Baldwin, B. S.; and Lomax, H.: Thin Layer Approximation and Algebraic Model for Separated

Turbulent Flows. AIAA Paper 78-257, 1978.

Obayashi, S.: Numerical Simulation of Underexpanded Plumes Using Upwind Algorithms. AIAA

Paper 88-4360-CP, 1988.

Anderson, W. K.; Thomas, J. L.; and van Leer, B.: A Comparison of Finite Volume Flux Vector

Splittings for the Euler Equations. AIAA Paper 85-0122, 1985.

Koren B.: Upwind Schemes, Multigrid and Defect Correction for the Steady Navier-Stokes Equations. Proceedings of l lth International Conference on Numerical Methods in Fluid Dynamics,

June, 1988.

Pulliam, T. H.; and Steger, J. L.: Recent Improvements in Efficiency, Accuracy, and Convergence

for Implicit Approximate Factorization Algorithms. AIAA Paper 85-360, 1985.

Steger, J. L.; and Warming, R. E: Flux Vector Splitting of the Inviscid Gas-Dynamic Equations with

Application to Finite-Difference Methods. J. Comp. Phys., vol. 40, 1981, pp. 263-293.

Liou, M-S.; and van Leer, B.: Choice of Implicit and Explicit Operators for the Upwind Differencing

Method. AIAA Paper 88-0624, 1988.

Mabey, D. G.; Welsh, B. L.; and Pyne, C. R.: A Summary of Measurements of Steady and Oscillatory

Pressures on a Rectangular Wing. Aeronaut. J. R. Aeronaut. Soc., Jan. 1988.

Chaderjian, N. M.; and Guruswamy, G. P.: Unsteady Transonic Navier-Stokes Computations for an

Oscillating Wing Using Single and Multiple Zones. AIAA Paper 90-0313, 1990.

12
o_

i

m

15

10

-5

15

10

i s

0

-5

15

10

5

0

n.

-5

15

UP-LU, 1440

...... UP-LU, 720

......... UP-LU, 360

i- Exp, Ref 12

-5

5

Semlspan

94%

1 I

I -5

5

77%

I f

t i -5

13

|

n0

-5 i -5

0 5 1.0 0

XIC r

-_o

• e

Root

I

.5 1.0

X/Cr

Figure 1. Effect of time-step sizes on inviscid upper-surface unsteady pressures for the rectangular wing

computed by using the upwi_d LU-ADI method: 91 × 25 × 25 grid, Moo = 0.8, c_m= 0 °, k = 0.27, and

oL0-- 1°.

13
15

10

5

rr

0

-5

15

10

i s

QC

0

-5

15-

10

eL

o

"6 5

(I) QC

-5

15

10

0

'3 S

0

I CD (0.01), 1440

CD (0.01), 720

......... CD (0.01), 360

O Exp, Ref 12

I I

L

o

m

c

I

_E

-5

S

Semispan

94%

I I

i o

I I -5

5

I I

f

= 0

-5

S

77%

I J

3

50%

I I

I

L

o

m

¢:

el

_s

I -5

Root

|

-5

0 .5 1.0 0 .5 1.0

X/C r X/C r

Figure 2. Effect of time-step sizes on inviscid upper-surface unsteady pressures for the rectangular wing

computed by using the central-difference method: 91 x 25 x 25 grid, Moo = 0.8, o_m = 0 °, k = 0.27,

and so = 1o.

14
15

10

a_

o

0

-5

15

10

0

-5

15

10

o.

o

-6 5

Q

Pr

-5

15

10

o.

o

"_ 5

0

UP-BL, 3600

UP-BL., 2160

......... UP-BL, 1080

[] Exp, Ref 12

I I -5

5

'=

Semispan

94%

I I

O.

o

I J -5

i 5

I I

I

-5

5

IL

(J

o

.=

_e

77%

I I

-I

50%

I I

Rool

-5 I -5 I J

0 .5 1.0 0 .5 1.0

X/C r X/C r

Figure 3. Effect of time-step sizes on inviscid upper-surface unsteady pressures for the rectangular wing

computed by using the upwind block ADI method: 91 x 25 x 25 grid, Moo = 0.8, a,_ = 0% k = 0.27,

and a0 = 1o.

15
15

10

o_

nO

-5

15

10

s

0

-5

15

10

s

m

0

-5

15

10

L

i s

rtUP-LU, 1440

...... UP-BL, 1440

......... CD (0.01), 1440

[] Exp, Ref 12

i 0

I I -5

5

L

U

I I -5

5

_ _ I"J__- "._"- ":.2o-----

z_l z..i ¢-i z..i z-i

_mispan

_%

I I

77%

I I

I I

-5

5

0

I I -5

5O%

I I

.....

Root

I I

-5

0 .5 1.0 0 .5 1.0

X/C r YJC r

Figure 4. Comparison of inviscid upper-surface unsteady pressures computed by using the upwind LUADI, upwind block ADI, and central-difference methods with 1440 time-steps/cycle: 91 x 25 x 25 grid,

Moo = 0.8, o_,.,,= 0% k = 0.27, and a0 = 1 °

16
15-

10

av

0

-5

15

5

UP-LU, 720

...... CD(0.005),720

......... CD (0.01), 720

[] Exp, Ref 12 0¢

I J -5

5

10

i s

nO

-5 I J

15-

IL

0

__ o

131

11l

_s

-5

5

¢L

o

0

'°+i .j

[] 0

0 r"

-5

,L-J

J

_'o

-5

"1 Semlsp-n

94%

I I

r

77%

I I

L LJ

50%

"+ I I

L

u

0

15

10

i 0

I I

.5 " 1.0

X/C r

i

i

I s

I1'.

1

0 .5

_C r

Root

-5 -5 ;

0 1.0

Figure 5. Comparison of inviscid upper-surface unsteady pressures computed by using the upwind LUADI method and the central-d:+fference method with two sets of dissipation coefficients: 151 × 25 × 34

grid, Moo = 0.8, a,n = 0 °, k =: 0.27, and a0 = 1 °.

17
15-

10

o"

i 5

0

-5

15

10

i s

0

IT.

-5 I I -5

15 5

10I

5 - ,,.. __ 0

0 -

-5 I 1 -5

15

10

0

Root

-5 I -5 I I

0 .5 1.0 0 .5 1.0

X/C r X/C r

Figure 6. Comparison of viscous upper-surface unsteady pressures computed by using 720 time-steps/

cycle: 151 × 25 × 34 grid, Re = 2 x 10 6, Moo = 0.8,o_,,, = 0°,k=0.27, anda0 = 1°.

18
15

10

o

s

re0

-5

15

10

rr

-5

15

10

L

0

-_ 5

0

re"

-5

15

10

O.

0

s

rr

5

UP-LU, 1440

CD (0.01), 1440

......... CD (0.02), 1440

I O__ Ex._p, Ref 1__2 _ 0

I J -5 I I 94% I

F

--L.J W

0

I I -5

["] 0 - --

I I

77%

I 1

O.

0

:1 o

Tr

_e

-5

50%

I I

,°f

-5 I j.l I -5

0 .5 1.0 0 .5

X/C r XIC r

Root

1 I

1.0

Figure 7. Comparison of viscous upper-surface unsteady pressures computed by using 1440 timesteps/cycle: 151 × 25 x 34 _,,-rid,Re = 2 × 10 6, Moo = 0.8, am = 0 °, k = 0.27, and a0 = 1 °

19
N/ A Nalion_ Ae_onau_cs end

SI_o_ Ad_in4l_,m_on

1. Report No.

NASA TM- 102800

4. Title and Subtitle

Report Documentation Page

2. Government Accession No. 3. Recipient's Catalog No.

5. Report Date

Extension of a Streamwise Upwind Algorithm to a Moving Grid

System

7. Author(s)

Shigeru Obayashi,* Peter M. Goorjian, and Guru P. Ouruswamy

9. Performing Organization Name and Address

Ames Research Center

Moffett Field, CA 94035-1000

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Washington, DC 20546-0001

April 1990

6. Performing Organization Code

8. Performing Organization Report No.

A-90101

10. Work Unit No.

505-60

11. Contract or Grant No.

13. Type of Report and Period Covered

Technical Memorandum

14. Sponsoring Agency Code

15. Supplementary Notes

*MCAT Institute, San Jose, California

Point of Contact: Shigeru Obayashi, Ames Research Center, MS 258-1

Moffett Field, CA 94035-1000 (415) 604-4982 or FFS 464-4982

16. Abstract

A new streamwise upwind algorithm h_sbecn derived to compute unsteady flow fields with the use

of a moving-grid system. The temporally nonconservative LU-ADI (lower-upper-factored, alternatingdirection-implicit) method imsl_cn applied for time-marching computations. A comparison of the

temporally nonconservative method with a time-conservative implicit upwind method indicates that the

solutions are insensitive to the conservative properties of the implicit solvers when practical time-steps are

used. Using this new method, computations have been made for an oscillating wing at a transonic Mach

number. The computed results confirm that the present upwind scheme captures the shock motion better

than the central-difference scheme based on the Beam-Warming algorithm. The new upwind option of the

code allows larger time-steps and thus is more efficient, even though it requires slightly more computational time per time-step than the central-difference option.

17. Key Words (Suggested by Author(s))

Unsteady flows

Upwind algorithm

Navier-Stokes computations

Moving grids

18. Distribution Statement

Unclassified-Unlimited

Subject Category - 02

19. Security Classif. (of this report)

Unclassified

20. Security Classif. (of this page)

Unclassified

21. No. of Pages

22

22. Price

A02

NASA FORM 1626 OCT86

For sale by the National Technical Information Service, Springfield, Virginia 22161