A Reproduced Copy - NASA Technical Reports Server (NTRS)

2 downloads 0 Views 2MB Size Report
Sep 1, 1981 - (NISI-CB-173254) TUDEE-DIMENSIONAL. EULER EQUATION SOLUTIONS USING FL~X. SPLITTING1Mississippi State Unlv.) He A04/MF AO ...
NASA-CR-173254 19840008789

A Reproduced Copy OF

Reproduced for NASA

by the

I'IASA

Scientific and Technical Information facility 111111111111111111111111111111111111111111111

NF01452

I'

SEP 1 1981 LANGLEY 1\·~cSfAf?CH GEN I ER 1.18RArlY, NASA

FFNo 672 AU~J 65

tlAJYIf':ION,

VlltGII':II"

, ,I

.-'

~

.

'. t

THREE-DIHENSIONAL UNSTEADY EULER EQUATION SOLUTIONS USING FLUX VECTOR SPLITTING

(NISI-CB-173254)

TUDEE-DIMENSIONAL

EULER EQUATION SOLUTIONS USING

FL~X

SPLITTING1Mississippi State Unlv.) He A04/MF AO 1

Ui~STIHWY

N84- J 605"1

VECTOR

58 p cseL t2A

Unclns

00555

.. ,..

,

.

.;

.,/ uect::mber 1983

.

'

,

~:.

....

'

-. ,.,,"

.

., ' .. "

.

r,.' ,

.

, iii'"

FOREWARD Thls research was sponsored by the NASA Langley Resenn:h Center under G~ant ~A~A NAG-1-226 with Dr. Frank Thames as TecJmifal Honitor,

and the United States Air Force, Eglin Air Force Base. Florida under Grant F08635-82-K-0409 wi th Dr. Larry Lij m.rski as Tt'chnical:. Officer. The support of NASA and the Air Force, as "Jell as tedmical. interactions with Drs. Thames and Lije\"s1d. in gratefully

"

acknotvledged~

. \

".,

," .

.... '

,./

" " I

,.\

"

..

t '.: . ~' ,

,

"

, ....

ABSTRACT A method for numerically,.soJving the three-d:f.mensional unsteady Euler equations using flu;cvt!"etuj: split.ti.ng is developed.

The equations

are cast in curvilinear coordinates and a finite volume discretization is used.

ru.

explici.t upwlnd second-order predictor-corrector scheme

:L8 used to solve the discreti:r.cd equations.

GFL

num~er

The scheme is st[tble· for a

of 2 and local time stepping is used to accelerate convergence

..

,,",,;

for steady-state problems.

Characteristic variabl~'boundary conditions

lIre developed and used in the far field and at surfaces. dissipation terms are included in the scheme.

No additional

Numerical results are

compared with results from an existing three-dimensional1!:uler code

and experimental data,

.....;

.,

.\

I.

INTRODUCTION

In support of NASA's interest in the use of prop-fans as a propulsion device, acomput'll1,t'l:aliitl mert,Jlod was recently developed for numerically solving the flowfield about a swept. tapered, supercriti.cal '''ing with a propeller producing thrust and

s\~irl

(Ref. 1).

The equations solved

were essentially Euler equations with body force terms included to simulate the propeller.

The Euler solver used was an extension of the

method of Jameson. et. aL (Refs. 2 and 3) xcferred to as FL057.

Although

good agreement was obtained with exped.mental data (Ref. l). some difficulty in convergence, particularly the energy equation, was encountered using the FL057 central-difference scheme.

Recently. convergence

difficulties Here also reported by Swafford (Ref. 4) in solving a sct of hyperbolic equations \'Jith source terms (ltlh:tch can be consj.dered similar to the force terms included in the Euler equations in Ref. 1) using a similar central-difference scheme. the convergence problems in Rcf. 4.

An upwind,. scheme eU.m:l.nated

Horeover, any additional smoothing

added to the upwind scheme was found to be detrimentaJ. with regard to convergence •. (Additional smoot.hing was found to ahmrys be necessary in usi,ng the central-difference scheme in Ref. 4,)

Tlt~e):efore.

it

seemed appr9priate to consi.der an upwind scheme for Sl:!lv:!.ng the thrcedimensional Euler equations. The flux-vector-split form of the equations used is developed in the following section.

The motivation and bacl,ground' of using splitting

is given in the literature {see, for example, Steger and \varming (Ref. 5» and is not repeated here. given in the literature

Horeover. many of the ma.tr:kes needed are (see Refs. 6 and 7); hot-lever. all equations

~

.

"

and matrices needed are developed and included in Section II fo·r clarity and completeness.

Formulation of

t~e

equations for numerical solution

a.nd the algorithm used to solve the discretized equations are discussed in Section In,' 'i1!f~~;; 'Ciu.:ract"~'· ~ siic variable boundary conditions used :l.n the far field and at surfaces are developed in Section IV.

NUmerical

results, including compa:dsons with FL057 solutions and experimental data, are presented and discussed in Section V.

"r'

.,

.,

II.

SPLITTING

The conservation law vector form of the

I:U))~$.!·Jr

equations in

Cartesian coordina:t.#asx .• y. and z are

(2.1) l~here

-

q '" [p, pu, pv, P\'1, el

T

f :: Cpu, pu 2 + p, puv, puw, u(e -7i' p),l g '" [pv, puv, pv

2

h '" [pt., • puw, pvw • P'd

P

=

.

+ p,

(y - 1) [e - JliP(u

pvw, v(e +. p}}

2

T T

+ p. wee +. p.~, JT

222'

+

+ w )1

v

.

'..

Using curvilinear coordinates defined as E:; ::

..

E:;(x,y,z)

.1

n '" n(x,y.z) I; '" l;(x,y,z) T '" t

,. ""

i t is straightforward to transform Eq. (2.1) t.'o

(2.2) "

,,,here

.

. Q" J [p, pu, P\', pw,

el

T

"

,

F" J[pU, puU

+

G = J[pV. puV

+ nxp.

H" J[pW, puW

+ C;)~p, pvW + C;yP, pwW ...

l;xP , pvU+E:yP, pwU pvv

+ nl'

pwv

+

~zP. U(e

+ nzp.

+

p)]

T T'

V(e

'+

p)]

l;zP, W(e

+

p)}T

, J = xE;(ynzr; - znYI;) -YE;(~nzz; - znxr.)

+ z~(xnYr.

- Ynx~)

z y ) n r;

nx 1;

E;

x

Y

.. J .. J

'" J

ny '" J /;

Y

E;,

z

n

, Z

'" J

-1

(zf;Yl; - y :r. )

t:

-1

(z x

n /;

-1

-1

/;z .. J

(Zt;X

n

-1 (Y

-1

- x z ) n r.

- x Z ) t; n

- Yn,cr;)

t;'\ -

(xE;,Y

n

x,y 1;)

- Y x )

t;, n

U .. t; u x

+ t: v + t: zv.' Y

V'" nx u

+ nv + n

W ..

ORiGiNl~t. rt'~c!S r~

OF POOR QUALITY

(Xt;Z I;; - z t; x c; )

= J -1 (x nYr.• = J

zE;Y ) n

(Y E;zn

-1

I;;

'

~

Y

r; x u + r; v + y

.'., ~

\oJ

'(,zw

The strong conservation law of the Euler equati.ons i.n curvilinear coordinates (Eq. (2.2»

can be wr.:ltten in the quasilinear form

.. ,

.

"

j.

, -I' . " ..

(2.3)

"

. ~ ~ .' I

,

"',.

i \ .. "

.

where the matrices At B, and C are given by

'.'.3k '" k xu +

It v Y

+ k zW

'"

Ok

(2.13)

A5 ., 0 I, It. !~

-

C

I17k I

/

ORlGfN.4L P,l~Cf! ffJ

OF POOI~ QUALiTY 1

2

345

,,,here A , A , A , A , and Ak are the eigenvalues of a for k = E;, b for k k k k k = 11. and c for k ::: r;.

Correspondingly,

O]T [k;x' 0, k z • -ky ,

Xl

O)T [It • ·-k z' O. kx ' y

x2

::

x3

= [It , k ,

z

1

x .. 4 x5

the eigenvec tors are

y

[l~

e

Ii

-kx • O. ('IT

(2.14)

-

T Kx ' ky • kz • pc]

1 = -12

P pc] 'I." ['(! • -kx' -ky • -k z'

_

k

where kx =

x

TViZr =

k

x

k

ky

=~T

(2.15)

Sufficient equations have now been d,~veloped to carry out the splitt:l.ng.

The integral conservation

18\11

form of the Euler equations

will be formulated and discretized for numerical solution :i.n the following section.

This formulation requires evaluation of the

vectors F. G, and H at faces of the f:Lni.te volumes.

By splitting

the vectors F, G. and Ii into thE! sum of separate vectors. each having an eigenvalue as a coefficient. the evaluation of each separate vector by extrapolating from the .!lppropdate direction indicated by the sign of its eigenvalue can be carried out.

OH!GiN;\~. i:;:.,~l~ ~r~

OF POOl1 QUt>.U'f'l The ve..:tors F, G, and H are homogeneous func.t:h:vns

(I.e

degree one

in Q; therefore, by E41er's Theorem (Ref. 9) ;2.H)

K = KQ

where K == F and K == A for k K = Hand K '" C for k

=

.=

E;. K

=G

and K

'C

B for It == TJ. and

The matrices K can be· t;rrlt.ten

1;.

(2.17)

where Ak is the diagonal w.atrb: whose diagonal elem.€nt·"I are the eigenvalues of

K

(cr, also, K) given by Eq. (2.13)

The matri.ces K given

by Eq. (2.12) can be Hrittcn

K ==

where .t·he columns of P

k

respective eigenvalues.

P II p-1

are the

(2.18)

k k k

(;.~gel1\

eCLot!:; of

K

corresponding to the

From Eqs. (2.8) -

K '" MKM

-1

(2.19)

Using Eq. (2.18) in Eq. (2.19) (2.20)

Th,en from Eqs. (2.17) and (2.20)

(2.21a) ··1 -1

Pk M

The matrices M and M- 1 are given by Eqs. (2.9) and (2.10). matrix P

'\

k

is

(2.7.1b)

The

. ,

.;

.' , \\

.....:-

, ,

-~ ~,

\

:~ ORIGINAL PI\[~2 ;r:J POOR QUAtnv

OF

fI i \;.

k )C

1,

y

-k Z

k

k

0

-k x

k

0

Y

Ii

'-..... ~'" i

1' • ...,

\'L

P

I

,

,,,/~' /

k

..

z

k _'L

If

"I.

,.>

k

-k

Y

x

Z

If ac

0

0

0

~

k -2S.

0

! !

a

z

2

.

__k!S

Ii __It.l:'

(2.22)

.,If k

--1'2 &

ac

2

where

... ]

Rather than inVel"t P cUrcctly. the mntrix P ' can bu dt\t(~rl\linerl casily k k from the left: eigenvectors of

of P

-1 k

left (~igmlV(lctors of

K

are the

rO\,6

• Le.

k

k

o

.

-kx

0

Y

k _l£.

-..Y-_

k

k

If

fi

Ii

z

k

-

-ky

It

k

z

z

0

-k

Y

o

k

0

x

..

k

It'

'."lli;!

Ko

x

----

/2

x

(2.23)

...' "

k

k

12



.. _'1-

. ~

z

z

I"~

.-:. 4

:,

"., .

.!/. \

ll~ ~

,

'.

OmG~NAl ~AGI! V~ Of' pom~ Q"AU1'Y

\-Ihere P ""

1 .-=---

Ii

~;i~

The mat.rices Tk and T;l can now be determined using Eqs. (2.21). These matrices are k

TIt ""

k It

k Y

uk x

uk

vI-k

0

0,

fJ

0

0

;.,3

0

0

0

0

0

>-,.t.

0

0

0

or

,1:.-

. 2

0

II

It

..

It

I~·

(2.26)

Q

~~•. i

. •.

;

,can be writt:en II

k

c:

,,1 I k

'.l

;It ,r ).4 I., 1,2,3 + k 1,+ + 'fu. ,)..

(2.27)

.

1iJhere 11 2 3 is a matrix \vhich has unity lAS t:he. £:h:st three diagonal

,

elCluents and Clll other elc:mcnt$ zero, 14 has !)A1ity as the fourth diagonal .~lement

and all other ele;nents zero, and IS i\13:'", unity as the fifth

diagonal element and all other elements zerO'.

lirO!t.e that only three

terms are needed on. the right hand side of Eq: •. (2 •.27) instead of five because the first three eigenvalues are I'he same '/;;' and ASr,; • , (. ':I eigenvalues are interrogated. and Similarly for Qi . ~n 1 Qi j f k except A •

, . k' -!2.J.

~n

i.

.. OR/(~/N I..\./

/.,,:.... .. " .. rr)f'.t.:d .... :.-3

. OF POOR QUALITY

Corrector:

(3.10) " "

'::n+l where the Qi+~,j,k' etc, variables are determined: in the same way ~n

as the Qi+'ii,j ,k' etc, varia.bles as described belo\Yi Eq. (3.9

n Qi.l,k' etc, are used in place of Qi,j,k' etc.

)

except

"n

The Qi+lti.j.k' etc, are

determined by

n Ql-1,j,k

,1 4 \5 if the corr:iEsponding A~ • A~, or At; eigenvalue evaluated at i+~,j.k

is :;. 0 or "n n Qi +~, I j , k" 2Qi+l

1

4

if the conr.e:sponding

A~,

eigenvalue evaluated

~t i+~,j,k

5

'.)

I.E;' or I.E;

is < 0

"n

and sim:l.1arly for Q,

I.

l.-'~.

j

,

k'

...n

Again, the same is oY:r.ne for Qi j+l_ k ,

-"'2,

n and Qi,j ,k±.lli by interrogaUng the appropriate nor. t eigenvalues. Although the algorithm given by Eqs. (3.9) and< (3.10) is an extensi.on to thr.ee dimcnsiona of that used in Ref. 13 for two dimensions, the interrogation of ei.genvalues differs significD.l'1.t1y •

Three methods

were tried in Ref. 13 to handle the computation of flow variables at cell faces Hhen e:lgenvalues were of different sign on either- side of

..

~

,/

ORIG1rJP.L FJ\G;;." ~'.; OF POOR QUALn:y the cell face.

The approach taken in the present

sdH~me

is to compute

e:Lgenvalues at cell faces rather than at cell centers as illustrated in F:lg. 2 on a z; ..

constlmt':~lan"',

changing sign is thus eliminated.

The difficulty ass.-r"""'_·""".~~""""'~""'''_''' __~-''''·''_'~'' _ _ _''''''-'''---_''--_'·-

,

.•~-----,.,---.,.... ...

G, I.

Ci~IGli~f"tl P!:.C!-~ tS OF POOR QUALITY

.

a(p-l q) h.,o

aT

a(p-l -I- It

q)

k,o

ak

k

+ s

k,m

.,

0 .

(4.4)

Defining the characteristic vector, W , as k ,I'

.

(4.5)

Eq, (4.4) becomes

(4.6)

The elements of the characteristic vecto.r, H , are -::.alled characteristi!! k variables, and are denoted hy wk,:i.'

., are determLed using Eq. (.f. 5). k ,~

The charac-teristic variables, w The vector q is

q = J [ p, -1 and the matrix PI

( ,0

U,

T v, w, p]

is given by Eqs, (2.23) ,.here the variables p and

c (speed of sound) in Eq. (2.23) are denoted rand c o

a reference

~onditil)n.

Using q and P

-1

k,o

c

to indicate

in Eq. (4.5), the elements

of the characteristic vector (4.7) "

.' are

.

~,

wk,l .,

J

IVkl

[kx(p -

-1-) c

+ kz v

k \~'l y

°

J

Il1k/

(I•. 8b)

c0

··,. ·..... ., ,,,

"

',.'

wk , 3

=_J_. [k (p

IVkl

z

- _1.'_) + k u 2 y c

°

.

"

wk 2 == - - - [ky (p - .J?i) - k z u + l~x 'v]

,

\

(4, Ra )

- kx v)

(4.8e).

'.

,,'

'.

"

,

.

1 . -J - -

(4. 3d ),

.. ---IVkl 12

(4.8e)

IVkl

Ii 1

J

The chat'llcteristic variables correspond in o'l"der to the cis111 in

3a with the flow in '.i'! dir0ction of increaslng

(~omputntl.onal

sign

l~ig.

cO(I:!.'dinHte

k, the first four eigenvalues ure posi.tivo nnd the fifth 113 ncgnth,"Iil.

,.

.,

'"

~

For the subsonic inflow case 5hO\m in Fig. 3b ,·lith the flow in

t.ht~

directioll of decreasing computational coordinate k, the first thl'.'tle nnd fifth eigenvalues are negative. and the fourth eigenvalue is positivu.

For a totally general three··dimensional code eithfell' situation in rlg. 3 Each possiblity in ! _.

1'ho subscript a refers to approaching the

the situatIon in Fig. 3b.

lind

r·· ',..

ilK€< defined by Eqs. (2.15).

Equat:ions (4. 1 2a),

(4. 12b) , (4.12

P.e

could be the ambient static pressure,

, ..

.,..''

·t· . ~ '1

~.: '.

,-' 'I.",'

\

ror a boundary across which there is no

flOlli

\

the first three

eigenvillues given by Eq,AA':'ff) are zero. the fouz:th is positive. and 'the f,ifth is negative.

One' condition must, therefore,. be specified.

'rhe condition specified is that there is no flow across the boundary. The following relations among the characteristic variables are used to determine the boundary conditions

-1)

- k w] + kv z y b

[k (p - ..Q_)

[~(p -

c0

y

[k

._----.

c

(p -

g'

2

=

[k (p •. .....E.) +k.v - k w) z x 2 y r

(4.158)

k u + kxw]b z

[k (p •. _2...) - k. u y 2 z

+ k xw] r

(4.15b)

k u - k v]

[k (p - L) +. 11 u - k v] z Y x r c2 o

(4.15c)

c0

0

.l?.) 2

+

Y

x

b

Co

/

.....

c0

'.

[kxu + ky v + kzw]b = 0

,'

kw) ) z

b

..

[.p IVk I. -+ Pc

o

0

~Jf.",~U '"

... +

k v Y

+ It w») z

r

(4.15c)

.,"f

.. .(

The subscript r refers to a reference vahlc. which is selected as the center of the first cell from the boundary.

The'minus and plus signs

-., ~.

:In Eq. (4. 15e)

corr~spond

to the locntion of the point r.

If the point

, "

lC

is in the positive k direction from the boundary

then the minus

:;1gn is used in Eq. (l1.15e). and i f it is in the In:bt.us direction then

"

.'

.1

the plus sign is used.

""

.

, \

Finite volume codes only require the pressure at an impermeable boundary and consequently Eqs. (4.15a). (4.l5b), and (4.l5c) are not needed.

However, to

fad:;i.~;',;t:;lW.,"','Lhe

handling of

point~'

near boundaries I

and aid in code vectori.zllItion, phantom points are used in the present' version of the code.

The use of phantom points r('!quires information of

variables other than pressure to ensure, for eJtal1lple, an impermeable boundary.

.

l/!Cl'O

flow across

Such information can be obtained from Eqs.

(4.15) • '~

Equations (4.lSd) and (4.1S e ) can be solved for P • b

Equations

(4.15a) - (4.1Sd) can then be solved for the remaining four variables. The solution of Eqs. (4.15), is ~

P b '" P r :j:.' p oc 0 (kx u r

P b '" P r

ub

.. ur

;.

Pb - P

c

r

(4.J.68)

(4.l6b)

2 0

- k (k u

x

+ kv + k Zwr ) y r

x r

+kv + k ".. ) Z r y r

(4 .16c)

.

4,

''..

"

v = v

,- k (k u

+ k Z wr ) y x r +kv y 1'-

(4.16d)

+ k Zt"r ) wb '" wr .- k Z (kxu r +kv y r

(4.16e)

b

r

-. .

. '

"

where the point r 1s the center of the first cell fr,?1n the boundary and the minus sign in Eq. :I '

(4~16~

is used if r is in the positive k

direction from the boundary, and the plus sign is used if r is in the negative k direction from the boundary.

,......

"

. i , I '.

/."

I

I

......

4.3

Phantom Points Phantom points are denoted by the subscript p.

obtained from the relations

/

/ ,/

The points are

",.

:,.'"

(4.17a)

p

U

p

p

:: 2p

::

2u

b

b

- p

- u

(4.17b)

in

(4.17c)

in

(4.l7d)

(4.17e)

where th(,: subscript in refers to the center of the first cell inside the computatj.onal domain and can be any of the points a, in this section.

./ .

For

.t, 'or r used

the phantom points for an impermeable

~xample.

. surface are

+ 2p oc 0 (Itxu r . + Iey v r ~

p

p

Pp

\

\

u •

.

Pr

2(P

. ,.\ .

=:

p

.. p

·r

., u

r

+

k2\~r)

- P )

b +--r2

c

(4.18a) (4.1Ub)

"

.,

"

0

.

x r + kv y r + k:e.\-lr>

(4.18c)

- 2k (k u +kv + k z wr ) y x r y r

(4.1Bd)

2k (k u

x

..

I

v

/

p

.. v

r

./

w .. w p r

/

"

:>'~-~' .I

z

x r

+kv

y r

+kw) z r

.. ~

(4.18e)

The velocity vector components are the same as those used by Jacocks and Kneile

.I

2k

(k u

/:

(R~!f.

14).

,',

..1~

' ,

v (t~

' ....... r.-:.

E t" -< '.:.~

..

ltil~~I~~ ~ I·

I

'It 0 •• :"

'~

-cp

~ ".~

. \

1~ ........

--,.... "'-.~......

,"

,".

---."~,

.....

/_

-

.....

.

'-,

..

'~

\~:

\, .\....... ,

\-

'.\," ..

\ \;

'"

.\

• F ... • •

..

'-,

\,

...~,,:-

~-=-

-.

..........

,.-

.:.!

~~'--~.; --. ...

'>-.

---.~.

.......

!

.\

•• _k

J"'-.. .', • . ,'\. '~.:'

~

".

~

....

'

~

",

-

'-.

"0

"

_

~.~~--~------

.

.....

!!. 1.4

r

ONER4 ~16 WING' M", = 0.84

.1~2 ~ A I.D

c.~

~~~~ I

!

HI

\

I

=

"'1

11

=

0.45

~

---:.

» \~-'e.-~/.

-cp

0

v...,

\

c.~ i 11 13.4

nh

~

-...:.

"\

H-Il

ti

.



:.:'

.

.

:

::::

"""" x LO

1.O

r-i

8

tJ

I='

r-I

/'

,/

;

/'

. j'

'-.J"

'-'

I'..

-' .:::r;

< 1.O

O'l

/

LO r-f

r-i

>< 1.O

;.

II

II ,~

'..0

-

.-...

II

r,

1.O

1.O

X

~

rl

s::

roof

1.O

rl

.-f X

X

,

~O

l.O

cn

cn

'-'

"-"

, I'·

'-. N

• ;.,1

-

t.:I



.

100

,

",

\

'"

,

),

'"

0.1

o

o

CD(TOTAL.)

I

\/ \'

I;j

\ \



ONEPJ\ N6 WING ~l >: 0.84 01 ~.:.: 3.06 ---- fl057 (48xS)