A compressible vortex method with application to the interaction of an oblique shock wave with a boundary layer

Applied Numerical North-Holland

Mathematics

257

8 (1991) 257-273

A compressible vortex method with application to the interaction of an oblique shock wave with a boundary layer Gary

A. Sod

Tulane University, Department

of Mathematics,

New Orleans, LA 70118, USA

Abstract Sod, G.A., A compressible vortex method with application to the interaction boundary layer, Applied Numerical Mathematics 8 (1991) 257-273.

of an oblique

shock wave with a

A numerical method is presented that can resolve a compressible flow field at high Reynolds number. It is used to study the interaction of an oblique shock wave with a boundary layer. These results are shown to be in excellent agreement with experiments.

1. Introduction supercomputer capabilities have made it possible to study fluid dynamics not just through laboratory experiments, but also by numerical experiments. A careful numerical computation can often provide new insight into difficult physical problems. With these numerical experiments we have a greater control over physical effects, such as the presence (and amount) or absence of viscocity and the choice of the boundary geometry. Also, we have greater control over the data, that is, we can extract a great amount of data from the numerical experiments (in principle the values of the dependent variable are obtainable at all points on space and time), far more than is possible from any laboratory experiment. The compressible form of the Navier-Stokes equations is a coupled nonlinear mixed hyperbolic-parabolic system of partial differential equations. The viscous terms which cause the system to be parabolic are of the order of the reciprocal of the Reynolds number in magnitude. At high Reynolds number, the system is almost everywhere hyperbolic; the viscous terms are negligible except in their layers. Within these layers, however, the viscous terms are significant and control the important phenomenon of boundary layer separation. We present a numerical method that can resolve a compressible flow field at high Reynolds number in two space dimensions. The method makes a judicious combination of a hybrid random choice method and a vortex method. We present results of the interaction of an oblique shock wave with a boundary layer. If the shock is sufficiently strong the boundary layer will separate. This provides a severe test of the method. The numerical results are in excellent agreement with experiments. Current

016%9274/91/$03.50

0 1991 - Elsevier

Science Publishers

B.V. All rights reserved

258

G.A. Sod / A compressible

vortex method

2. Governing equations Consider

the equations a,U+ a,F(U)

of gas dynamics + a,G(U)

in a domain

= p AH(U)

in

Q with boundary

a&I?.

P-1)

$2

where

lJ=

I

‘P\ PU p”

F(U)

5

\ e

=

puv

\

)

p(e+d, /

G(U)

PU Pu+P

=

PU PUU

/o\

\

U

puv*+p

w4=

’

v

>

0,

,h+P),

is the density, u is the velocity in the x-direction, v is the velocity in the y-direction, p is the pressure, e = PE + $p( u* + v*) is the total energy per unit volume, E = p/( y - 1)p is the internal energy per unit mass, y is the ratio of specific heats (a constant > l), and p is the shear viscosity associated with the flow. The Hodge decomposition theorem [4] is used to decompose the velocity field u = (u, v) into a divergence-free part, U, (V . t.$ = 0, i.e. incompressible part), and a curl-free part, v, (V X u, = 0, i.e. irrotational part), that is u = ud + v,. Both U, and u,are required to separately satisfy the boundary conditions p

ud.nl,,=O

(2.2)

and 0,. n I aa = 0,

(2.3)

where n is a unit normal vector satisfy the boundary condition

to the boundary

afi. However,

only 0 = V, + 9 is required

?J*sIaa=O,

to

(2.4)

where s is a unit vector tangent to the boundary a0. Embed the region fi in a grid of spacing h in both the x- and y-directions. Denote the numerical domain by fi2, and the numerical boundary by a&Z,,. Divide time into intervals of length k. Let UjIf approximate the solution U(x, y, t) = U( ih, jh, nk) to (1.1). Let UCdenote U where 2) is replaced with 0,. At high Reynolds numbers, the viscous terms are negligible except in thin layers. Exploiting this, V, is determined by considering the second and third equations in (2.1) which, using the first equation in (2.1), may be written in the form p(&u+

(IJ.V)U)

=p

Av-

vp

in 9.

(2.5)

Define the vorticity by 6 = v X v = V X v, and [ = 6. k. Taking scalar vorticity transport form of (2.5) a,<+(v.v)[=Re-‘A<+

-$(vpXVp).k-

the curl of (2.5) gives rise to the

y(vpXAl;).k

inO,

(24

259

G.A. Sod / A compressible vortex method

where Re = @L/p is the Reynolds number with U representing a characteristic velocity and L representing a characteristic length. We consider the equation (2.6) along with the condition V.u,=O and boundary

in52

(2.7)

conditions

v,.n=O

(2.8a)

onaD

and v-s=0

(2.8b)

ona0.

Away from these thin layers, we consider a,U+ a,F(U)

+ a,,G(U)

the inviscid

equations

of gas dynamics,

in fi.

= 0

(2.9)

The numerical method for solving (2.6))(2.8) is a grid-free method discussed in Section 4. The numerical method for solving (2.9) is a grid-based method, discussed in Section 3. The coupling of these methods is discussed in Section 5.

3. Random choice method The random choice method was introduced by Glimm [16] for the construction of solutions of systems of nonlinear hyperbolic conservation laws. It was developed for hydrodynamics by Chorin [6,7] and further developed by Colella [lo], and Sod [26-291. Consider the hyperbolic conservation law

a,u+ a,F(u) = 0 along with the step function

(3.1) initial condition x < 0, x > 0.

Such an initial-value form U(x, Suppose function

problem

t) =R(x/t,

is called a Riemann u,,

u:,

XE

[i-

which has a similarity

solution

of the

&)*

that u,‘, the approximation of x,

u”(x) =

problem,

:h,

to U(ih,

(i+

nk),

is given.

Consider

the piecewise

+)h].

To advance the solution from time c = nk to t = (n + l)k, consider the initial-value given by (3.1) and (3.2). On each interval, [ih, (i + l)h], this initial-value problem Riemann problem. If the CFL condition

constant (3.2) problems defines a

where pj(U), j = 1, _. . , p, are the eigenvalues of A(U) = aF/XJ, is satisfied, the waves generated by the individual Riemann problems will not interact. In this case, these solutions to the

260

G.A. Sod / A compressible vortex method

local Riemann problems for nkdt<(n+l)k, U”(x,

can be combined

t) =R i

by superposition

x - (i + $)h 3 u:, ui”+l t_nk

1

into a single exact solution

V(x,

ihgx<(i+l)h.

3

t)

(3.3)

point value of a locally In the random choice method, to obtain ,:+I by a representation exact solution, this representation point value is obtained by sampling (3.3). Let 5, denote a equidistributed random (or quasi-random) variable in the interval ( - i, i). Sampling the locally exact solution (3.3) uin+l = U(i(i R =

+ &)h,

(n + l)k)

(5A)h i

R I

i

(6.i:)h k

> u:-1, u7 , 1

E, (0,

9 u;, u,‘+, , 1

E,,>O.

For details of the solution of the Riemann problem see [6,10,26,29]. The two-dimensional hyperbolic conservation law (2.7) is solved by operator splitting which reduces the problem to two one-dimensional hyperbolic conservation laws (see [26]). Although each of the one-dimensional fractional steps models the resulting one-dimensional gas dynamics well, the combined steps do not model the two-dimensional gas dynamics well. In fact there is an 0( 1) error. To remedy this, near a region where large gradients exist, some dissipation is required (see [lo]), which is due to Godunov [16]. This is accomplished by using a diffusive finite difference method. Godunov’s method follows the construction of the random choice method through the solution of the Riemann problem and the construction of V(x, t) in (3.3). To obtain u:+i, Godunov projected V(x, t) onto the space of piecewise constant functions by n+, _ 1 u, - h

(i+i/W J (i-l/2)/I

V(x,

(a

+

1)k)

dx.

Since ZP(x, t) represents the exact solution to (3.1), integrating rectangle [(i - i)h, (i + ih)] X [nk, (n + l)k] gives

J

(i+1’2%‘(~,

(n + 1)k) dx -

(r-1/2)h

+

J

(n+‘)kF(Ue((i

+

:)h,

i;+‘)*P(U’((i&:)h,

J(i-l/2)/l

t)) dt -

t)) dt=F(R(O,

:)h,

(“+‘)kp(V-((is nk and

u:)

nk) dx

V(x,

U”((i

-

:)h,

u;, uj’J)k,

so that (3.4) becomes .;+I

= u: - X(I;(R(O,

u:, u:+J

law over the

(i+1/2)h

nk

Using (3.3), U”((i - t)h, t) = R(0, u:_~, independent of t. Substituting into (3.4),

the conservation

- F(R(O,

u:-1,

u:)).

t)) dt = 0.

t) = R(0,

ur, u:,,),

(3.4) which

are

G.A. Sod / A compressible vortex method

261

In the situation where (3.1) represents the equation of gas dynamics, the pressure is used to construct a switch to determine when the Godunov method is used. In [lo], Colella presented one type of switch. We present a slightly different pressure switch. For each grid point i, define, for a scalar grid function q:,

then

where P,?_,,~ comes from the pressure component of R(0, ~:_i, u;) and P:_,,~ comes from the pressure component of R(0, u:, unr+l), the pressure in the region separating the two sonic waves in the Riemann problem. The switch to the Godunov method at the grid point i where Q, > E, where E is a parameter that is a measure of the weakest sonic wave that is treated as a discontinuity. Note if e is chosen too small, then the Godunov method is used unnecessarily which results in smoothing of the wave. While if E is chosen too large, then the failure of the switch to detect a pressure jump results in noise.

4. Vortex method A grid-free method, the vortex method, introduced by Chorin [5,9] is used to solve (2.6). The method numerically simulates the creation and development of vorticity through the introduction and appropriate interaction of discrete elements of vorticity. Equation (2.6) is solved in a series of three steps. Consider first the inviscid limit obtained as Re -+ cc in a domain with no boundaries and where p = const., a,[+

(0.0)(=0

which may be written

in 3,

(4.1)

in the form

a,<+ ((2)d+2),).V).$=0 A stream function u= Then

[=V

in L?.

Ic/ is introduced

-a,,+ Xv=V

and Xv,=V

a$.~= -.$

by

v=aX$. x(-a,~J,d,p!~) givesriseto

ina.

(4.2)

Assume at time t = nk there are N point vortices in L?. Then the jth point vortex has center xi and strength <,. The vorticity field at point x due to the N point vortices is 6(x)

= :

5,8(x

- x,),

J=l

where 6(x)

denotes

the Dirac delta function.

e

262

G.A. Sod / A compressible vortex method

The solution of A$ = - tjS( x - x,) gives rise to the stream function jth point vortex,

qj( X) associated

with the

-6 +,b) = ~lnII~-qI. This stream function is singular at x = x,. Thus the velocity field induced by the jth point vortex becomes unbounded as x -+ x,. To remedy this, Chorin [5] chose a typical eddy with vorticity Ej spread over a positive area. A vortex “blob” at a point x = 0 is defined by

where u is a parameter to be chosen. Thus Es is a smooth field induced by the N vortex blobs becomes

approximation

to S(x).

The vorticity

N 6$(X)

C

=

tji56(x-xJ)e

j=l

,$a generates

its own stream function

W,(x)

+Ls(x), obtained

by solving (as in the point vortex case)

= -6&x>.

Solving,

&ll~

YM-4 = i

&lnllxII3

IIXII
IIXII 20.

The velocity field induced by q6( x ) is Us point x induced by the N vortex blobs is

= (- a_,$,( x), a,+,(x)).

Thus the velocity

field at a

Equation (4.1) states that the vortex blobs are advected by this velocity field plus the velocity field 0,. Given the position of the ith vortex blob at time t = nk, its position at time t = (n + 1) k is obtained by Euler’s method ,;+I

= x,? + k ( us”,+ o,, ) ,

(4.4

where

j=l j#i

is the velocity at x,! due to the remaining N velocity at x,?. The determination of & will be This method is O(k). Sethian [22] has used boundary is now introduced, the velocity field

1 vortex blobs, and oci denotes the free-stream made in Section 5. Heun’s method to obtain 0(k2) accuracy. If a must be such that there is no flux across the

263

G.A. Sod / A compressible vortex method

boundary, that is, V, . n = 0 on i3fi. This potential function I#Isatisfying A$=0

inQ

V+*n

= -v,(x)

.n

boundary

condition

can be satisfied

by finding

a

on ati.

Let ur, = v+, then U, + u,, is a velocity (4.4) now becomes

field which satisfies

the boundary

condition.

Equation

The above Neumann problem may be solved in a variety of ways. Chorin [5] used a single-layer potential distributed along the boundary, satisfying an integral equation. The discretization of the distribution leads to a set of linear equations which is solved to find the potential strength at points on the boundary. Davari [ll] uses the method of images to satisfy the boundary condition for flow past a cylinder by introduction of image vortex blobs, centered inside the cylinder and thus satisfying Laplace’s equation in the flow region. Cheer [3] uses conformal mapping to consider flow past a Joukowski airfoil by mapping it onto a cylinder. Ghonlem, Oppenheim, and Chorin [14] consider flow past a backward forcing step in a channel by mapping the region onto the upper half complex plane. In this geometry, complex conjugate image vortex blobs are introduced to satisfy the boundary condition. All of these methods have the advantage that no grid is introduced, the potential velocity is only calculated at the position of the vortex blobs, where it is needed. However; the scope of the problems accessible to these methods is limited by geometric constraints. Sethian [23] uses a package of elliptic partial differential equation solvers to obtain the potential flow in a bounded region and in a region with prescribed inflow and outflow boundary conditions. These solvers use grid-based methods, so interpolation is required to obtain the solution at the vortex blob locations. However, the methods are not applicable to unbounded regions. In the application presented in Section 6, the domain in the upper half-plane ( y 2 0) both the boundaries are along this x-axis. The method of images shall be used. The normal velocity on the boundary produced by a vortex blob of strength tk and location ( xk, yk) is balanced to produce a no flux condition by adding the effect on an image vortex blob of strength -<,, and location (x,3 - yk) (see Fig. 1). To modify the vorticity field due to the effects of diffusion, we consider the diffusion equation $5 = Re-’

(4.6)

A$.

>x

“\\\\\\\\\\\\\\\

($‘-Yk)

13 l

-5,

Fig. 1. Vortex blob and its image.

264

G.A. Sod / A compressible vortex method

In one space dimension

a,<= Re-’ The solution

this reduces

to

a,‘,$.

to this one-dimensional

(gy2 exp(

((x, t) =

diffusion

equation

with initial condition

<(x, 0) = 6(x)

is

- x2/(2t/Re))

which is the probability density function associated with a Gaussion distribution with mean 0 and variance 2t/Re. The effect of diffusion on the vortex blobs from time t = nk to time t = (n + l)k is modeled by allowing each blob to undergo a random walk in both the x- and y-directions. The combined random motion of the vortex blobs will approximate the solution to (4.6). To this end, let di = (n,,, q2,) denote a random variable drawn from the above Gaussian distribution. Then for i = 1,. . . , N, ,,y+’ =I,: + k(us”; + tf;) + q;.

(4.7)

It remains to satisfy the other boundary condition (2.8b), the “no-slip” condition, first consider (v, + v,) . s = 0 on as2, where s is a unit tangent vector. We assume that vorticity within a boundary layer is primarily produced from large gradients of u in the y-direction (assuming the boundary lies along the x-axis). Vorticity is produced as a result of the moving fluid adhering to the boundary (due to the “no-slip” condition (2.4)). The vorticity diffuses from the still fluid at the boundary to the main flow in the interior. With this in mind, Chorin [8,9] introduced the vortex sheet method. Consider the Prandtl boundary layer equations in vorticity transport form, with 5 = - $u, a,,$+ (v.v)$=Re-‘aj5 v.v=O u(x,

in fi,

in52, cc> = U,(x), on

v=O where u is tangent

ati,

to ati and v is normal

u(x, t) = U,(x) +

to a52. Integration

lm,$(x, s, t)

of 5 yields

(44

ds.

Y

Since the flow is incompressible, u(n,

t) =

=

-ax(L(x,

-ax/[0’

=

-a, [

that is, v . v = 0, we obtain

2,

t) dz

u,(x)

U,(X)Y

+

Jm<(x, s, t) z

+

u,(~)~+

1

ds dz

~y~mE(~, t) ds] dz S,

lo"Jmi"'"“')~(~,

S,

t) dz ds]

0

zz

U,(X)Y +

so%t)(mh(s, X,

s,

Y>)

ds . 1

(4.9)

G.A. Sod / A compressible vortex method

A vortex sheet of intensity 5 centered at (x0, YO) of length to jump an amount 5 as Y crosses Y =Y,,, that is, u(x, Y,‘) Consider a vorticity distribution within a boundary layer strength .$‘, and center XT at time t = nk. The approximation u: = U,(x;)

+ $$, +

c

265

h causes the tangential velocity u u(x, Y;) = -E. consisting of N vortex sheets of to (4.8) yields

tjdj,

J

y; ” y,”

where dj = max(O, 1 - ) x7 - xJn I/h). The number of interactions among since for a segment S, to influence a given segment S, it must intersect R;=

((x3 Y)l x,-

:/z
the segments is small the semi-infinite strip

i/r},

when x, + ih are the ends of the segment S,. The term d, is a weight which corresponds to the percentage of the segment the strip R,. If S, does not intersect the strip R,, then d, = 0. Below the sheet S, the velocity is increased by 6;. The velocity of a segment increased by all segments of sheets which lie above Y. Use of a centered difference approximation to 8, in (4.8) yields $=

--

1 h

(J:l’u(x,Y+ +h, y) dy-

so”‘u(x;-

;h,

S, which lies in at a height

Y is

y) dyj,

and using (4.9),

where y,* = min( y,!‘, y:) The integral

and d,’ = max(O, 1 - 1x: + :h - xs 1/h).

II

-’ u x,: + +h,

J (

y) dy

0

is the flux of u across a segment of the vertical line at x,: + ih between 0 and y:. Each segment which intersects this vertical line segment affects the velocity of this line and hence the flux. If a segment Sj lies above y:, that is y,” > y:, then S, increases the velocity of the entire vertical segment from 0 to y,” and hence the flux is increased by yytj. If a segment S, lies below y: that is yJy
= x; + ku;,

Y*n+l =y;

+ ku: + n,,

(4.10a) (4.10b)

where vi is drawn from a Gaussion distribution with mean 0 and variance 2k/Re. To satisfy the “no-slip” condition 0. s = 0, that is, within the boundary layer u = 0 at y = 0, vorticity is created. At each time step compute u at discrete points a distance h along aL?. If u # 0 at any of these discrete boundary points, vortex sheets are created with total intensity equal

266

G.A. Sod / A compressible vortex method ti(xi*Yi) ,_______________ I

:$

Y;

--------I ! ,

i___________________----2 x.1

f.l

YTI x.+h 1 2

2

Fig. 2. Circulation

around

a vortex sheet.

to u at the boundary points (i.e. u(x, 0’) - u(x, O-), where u(x, O-) = 0 by the “no-slip” condition). This provides a transition from the “no-slip” condition to u(x, 0’). These newly created vortex sheets are diffused off the wall using (4.10b). To connect the interior flow calculation from the vortex method to the boundary layer calculation, the component of the velocity field tangent to the boundary (i.e. (0, + u,) . s) is taken to be the free-stream velocity U,(x) in the boundary layer. In this way the interior flow determines the creation of vorticity. Where a vortex sheet leaves the boundary layer, it is converted to a vortex blob, where the strength of the vortex sheet 5 is multiplied by h (the length of the sheet) to yield the circulation of the newly created vortex blob. To see this, take a segment with strength [, at location (xi, y,) and consider the circulation around the segments, that is, in the region [xi - +h, xi + +h] x [v,-,

y,‘l (see Fig. 2),

Lx2;[,+((x,y)

dy dx = /x.-h’2p* x,-h/2 x, +

=

- a,u dy dx Y,

h/2

(u(x, J Jx’ih’2& -

x,-h/2

=

y,‘) - U(X>Yi-1) dx

dx

x,-h/2

=&h+O(h2). In this way vorticity is transferred from the boundary layer to the interior flow. Similarly, if a vortex blob enters the boundary layer, it is converted to a vortex sheet with strength t/h, where 6 is the circulation of the vortex blob. Even in separated flow there is a thin viscous sublayer near the boundary where the boundary layer equations remain valid. This is an essential fact. For this allows us to use the vortex sheet method near the boundary even in a region of separation. To determine when a sheet becomes a blob, let L be a length chosen so that P[ 1771 > L] -c E, where E is some small number. Using Chebyshev’s theorem (see Lo&e [IS]), L will be a multiple of the standard deviation \/2k/Re. Any vortex sheet which is greater than a distance of 3L from the boundary will be converted to a vortex blob (see [26]). In computing the velocity field there are three possible interactions. A blob-blob interaction is treated using the vortex blob method, a sheet-sheet interaction is treated using the vortex sheet method, and a blob-sheet interaction is treated as a sheet-sheet interaction and as such uses the vortex sheet method. Cheer [3] did not consider blob-sheet interactions. The latter interaction allows for a consistent transition when a blob is near the edge of the boundary layer or the

G.A. Sod / A compressible vortex method

I I

(Xb,!db)

I I

267

blob treated

I ___________________-___

as sheet

j---

edge

of visious

sublayer

I I

-

I

’

(Xs,Ys)

I

I I I

I I

\\\\\\\\\\\\\\\\\\\~~~

Fig. 3. Effect of a vortex sheet on a vortex blob.

viscous sublayer. In a blob-sheet interaction, the tangential component of velocity of a blob at (x,, yb) is not affected by a sheet at (xs, y,). This can be seen by the discretization of equation (4.8) where yb > y, (see Fig. 3). However, the normal component of velocity of a blob will be affected by a sheet that intersects the semi-infinite strip S, = {(x, y) 1xb - h < x < xb + h} (using discretization of equation (4.9)). This provides a continuous transition when yb approaches the edge of the viscous sublayer. Similarly, in a blob-sheet interaction, the tangential component of velocity of a sheet at ( xs, y,) will be affected by a blob at (x,, yb) if the blob treated as a sheet intersects the semi-infinite strip S, = {(x, y) 1x, - :h -c x < x, + ih}, since yb >ys (see Fig. 4). Also the normal component of velocity of a sheet will be affected by a blob if the blob treated as a sheet intersects the semi-infinite strip S,’ = {(x, y) 1x, - h < x < x, + h }. The cut-off parameter u, used in defining the vortex blob, is taken to be h/q (see Cheer [3]). Also a long-range cut-off, S,, is used with the vortex blobs, so that if the distance between two blobs exceeds S, in (4.3) that particular blob interaction is neglected. In an exterior flow problem, when a blob travels downstream beyond the edge of the computational domain it cannot be deleted. What we have chosen to do is retain the blob but keeping its velocity due to other blobs fixed for all time at the value just before crossing the edge of the computational domain. These “old” blobs are allowed to undergo a random walk.

[As;-1 I I I

blob treated as sheet 1 affecting both tangental and normal ’1 components of ----I--_+__ velocity of sheet

1

I

;ws,+1

I I (X&b)

I I

______

_I-__~

1 $,Y,) I

I I I \\\\\\\\\\\\\\\\\\\\\\\\\\\\

(xb,y b)

I -

at h,, Y,) /

I1

I l I

I

____--

blob treated as sheet affecting only normal component of velocity of sheet at (x,,y,)

edge of visious

I I I l

I I

Fig. 4. Effect of a vortex blob on a vortex sheet.

sublayer

268

G.A. Sod / A compressible vortex method

Physically this is reasonable, since far downstream the flow can be viewed as fully developed independent of the flow upstream. So the “old” blobs can be advected with a velocity independent of time.

5. Compressible

vortex

and now

method

It remains to describe how the methods are coupled to produce a compressible vortex method. Suppose Uj[: and v:,~ are given. To start, at .vl= 0, U1y represents the initial conditions and in the absence of an initial vorticity distribution vz,, represents the initial velocity field. At time t = nk consider the piecewise constant function U”(x)=qY,

(i-:)h
(5 .la)

v,n(X)=v:,j’

(j-#
(5 lb)

To advance the solution to time t = (n + l)k, we must compute vj+r. For the jth vortex blob at center _Y,“,compute the local Reynolds number Re = Re,! = p”( x,“)uL/p in (4.6) where p”( 1,“) is obtained from (5.la) and the local free stream velocity vi, = v:( xJ’> in (4.4) using (5.lb). Once the vortex blobs have been advected and diffused, the third fractional step in (2.6) is performed. This consists of solving the ordinary differential equation a,,$= l(vp

x VP) -p

P2

Re-‘(vp

x Av).

(5.2)

This term represents the creation of vorticity along a front where there is a jump in the density. This is solved at the grid points (ih, j/z) E Qnh, of the random choice method using Euler’s method

[

5,, = k $(vp where vp

and vp

x vp)

- p

1

Re-‘(vp x Au)

ij

are approximated

using centered

‘xPij = (Pi+l,j

- Pi-l,j/2h

+ O(h’)>

ayPi, = (Pi,j+l

- P;,j_,/2h

+ 0( h2)

differences,

that is,

which introduces no numerical diffusion. In most cases, the second term on the right of (5.2) is O(Re-‘) and is negligible. However, where it is not, Au is approximated by the five-point difference, Avij=

(‘l+l,j

+ ‘i,j+l

- 4U,j + ~,-l,j

+ u,,j_l)/h

+ 0( h*).

At first glance it might seem that this would result in a large number of vortex blobs being introduced. However, vorticity is only produced where the gradients in (5.2) are large, that is, across a front. At each grid point where 1tij 1 > tmi,,, the vorticity is divided among an integer Mjj number of vortex blobs, each with circulation &,/Mjj where 1tij 1/Mij G (,,, and with sum of circulations at grid point (ih, jh) totaling Ei,_ These vortex blobs participate with the other vortex blobs.

269

G.A. Sod / A compressible vortex method

Compute

the velocity OiG’ = F

field induced

by the vortex

blobs on the grid fi2,,

tj~~(s(x~j-x~)~

I=1

where xlj = (ih, jh) E L?,,. This is added ,I;+”

to v:,~,

= Vfi, + v;;;r

which replaces the velocity field v:, in Ui7. Let U.(ni1) denote The random choice method is used to solve (2:4) with Ql” the boundary condition (2.3) 4,. to obtain

subject to

n I cm,,= 0,

with velocity

u.7” V n+l ClJ

this modified vector. as the initial condition

=

vyj+l

-

n+l ‘d,j

field v:,+‘. From which we obtain, .

6. Application to the interaction of an oblique shock wave with a boundary layer The essential features of the interaction of an oblique shock wave with a laminar boundary layer are well understood. The early experiments of Chapman et al. [2] and Hakkinen et al. [16] provided surface data. More recent experiments of Degrez et al. [13] provide both flow-field and surface data. This problem involves two of the most difficult features to model numerically, shock waves and separated flows. (A schlieren photograph is shown in Fig. 5 along with a line drawing depicting the essential flow features in Fig. 6.) For this reason it has been chosen as a test problem by Beam and Warming [l], Dawes [12]. Hanin et al. [18], and MacCormack [20,21].

REFLECTED SHOCK / i/

EXPANSION

POINT

Fig. 6. Line drawing

of an interaction

of an oblique

POINT

shock wave with a boundary

layer.

270

G.A. Sod / A compressible vortex method

M’rlO

,--____ i ‘\ I L\,, I

P = PO

I-+

p=po

I -

(a) upstream

Cd) y = t-l=rl,>

P-P,,

_____\---____

p=

Pl

-L--’

\L\

; (b) downstream I

outflow

I

condition

\ ‘1

\

\

\

I-

I L__,

\

\

\

\\

e<,\; \ \ \ \ Cc) Flat plate no-slip y = 0 adiabatic wall

Fig. 7. Computational

\\\

domain and boundary conditions.

We choose Hakkinen’s data because the experimental data includes both plate surface pressure and skin friction. The Mach number is 2.0, the oblique angle is 32.585 O, and the pressure ratio (final pressure after reflection to initial pressure) is 1.4. The free-stream Reynolds number is 2.96 x 105. The boundary conditions are: (a) the upstream values were initially set at the uniform flow and thereafter held fixed; (b) the downstream boundary was treated as an outflow boundary; (c) the boundary containing the flat plate satisfies the “no-slip” condition and it is treated as an adiabatic wall, and (d) the values of the boundary opposite the flat plate ( y = cc) are initially set using the inviscid oblique shock wave theory such that a shock of given strength and point of incidence on the plate is obtained, and thereafter were held fixed (see Fig. 7). A uniform grid spacing h = 0.05 was chosen, which was also the length of the vortex sheets. Let X,, denote the point of incidence on the plate, X, denote the point of separation of the boundary layer, and X, denote the point of reattachment of the boundary layer. The points of

‘1 0.2

0

I

I

I

0.4

0.6

0.8

I

I

1.o

1.2

X/X%

Fig. 8. Surface pressure.

I

1.4

I

1.6

271

G.A. Sod / A compressible vortex method 1.4

1

0

oooo~~~~

0

o”ooooooo 0.2

0.4

0.6

I

I

I

I

0.8

1.0

1.2

1.4

I 1.6

XIXCJ

Fig. 9. Skin friction. 2.0-

: 0

(a)

: : :

1.5-

0

: 0

F

: 0 0

i.o-

F

,” 0 0 : 0.5-

0

0 0.0-f

O

0.0

I

0.2

0

0

0.4

0

0

0

0”

I 0.8

I 0.6

r 1.0

2.0- (c)

0.0

0.2

0.4

0.6

0.8

0.6

0.8

0 0 0 0 : :

1.5-

: 0" 0 0

F

l.O-

0.5-

0

0

0

0

0

0

0

0

0

,” 0 0 0

0

0” 0 0.0-r

Oo

0.0

,

I 0.4 u/u0

I 0.8

0.0

0.2

0.4

u/u0

Fig. 10. Velocity profiles at (a) X/X,, = 0.5, upstream of separation; (b) X/X,, = 0.89, near point of separation; (c) X/X,, = 0.97, within region of separation; (d) X/X,, - 1.4, downstream of reattachment.

272

G.A. Sod / A compressible vortex method

separation and reattachment of the boundary layer correspond to the points when the skin friction cr vanishes. The computed values of these points (using linear interpolation) are = 1.181, which are in excellent agreement with Hakkinen’s results, X,/X,, = 0.858 and X,/X,, X,“/X,“, = 0.881 and X,“/X,“, = 1.14, respectively. The surface pressure and the skin friction are plots versus X/X,, depicted in Figs. 8 and 9, respectively. In Fig. 9, the solid line represents the skin friction as predicted by boundary layer theory (see [22]). The results are in excellent agreement with Hakkinen’s data, except that the plateau in the pressure profile is not as pronounced as in the experimental results. This possibility is due to the relative coarseness of the grid. The accurate prediction of separation and reattachment, despite the coarse grid, is due to the vortex elements (both sheets and blobs). In Fig. 10 we present the flow-field results, y versus u/Q, where y is the similarity variable ylTu,/xv and U, is the upstream velocity. These velocity profiles are presented at four locations: (a) X/X,, = 0.5, upstream of separation; (b) X/X,, = 0.89, near the point of separation; (c) X/X,, = 0.97, within the region of interaction; and (d) X/X,, = 1.4, downstream and reattachment. The backflow is clearly visible in the velocity profile at X/X,, = 0.97.

References [l] R.M. Beam and R.F. Warming, Alternating direction implicit methods for parabolic equations with a mixed derivative, SZAMJ. Sci. Statist. Comput. 1 (1980) 131-159. [2] D.R. Chapman, D.M. Kuehn and H.K. Lorson, Investigation of separated flows in supersonic and subsonic streams with emphasis on the effect of transition, NACA TR 1356 (1958). [3] A.Y. Cheer, Numerical study of incompressible slightly viscous flow past blunt bodies and airfoils, SIAM J. Sci. Statist. Comput. 4 (1983) 685-705. [4] A.J. Chorin, On the convergence of discrete approximations to the Navier Stokes equations, Math. Comp. 23 (1969) 341-353. [5] A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid Mech. 57 (1973) 785-796. [6] A.J. Chorin, Random choice solution of hyperbolic systems, J. Comput. Phys. 22 (1976) 517-533. [7] A.J. Chorin, Random choice methods with applications to reacting gas flow, J. Comput. Phys. 25 (1977) 253-272. [8] A.J. Chorin, Vortex sheet approximation of boundary layers, J. Comput. Phys. 27 (1978) 428-442. [9] A.J. Chorin, Vortex methods and boundary layer instability, SIAM J. Sci. Statist. Comput. (1980) l-21. [lo] P. Colella, Glimm’s method for gas dynamics, SIAM J. Sci. Statist. Comput. 3 (1982) 76-10. [ll] B. Davari, Numerical study of viscous flow past a circular cylinder: application of Chorin’s method, Lawrence Berkeley Laboratory Rept. LBL-2480, Berkeley, CA (1973). [12] W.N. Dawes, Efficient implicit algorithm for the equations of 2D viscous compressible flows: application to shock boundary layer interaction, Znternat. J. Heat Fluid Flow 4 (1983) 17-26. [13] G. Degrez, C.H. Bocca Doro and J.I. Wendt, The interaction of an oblique shock wave with a laminar boundary layer revisited: an experimental and numerical study, J. Fluid Mech. 177 (1987) 247-263. [14] A.F. Ghoniem, A.J. Chorin and A.K. Oppenheim, Numerical modeling of turbulent flow via combustion tunnel, Philos. Trans. Roy. Sot. London Ser. A 304 (1982) 303-325. [15] J. Glimm, Solution in large for nonlinear hyperbolic system of equations, Comm. Pure Appl. Math. 18 (1965) 697-715. [16] S.K. Godunov, Difference methods for the numerical calculation of discontinuous solutions of the equations of fluid dynamics, Mat. Sbornik 47 (1965) 271-306. [17] R.J. Hakkinen, I. Greber, L. Trilling and S.S. Abarbanel, The interaction of an oblique shock wave with a laminar boundary layer, NASA Memo 2-18-59W (1959). [18] M. Hanin, M. Wolfshtein and V.E. Landau, Numerical solution of Navier-Stokes equations for interaction of shock wave with laminar boundary layer, ICAS 74-17 (1974).

G.A. Sod / A compressible vortex method

273

[19] M. Lo&e, Probability Theory (Princeton University Press, Princeton, NJ, 2nd ed., 1960). [20] R.W. MacCormack, Numerical solution of the interaction of a shock wave with a laminar boundary layer, in: Proceedings 2nd International Conference on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics 8 (Springer, Berlin, 1971). [21] R.W. MacCormack, A numerical method for solving the equations of compressible viscous flow, AZAA J. 20 (1982) 1275-1281. [22] H. Schlichting, Boundary-Layer-Theory (McGraw-Hill, New York, 6th ed., 1968). [23] J.A. Sethian, Turbulent combustion in open and closed vessels, J. Comput. Phys. 54 (1984) 425-456. [24] G.A. Sod, Numerical study of a converging cylindrical shock, J. Fluid Mech. 83 (1977) 785-794. [25] G.A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 27 (1978) l-31. [26] G.A. Sod, Automotive engine modeling with a hybrid random choice method I, SAE Paper No. 800288 (1980). [27] G.A. Sod, Automotive engine modeling with a hybrid random choice method II, SAE Paper No. 820041 (1982). [28] G.A. Sod, Numerical Methods in Fluid Dynamics, Vol. 1 (Cambridge University Press, New York, 1990).