Stream function vorticity revisited

COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

STREAM

FUNCTION

C.A.J. ~epar#~ent

MECHANICS

FLETCHER

ofmechanical

Engineering, Received

Revised

manuscript

AND

ENGINEERING

VORTICITY

41 (1983) 297-322

REVISITED*

and K. SRINIVAS Universi~ 4 April received

of Sydney, NSW ZUU6, Australia

1983 29 July 1983

A three-level consistent time-split group finite element formulation is described for a stream function vorticity representation of incompressible laminar separated flow. A consistent time-splitting is achieved by explicitly extracting d~rectiu~a~ mass and di~erence o~rators. It is shown that application of linear rectangular elements to the convective terms with a non-uniform grid introduces an artificial viscosity, but that a simple modification removes this effect. The three-level fully implicit marching algorithm is shown to converge to the steady solution faster than a Crank-Nicolson scheme when combined with a banded Gauss elimination procedure along each grid-line in turn. The present method is applied to the flow past rearward- and forward-facing steps for which the vorticity solution is known to be singular at the sharp corners. To isolate the sharp corners from the computational domain a surfuce layer is introduced across which the computational solution is matched to an analytic solution for the corner flow. A one-dimensional integration of the energy function and transverse extrapolation is used to obtain the pressure distribution on the body. Results, for a range of step-height Reynolds numbers of up to 250, demonstrate good agreement with corresponding experimental results.

1. Int~du~tion Since the early 1970s there has been a noticeable shift of interest from vorticity/strearr function (l/JI) formulations to primitive variable (u, u, p) formulations for incompressible viscous flow. This shift has been apparent for both finite difference and finite elemen methods. The early status of both formulations, using finite difference methods, is well documentec by Roache [l]. Since that time finite difference primitive variable formulations have beer developed considerably, particularly for three-dimensional internal flows [2-51. For three. dimensional flows vorticitylvelocity potential formulations require six equations per grid point [6]. In contrast, four equations are required by primitive variable formulations. However vorticity/velocity potential formulations have been used in three dimensions f7]. Most finite difference formulations for incompressible viscous flow treat the problem a! being governed by parabolic (time-like) partial differential equations. This occurs naturally fat unsteady problems and is deliberately induced in the pseudo-transient [7] and parabolisec Navier-Stokes formulations [8] for solving steady flow problems. In contrast, finite element formulations for steady viscous flows [9, lo] have preferred to deal directly with the elliptic form of the governing equations. For steady flow finite element vorticity/stream function formulations an iteration is usually *The authors Scheme.

are grateful

0045-7825/83/$3.00

for the continuing

@ 1983, Elsevier

Science

support

Publishers

they have received

B.V. (North-Holland)

under

the Australian

Research

Grants

298

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

introduced to handle the vorticity boundary condition at a solid surface; although it appears that this iteration can be avoided [ll]. This iteration is an addition to the iteration required to solve the nonlinear global stiffness matrix equation. The vorticity/stream function formulation does have the major advantage of avoiding the continuity equation and the explicit appearance of the pressure. More recently, primitive variable finite element formulations have introduced penalty functions to avoid explicitly treating the pressure and the continuity equation [12, 131. For external flow problems, i.e., where part of the computational domain is bounded by the freestream, the vorticity goes to zero as the freestream is approached if there is no source of vorticity in the freestream. For the representative problem of the flow about an isolated body, the region away from the body, where the vorticity is zero, is governed by the Laplace equation for the velocity potential. Consequently, it becomes computationally attractive to solve the Navier-Stokes equations (in vorticity transport form) in the small region close to the body where vorticity is non-zero and to match the solution to the velocity potential solution in region can be much smaller when the vortithe outer region. The inner computational city/streamfunction formulation is used than when a primitive variable (u, II, p) formulation is used. This approach is followed in the integro-differential formulation of Wu [14]. The pseudo-transient approach has been used previously with a finite element formulation for compressible, viscous flow [15, 161 and the extension to incompressible flow using a pseudo-compressible prescription for the density or pressure [17, 181 is available. Here we re-examine the vorticity/streamfunction approach within a framework that will embrace both finite element and finite difference interpretations. A primary aim is to develop a computationally efficient algorithm without sacrificing the accuracy that can be obtained with the finite element method. An efficient pseudo-transient time-split method is applied to both the vorticity transport equation and the streamfunction equation. However, the stream function equation is only partially converged at each time step. The velocity components are retained explicitly in the vorticity transport equation and the velocity field is obtained from the streamfunction solution using an accurate implicit algorithm. Once the solutions for the vorticity and velocity fields have converged, the pressure solution is obtained by a one-dimensional integration of a total energy function parallel to the surface and an extrapolation to the body surface. The present formulation will be demonstrated on flows about bodies with both forwardfacing and rearward-facing corners. However, for this situation the vorticity becomes singular at the sharp corner. This difficulty is overcome by isolating the sharp corner from the computational domain. This is brought about by introducing a thin surface layer, i.e., thin compared with the local mesh, adjacent to all solid surfaces. The computational domain terminates at the outer edge of the surface layer. The solution at the boundary grid point adjacent to the sharp corner is matched to the analytic solution at the sharp corner [19]. Physically the presently considered flow problems may be interpreted as ones in which vorticity is generated at the solid surfaces and diffused and convected to the rest of the flowfield. However, close to the solid surface the convection is primarily parallel to the surface and the diffusion is normal to the surface. Consequently, the solution is readily extrapolated across the surface layer. The new vorticity/velocity/streamfunction formulation is described, in detail, in Section 2 and applied to flow problems with forward-facing and rearward-facing steps in Section 3.

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

299

2. New formulation

The governing equations for two-dimensional dimensional, unsteady conservation form)

incompressible

viscous flow are (in a non-

(1)

a251,a’4 _ ax*

1

5

(2)

’

ay*-

where

In (1) Re is the Reynolds number. Boundary conditions for the above equations will be considered in Section 2.5. The general approach to discretising (1) and (2) will be to apply a Galerkin finite element formulation with linear Lagrange rectangular elements. Neighbouring elements may vary in size which produces the nonuniform, rectangular mesh shown in Fig. 1. By introducing a group [20] finite element interpretation, the nonlinear convective terms in (1) can be handled accurately without the usual loss in economy associated with the conventional finite element treatment of convective noniinearities 1211. I-l J+ 1

1 J+

1

I1

I

I+1 J+l

‘y *Y

l-1

I

J

J

I+1 J 1

*Y

I

I+1 J-l

J-l

i

1

rxAx-

-AX Fig. 1. Global

mesh definition.

300

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

Application of the group finite element formulation with linear elements to (1) produces the following equation (for internal nodes):

Lagrange

M,OM~4+M,OL,uT+M,OL,ut-~(M,OL,+M,OL,,)f.=O

rectangular

(4)

where Mx = 6,

$1+ rx), kx> ,

Lx=~{-l,O,l},

Lx, =

&

M; = {:r,, f(1 + r,,), $} ,

(5)

(6)

L:+{l,o;-l},

(1, -(l + l/r,),

l/r,] 1

Lb, = &

{l/r,, -(l + l/r,),

l>

and where @ denotes the tensor product and r,, r, are defined in Fig. 1. The above operators have been obtained by applying the group finite element in the conventional manner and then dividing the resulting ordinary differential AxAy. If the muss operators, M, and MY, are lumped, i.e.,

(7)

formulation equation by

centered finite difference schemes result. For the difference operators, Lx, and LyY, this is more easily seen for a uniform grid, r, = r, = 1. The mass operators effectively transform three-point ‘finite difference’ formulae for the spatial terms into nine-point finite element formulae. In three dimensions (1) would be replaced by

at

dF aG at+ax-ay=Re

1

(8)

and u, v, where F = u& - wlx and G = wf; - v&, lx, &,, lz are the three vorticity components w are the three velocity components, Additional vorticity transport equations would be required for lx and 5;. If the group finite element method is applied to (8) the result would be (for internal nodes)

=~{M’oM,oL,,+M,OM,OL,,+M,OM,OL==}~.

(9)

In (9) the operators are as given by (5~(7) with an equivalent form for M,, etc. The mass operators M,, MY and M, have a similar role to that in two dimensions. For the spatial terms, three-point ‘finite difference’ formulae are transformed into twenty seven-point finite element formulae. The structure that permits mass and difference operators to be separated is inherent in Lagrange shape functions. However, such is not the case for serendipity shape functions [22].

C’. A.J.

Lagrange

shape functions

Fletcher,

K. Sriniuas.

can be written

Srream

function

as a product

corticity

30

rerisited

of one-dimensional

functions,

1

i.e..

N,(x. y. z) = rv~r’(x)lvg~‘(~)ru~=~(z).

(10)

For linear clcments. N!.“(X) = A(1 + x,x) . with - 1 s x d 1 and x, = kl. Equivalent formulae can bc written Components of the various operators in (S)-(7) are evaluated from

down

for N!y) and

IV:“.

(11)

where 1, denotes summation over all elements. Thus, if quudratic Lagrange elements are used on a uniform expressions are obtained: (1) Corner

grid. the following

operators

nodes:

M,-A{-1.3.8.2.-l},

L, = &1.-4.0.4.-l}. (2) Midside

L,, = &-{-1.X.-14.8.-1).

nodes:

M, = i${O,I, 8. 1. 0) , L, =

&

(0. - 1, 0. 1, 0) .

L,, = g&

(0. 1. -2, l,()I.

In practice there may not be too much interest in higher-order elements since one expects them to be significantly less economical due to the increased connecriciry associated with the mass operators [20]. This problem is aggravated in three dimensions. The question of whether the increased connectivity (and hence increased execution time) is justified by the extra accuracy will be considered in a associated with the mass operators related paper 12.31. 2. I. Vorticity equation The explicit

extraction

of the mass operators

also facilitates

the construction

of the efficient

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

302

time-splitting scheme that is at the heart of the new formulation. The time-split scheme will be developed for the vorticity equation. In the conventional finite element approach to unsteady problems, e.g. [24], the marching in time of the system of ordinary differential equations, (4) requires a single factorisation and storage of the mass matrix M, @M,,. At every time step matrix multiplication is required involving the factored form of the mass matrix, if the spatial terms are treated explicitly. However, such a scheme has a severe stability restriction on the maximum time step. If the spatial terms are treated implicitly, the stability restriction is avoided but an augmented mass matrix (i.e. including contributions from spatial terms) must be factorised at every step. If a global mass matrix must be factorised, then it is preferable to apply a variant of Newton’s method, e.g. as in [lo], directly to the steady equations. Even so, execution times are typically dominated by the repeated factorisation of the Jacobian (effectively the augmented mass matrix without the time-derivative terms). Here we construct an implicit time-marching scheme but introduce a time-split procedure that avoids the appearance of global matrices. Instead, subsystems of equations, associated with each gridline, are solved sequentially at each time step. Since the subsystems are banded close to the diagonal, Gauss elimination is particularly efTicient; in fact it is linear in the number of equations in the subsystems. For linear elements, the subsystems are tridiagonal; for quadratic elements the subsystems are pentadiagonal and so on. Strictly, for higher-order elements the structure of the subsystem equations depends on whether the equation originated at a corner node or midside node. It is possible [25] to take advantage of this regular structure in the implementation of the Gauss elimination process to reduce the execution time further. To facilitate the introduction of the time-splitting and to show how the time derivative is handled, (4) is written as follows, ~~~~~~~~+(l-*)~~=~RHS”+‘+(l-~)RHS~

(12)

where Al n+l_ - {“+‘-in, RHS=~(M,OL,,+M,OL,)f_MOL,ui.-M,OL,o5.

(13)

In (12) cy and /I are parameters introduced to weight the different time levels it and n + 1 at which Al and RHS are to be evaluated. Everything up to time level n is assumed known. After rearrangement, (12) will be used to obtain the correction Ap+l. Since Are1 will be obtained through solving a linear system of equations, it is necessary to linearise RHS”“. Therefore, a Taylor expansion is made about time level it, giving RHS””

Equation gives

=

RHS”

(14) is truncated

+

(14) after the terms shown with an error O(A9). Substitution

into (12)

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

a M,@M,-At; (

P a(;w)Aln+l

= At RHS”s@ - (1 - a)M, @ MyAl”.

In RHSWP (13) all terms are evaluated at time level n except u and o which are evaluated time t” + /3At. By projecting u and 2, and evaluating as part of RHS, a scalar subsystem equations for AC”+l can be obtained. Equation (15) is approximated as follows, a M,-At; I

P .,L,,-Lxu 1 (

has introduced

(15) at of

)]0[M’-At~(~L,,-L,v)1A5”+1

= At RHS”@ - (1 - a)Mx@ k&AC”.

The splitting

303

an additional

(16)

term, - L,V Al”+l.

However, this term is of O(A’t) and is acceptable even for transient solutions if At is small. That the splitting is possible, for the finite element method, follows directly from the product nature of the Lagrange shape functions (10). The present splitting preserves the mass operators, M, and My, which one would expect [26] to contribute significantly to the accuracy of transient solutions. A previous time-split finite element formulation [15] lumped the mass operators. However, this was primarily to achieve acceptable stability with a two-level scheme ((Y = 1, p = 0.5 in the equivalent of (12)). Equation (16) can be solved efficiently as a two-stage algorithm. In the first stage the following equation is solved,

>I

Al* = e

and in the second

RHS”.’ -

(17)

stage,

(18) At each stage the solution of (17) or (18) involves only operators associated with a particular coordinate direction. This permits a decoupling of the equations and the efficient implementation of a banded Gauss elimination along each grid-line in turn. In (17) and (18) many choices exist for (Y and /3. Two, in particular (a) (Y= 1, p = 0.5 and (b) a = 1.5, p = 1.0, are compared in Section 3.1. The linearisation introduced by (14) can be treated as the first-stage of a Newton-Raphson algorithm. Recently, a tensor product splitting, which is conceptually similar to (16) has been applied with the finite element method to the Newton-Raphson algorithm and used to compute shock-tube flow [35].

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

304

2.2. Stream function equation

In previous vorticitylstream function formulations (2) has either been solved iteratively using perhaps an AD1 or SOR scheme [I] or directly using a Poisson solver 1271. If (2) is solved, with no previous knowledge of the solution, the use of a Poisson solver is computationally more efficient than the use of an iterative scheme. However, in obtaining transient or pseudotransient solutions of (l)-(3) the solution $J” is clearly a close approximation to tin+‘, particularly if Ar is small. Consequently, relatively few iterations should be required to obtain n+l 31 * Here we use the same algorithm described in Section 2.1 to solve a ‘transient’ version of (2) i.e.,

a# a’+ f..a2+

dt- ax=

--

ay2

5

(1%

.

This approach has been used previously in the method of false transients [7]. Following the same steps as outlined in Section 2.1, the final algorithm, equivalent to (17) and (18) takes the form

and

(21) In (20) and (21) AT is a pseudo time step that permits the iterative evaluation of (20) and (21) at each time step At. In principle, at every time step At the stream function solution is started from +O= 9” and iterated until the steady-state residual (Poisson equation for $), RHS = [MY@ Lx, + M, @ L,,]$”

- M, @ MJ*+l

,

is sufficiently small. In practice, faster overall convergence to the steady-state vorticity solution is obtained by initially extrapolating, i.e. 4fi”= 24~”- Jt’-*, and by only iterating (20) and (21) four or five times for every step of (17) and (18). 2.3. Evaluation of velocity components Some early applications of the vorticity/stream function formulation [l] have avoided the explicit appearance of velocity by substituting (3) into (1). This is economical of memory and slightly more economical in execution than the present procedure. However, it is not as accurate as the method described here. Application of the Galerkin finite element method with one-dimensional linear elements to (3) gives MYu = L,$

and

M,v = --L&I

P-4

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

305

where M,, MY, Lx and L, are given by (5) and (6). Equations (22) connect adjacent nodes along a single grid-line implicitly. Consequently, the same banded Gauss solver mentioned in Section 2.1 is used for each grid-line in turn. On a uniform grid, equations (22) have a fourth-order truncation error. 2.4. Truncation error on a non-uniform grid By using linear finite elements in the formation of (4) (20) and (22) we are assured of convergence to the exact steady-state solution like O(A*x, A’y), as Ax, Ay + 0. This result does not indicate the absolute magnitude of the error on a specific mesh. For convergence to the exact steady-state solution, we also require that the algebraic equations, formed by application of the finite element method, are consistent with the governing partial differential equations. For the rectangular nonuniform grid considered here (Fig. l), the consistency can be established by a Taylor expansion about the Galerkin node. For the steady vorticity transport equation, this gives

+(l+rX)(l+rY)[-$@+%$-$-(~+$) +

iAx(r, - 1) $

(US) + iAy(r, - I) y

+ O(A*x, A2y

)] .

(23)

As expected, the finite element expression is consistent with (1). However, the leading terms in the truncation error are O(Ax, Ay) on a nonuniform grid. If the u and v velocity components are locally frozen, i.e., removed from inside the second differential, the truncation error can be interpreted as a negative viscosity. Thus, if the finite element equation is assumed to be second-order accurate, then it is equivalent to the governing equation,

Clearly, computational solutions at any Reynolds number will approximately correspond to physical solutions at a higher Reynolds number. The approximation will not be very good since Ax, r,, Ay, r,, may vary throughout the computational domain. For high Reynolds numbers, the spurious terms may well be of comparable magnitude to the physical viscosity. If the mass operators in (4) are lumped to generate an equivalent finite difference scheme, the leading terms in the truncation error are

i.e. first-order dissipative and dispersive terms are introduced. The additional terms in expression (23) have been interpreted as a negative viscosity. This is appropriate to a mesh which grows in the positive x or y direction. However, in modelling

C.A.J.

306

Fletcher, K. Srinivas,

Stream function

vorticity revisited

flows about an isolated body, or the forward-facing step considered in Section 3.3, it is desirable to have a fine mesh adjacent to the body and a growing mesh in the negative x-direction ahead of the body. In this situation the first-order truncation error provides positive dissipation. The rapid change of a convected variable in front of a body with a coarse grid is precisely the situation where one would expect ‘large cell-Reynolds number’ wiggles. Clearly, the use of a non-uniform mesh in the negative x-direction tends to suppress such wiggles without, necessarily, providing the correct solution. The first-order truncation errors discussed above have the same origin as the first-order truncation errors associated with linear upwind differencing. Namely, that the difference operator Lx in (6) becomes asymmetric about the (i, j)th node. This suggests that the remedy is similar. That is, to introduce a Petrov-Galerkin [22, 281 formulation to compensate for the first-order terms in the truncation error. Here we achieve the same ends by adopting the tau strategy [22, 291. That is, we perturb the original problem in such a way that our approximate solution to the perturbed problem is a closer approximation to the exact solution of the unperturbed problem. In practice we apply the group finite element formulation to the following equation in place of (1)

$f+2

(~5) - iAx(r, - 1)

2

(uL)+ y-iAy(r,-I)$(v+-$$+$=O.

(25) The result is, in place of (4) M, @ Myi + My @ L,ul + M, @ -M,@L,,

L,vc

(,+6A~(~X-1)~5)-M~@LYy(~+iAy(~Y-1)v~}=0. ’ i

This equation is consistent with the original partial differential A Taylor expansion of the algebraic Poisson equation,

produces

a truncation

error

-i$(l + r:)l(l

+ rx) A2x $$

equation

(26) to order

(A%, A’y).

equal to - &(l + r;)/(l + r,) A*y 3,

i.e., the algebraic Poisson equation is consistent to order A Taylor expansion of (22a), i.e., M,u = L,t+b,gives [M,,u - L,I)] =,$(l + rY)[ u - $$ - &Ay(r, - 1) 3

A2x, A2y.

+ &A’y (r, - l)(r; + 1) $$

+ W”Y

,] .

(27)

307

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

It can be seen that (22a) is first-order accurate on a nonuniform grid but fourth-order accurate on a uniform grid. Following the tau procedure outlined above, (3) is replaced by and

~1= -$+iAx(r,

- l)$,

(28)

M,u = -(Lx - kAx(r, - l)L,,)cC,.

(29)

and (22) by M,u = (L, - iAy
and

The results produce equations that are consistent to second-order on a nonuniform mesh and consistent to fourth-order on a uniform mesh. The above treatment of (l)-(3) can be generalised by replacing first derivatives in the governing equations by $-iAx(r,

- l)$

and applying a conventional elements.

and

d-‘A dy 6 y(+)$

group finite element

formulation

with linear rectangular

finite

2.5. Boundary conditions The precise boundary conditions will depend on the problem considered. In Sections 3.2 and 3.3, respectively, the flow over a rearward-facing and forward-facing step are considered. The computational domains are shown in Figs. 2(a) and 2(b). Here we will consider the following classes of boundary condition: (i) solid surface, e.g. FE in Fig. 2(a), (ii) farfield boundary, e.g. AB in Fig. 2(a), (iii) upstream boundary, e.g. AF in Fig. 2(a), (iv) downstream boundary, e.g. BC in Fig. 2(a). The treatment of the corners at E and D will be considered in Section 2.6. 4

F’

““__r___“‘____~__“--_’

I

8

I I I

I I I I

I# i

I I I

I

I

I I

I I

6

A

c!

_;______-;-______+__--__

F

I IA’ I I I I I

I IF I I I I I

I I

I I

B

c’

Fig. 2. (a) Computational domain for flow over a rearward-facing step. (b) Computational domain for flow over a forward-facing step.

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

308

2.5.1. Solid surface A solid surface, u, v and J,!Iare known and l must be computed. If no surface layer (see Section 2.6) is used, u, v and Cc,are zero. If a surface layer is used, u, v and I,!Iare computed from the surface value of 5 (Section 2.6). Consequently, Dirichlet boundary conditions are provided for the calculation of $ (Section 2.2) and u, v (Section 2.3) and no special equation is required at a solid surface. To calculate 5 at a solid surface, (2) is used in the following pseudo-transient form,

(30) Application of the Galerkin finite element method and the time-splitting procedure produces an algorithm similar to (20) and (21). In the present situation this gives, for the first half time step, [a + PY AtIM, U* and for the second

= Y At[M, @Lx, + M, OL,,l+”

+ Y At[MJXv)+

M&u)1

3

(31)

half time step,

My Al”+’ = Al* .

(32)

The additional terms f(v) and g(u) arise from the application of Greens theorem to J2$/dx2 and a2$/‘ldy2 when the Galerkin node is on the boundary. For surface FE in Fig. 2(a), f(v) = 0 and g(u) = u,,/Ay. Equations (31) and (32) can be coupled directly to the system of interior equations (17) and (18). However, y must be less than unity to prevent instabilities from developing; typically, y = 0.1. At the surface, the form of the mass and difference operators is different from the interior form, (5) and (7). Thus, on DE in Fig. 2(b), M, = {a, 4, 0)

and

L,., = &

{I, -I,

O>.

(33)

2.5.2. Farfield boundary The boundary AB is placed far enough from the body that the vorticity is zero. In addition, u is set to the freestream value, i.e., u = 1, on AB. Consequently, Dirichlet boundary conditions are provided for the calculation of vorticity and velocity component u. The stream function 4 is calculated on AB using (20) and (21) except that only two elements contribute. As a result, the operators My and L,, have the following form, Mi={O,f,a}

and

L\,=-$--{0,-l,

l}.

Operators M, and Lx, are as in (5)-(7). Application addition of the term ArMXu,,/Ay to the right-hand on AB is calculated, conventionally using (29).

(34) of Greens theorem to J2+l’la2y requires the side of (20). After $ is calculated on AB, u

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

309

25.3. Upstream ~o~~da~ On AF in Figs. Z(a) and 2(b), u is specified as a boundary layer adjacent to F and as the freestream value, u = 1, from the boundary layer outer edge to A. The stream function rl/ on AF follows directly from (22). Thus, on AF a Dirichlet boundary condition is provided for the calculation of $ in the interior. On AF, (31) and (32) are used for the vorticity with f(v) = v,&Ax and g(u) = 0, and the following form for A4, and LX, (rX= l), K

= K4f, $1,

Lx

=&IO,-1311.

(35)

The operators M, and L,, are given by (5) and (7). A Dirichlet boundary condition for v on AF is required for the calculation of u in the interior using (29). In the present formulation vAF is obtained by a three-point extrapolation from the interior in terms of I& Thus,

(36)

V

with Ax = xi-i - x;_~. Alternatively,

a single element version of (29) could be used.

The stream function rt/ is calculated on BC using (20) and (21) with the addition of -h~M,,v~~/Ax to the right-hand side of (20). This arises from the application of Green’s theorem to d%,!~//lax*. Also, in (20) M, and LX, take the form (rX= l),

Mx = It, 4,01,

L,,-&-{I,-l,O}.

From the + solution, VBCis obtained using a three-point

(37) extrapolation

from the interior,

with Ax = Xi-1- Xi-z. The vorticity 6 is calculated on BC using (17) and (18) but with the L, term deleted, i.e., effectively ignoring the term #LJ~x~ in (1). It may be noted that d2[/ax2 is expected to be small far from the body. Close to the body a2L/ax2 will be small compared with J21/dy2. Consequently, any additional errors introduced through the treatment of the #l/ax* term will have a small influence on the l solution, even locally. 2.6. Surface

Layer arzd treatment of corners

For the problems shown in Figs. 2(a) and 2(b), the solid surfaces are sources of vorticity which is diffused and convected away by the fluid flow. Consequently, the vorticity has a maximum value at the solid surface and may have large gradients close to the solid surface. To

310

C.A.J. Fletcher, K. Sriniuus, Stream function uorticity revisited

obtain accurate solutions implies having grids elsewhere. However, very close to a solid surface magnitude analysis indicates that d{/an (Fig. 2(a)). For example, for surface FE

a fine grid adjacent to the solid surfaces and coarser the velocities will be relatively small and an order of is a constant across a small layer E next to the wall in Fig. 2(a),

(39) From the de~nition of vorticity and stream function it follows that ui = E& - 0.5ke2,

(40)

$!lj= 0+5&Uj

(41)

and

where node j is on the boundary of the computational domain, i.e., displaced from the solid surface by the width of the surface layer E. A typical value is E = 0.2 Ay. In practice, L$ is evaluated as in Section 2.5.1 and (40) and (41) give Uj and $j, respectively (Us= 0 for a horizontal surface), In (39) k is evaluated from a one-sided three-point formula for al/-/y. An advantage of introducing the surface layer is that the sharp corners, D and E in Figs. 2(a) and 2(b), are placed outside of the computational domain. This is important because it is wellknown [19] that the vorticity is singular at a sharp corner; whereas the velocity and streamfunction are regular. The analytic solutions adjacent to the sharp corner E in Fig. 2(a) can be obtained [19] as

cos mp ui,j -_-yi,j 51 --22’““-“B i 1

W%-1

cos(m

and

-

2)p c

’

(43)

where m = 1.5444837, /3 = $7~and B is chosen to match the solution away from the singularity, i.e., in the computational domain. In practice, li.1 is computed using the appropriate form of (3f) and (32) and Ui,j,Ui,jand +i,j are obtained from (42)-(44). A similar procedure is used for node D and for the forward-facing step shown in Fig. 2(b).

Once a converged solution for the vorticity and velocity fields have been obtained, the pressure field can be evaluated. Traditionally [l] this requires the solution of a Poisson equation for the pressure with Neumann boundary conditions evaluated from the velocity solution. However, often the pressure is only required on the surface of a body which forms the boundary of the computational domain, e.g., FEDC in Fig. 2(a). At first sight one could

311

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

integrate an ordinary differential equation for the pressure, starting at F and eventually arriving at C. However, one expects errors to accumulate with particularly large contributions arising at nodes E and D. Here we introduce the variable, H = ;+

;(u* + v*)

where H can be thought written,

aH -= ax

1 al ---

of as a total energy

function.

The momentum

equations

can then be

v(,

Re ay

dH

1

x

ay-- -Reax+

(47)

lg.

For flows parallel to solid surfaces at high Reynolds number, the viscosity is only significant in the thin boundary layer adjacent to the surface. Outside of the boundary layer, H is constant if the flow is homenergic [30] and P follows from (45). Therefore, the following strategy has been used to obtain the pressure at the solid surface. An integration of (46) is made along a grid line parallel to surfaces FE and DC, and (47) parallel to ED, from boundary AF to boundary BC. A second integration is made along an adjacent grid-line and an extrapolation made back to the surface. If the pressure is required in the domain, the following equation is solved,

as

I

a*H

a(4 _ s

w4

ax*ayZ=ay

ax

(48)

’

using the vorticity and velocity solutions to provide the Neumann boundary conditions for H. Due to the behaviour of H in the inviscid flow region, the computational domain can be reduced to include only the boundary layer, separated flow and wake region, at high Reynolds number. The same approach as indicated in Section 2.2 is used to solve (48). Thus a ‘transient’ adaptation of (48) is introduced,

dH _- a*H+ -a*H ata? a$

s

(49)

.

Application of the Galerkin finite element method with linear rectangular elements and the time-splitting procedure described in Section 2.1, produces the following two-stage algorithm. For the first half time step, P Mx-ArcyLxx

1AH*=~[M,OL,,+M,OL,,IH” “s-mAHm]. a cl

312

C.A.J.

Fletcher, K. Srinivas, Stream function vorticity revisited

For the second half time step,

i

M,-17$L,,lAH”‘“=AH*.

The transient integration is continued pressure P then follows from (45).

(51) until the solution

H””

is not longer changing. The

3. Results and di~ussion Here, the new formulation, described in Section 2, is applied to the flow past rearwardfacing and forward-facing steps (Figs. 2(a) and 2(b)). These two problems have regions of separated flow and locations where the vorticity becomes singular. Consequently, the problems are considered a severe test. Solutions will be presented for a range of Reynolds numbers for which laminar flow is expected. Although the geometries shown in Figs. 2(a) and 2(b) allow rectangular elements to be constructed, irregular domains can be handled by introducing generalised coordinates and transforming to a uniform domain before the time-split group finite element formulation is applied. The use of the group formulation is crucial if an economical algorithm is to be obtained in generalised coordinates. Early results, applying the present splitting in generalised coordinates for the flow adjacent to an aerofoil trailing edge, are reported by Srinivas and Fletcher [36]. 3.1. Time marching

algorithm

In Section 2.1 the following time-marching

MxOMy

I

\

i\Y*+l

a _+(l-“)T

U”

algorithm was introduced,

= p RHS”” + (1 - ,6)RHS”.

(12)

The same marching algorithm is used to solve the Poisson equation for the streamfunction Z,!J in Section 2.2. The Crank-Nicolson (trapezoidal) scheme (a = 1, p = 0.5) has been widely used in the past because it possesses a truncation error of O(A’t) and is unconditionally stable, in the von Neumann sense, for many implicit formulations. However, it also demonstrates rather slow convergence when used with a pseudo-transient formulation, particularly as the steady state is approached for stiff problems. A number of different three-level schemes are discussed by Beam and Warming f31]. A particularly effective choice, for the pseudo-transient formulation, is (Y= 1.5, p = 1.0. This scheme is unconditionally stable and has truncation error of O(A”t). The convergence of this scheme with the number of time steps is shown in Fig. 3. ~~RHS~~_~ is the rms value of RHS in (13), i.e., the steady state residual. The results shown in Fig. 3 were obtained for the flow over a rearward-facing step with a step-height Reynolds number,

313

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

Fig. 3. Convergence TRAP 5 trapezoidal

comparison of (Crank-Nicolson)).

alternative

time-marching

algorithms

(3LFI = three-level

fully

implicit;

Reh===50, V

and a 59 X 59 grid. It is clear that the steady-state solution is reached with fewer steps when the three-level, fully implicit ((Y= 1.5, p = 1.0) scheme is used, although the trapezoidal (Crank-Nicolson) scheme converges faster, initially. The three-level scheme is slightly slower (15%) per time step due to manipulating the previous correction A[” in (12). Also, the three-level scheme does require the storage of two levels of vorticity and stream function. However, the extra level of stream function permits a projection forward in time which reduces the number of iterations required to solve (20) and (21). The three-level scheme has been used to generate all the results presented in Section 3.2 and Section 3.3. 3.2. Flow over a rearward-facing

step

The computational domain for this problem is shown in Fig. 2(a). Solutions for various grid sizes and Reynolds numbers have been obtained; grid details, of examples that are presented in Figs. 4-8, are given in Table 1. All lengths shown in Table 1 are nondimensionalised with the step-height (Fig. 2(a)).

314

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

Table 1 Grid details for results presented

Case Grid (AB x BC) Length AB Length BC

AXmin = (Aymin) A.&lI~X Aym.a

Boundary layer thickness (FF’) Re (in (1)) Step-height Re

in Figs. 4-f I

(RAF) 111x67 27.45 8.67 0.063 0.372 0.482 1.266 250-1000 63-253

C

D

(R”Ft

(RF)

(RF)

111 x 67 44.94 6.69 0.063 0.810 0.273 1.266 300 119

59x 59 88.75 10.73 0.114 2.830 0.467 2.273 600-1000 145-242

59x59 28.89 8.91 0.114 1.390 0.340 1.136 1%-600 53-211

(& 59x59 34.43 6.43 0.051 3.850 0.630 1.010 300 89

(F& 59x59 24.54 12.30 0.05 1 1.260 0.478 0.505 100-400 40-160

“RF = rearward-facing step. ‘FF = forward-facing step.

The minimum grid size, Axmin= Ay,in, has been used adjacent to E (Fig. 2(a)). The mesh grows in both x-directions with a geometric ratio r, to FA and D’B’. Between D’B’ and CB the x-mesh is uniform (rX= 1). The y-mesh is uniform from DC to F’C’ (the edge of the boundary layer FF’). From F’C’ to AB the y-mesh grows with a geometric ratio r,. For cases A and B in Table 1, the truncation error correction (Section 2.4) has not been used. For these cases, r, < 1.04, r,, d 1.07. Truncation error correction has been used for cases C to E for which r, d 1.15, r,, 6 1.15. The results presented in this paper have been obtained on a P~rkin-Elmer 3220 computer. For case A, with a 111 x 67 grid, the execution time is 45s per time-step. The execution time per time step is approximately linear in the number of active grid points (i.e., excluding the step). Typically, 500 to 1000 time steps are required to reduce the steady-state residual, RHS efficiency of the present scheme and the in (13) to less than 1 x 10m5. The computational influence of the various contributions, like the directional mass operators, to the effectiveness of the algorithm will be discussed in a related paper [23]. For various stations upstream and downstream of the step, the u velocity profiles are shown in Fig. 4. On the left-hand boundary a nominal Blasius boundary layer profile, of height = 1.266, is introduced. Because of the low Reynolds number a preliminary computation over a flat surface was made with a precise Blasius boundary layer profile at the left-hand boundary. The computed u velocity profile ten grid intervals downstream was used as the boundary condition on FF’ (Fig. 2(a)) for the present results (Fig. 4). The results shown in Fig. 4 correspond to a step-height Reynolds number, Reh = 126. The flow field is characterised by a bubble of slowly moving reversed flow just behind the step. Downstream of the bubble, the flow reattaches and eventually, 14.8 step-heights beyond the step, a boundary layer starts to redevelop. The computational domain for the present problem extends sufficiently for downstream that the beginning of the downstream boundary layer growth is captured. A number of experiments, e.g. [32, 331, have measured the length of the separation bubble (reattachment length) for laminar flow over a rearward-facing step. A comparison with bubble

315

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

X/h=

X/ h=Od 6

X/h=-044

1”

i

r”

I.90

X/h=700

X/h=2500

i

U/UINF

Fig. 4. u-velocity

profiles

downstream

of the rearward-facing

step, Re,, = 126.

lengths computed using the present method is provided by Fig. 5; grid details for the various cases are given in Table 1. Given the scatter in the experimental results, the agreement is considered satisfactory. The present results show a tendency for the bubble length to be greater, at the same step-height Reynolds number, if the ratio of the step height to boundary

11 1X67

STEP-HEIGHT

Fig. 5. Variation

of reattachment

REYNOLDS

length

MESH

NUtlBER.Reh

with step-height

Reynolds

number.

316

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

layer thickness is less. This same trend is apparent in the experimental results of Goldstein et al. [33]. The flow solution in the neighborhood of the sharp edge (E in Fig. 2(a)) is shown in Fig. 6. These results correspond to Re, = 63. They suggest that the flow is still attached at the corner and, in fact, separates from just below the corner. Solutions at higher Reynolds number suggest that the separation point moves closer to the corner as the Reynolds number is increased. This supports the results of Roache and Mueller [34]. However it should be remembered that the computational domain (Fig. 2(a)) in the present study is separated from the physical surface by a layer of thickness, E = 0.2 by,,“. In addition, to resolve the detailed behaviour close the corner at high Reynolds number, requires a finer grid than was possible in the present investigation. Vorticity is generated at the solid surface and diffused and convected to the rest of the flow field. The surface distribution (strictly at the edge of the computational domain) of vorticity is shown in Fig. 7 (S/h is a surface coordinate), for two step-height Reynolds numbers. The absolute magnitude of the vorticity is primarily a function of the boundary layer thickness (FF in Fig, 2(a)); the thinner the boundary layer the higher the surface vorticity. But for the higher Reynolds number case the edge of the computational domain is closer to the boundary layer edge. The qualitative behaviour is the same for both Reynolds numbers. The vorticity rises to a maximum positive value just ahead of the sharp corner E. It drops very rapidly as the corner is traversed and is negative adjacent to ED. This corresponds to an upward velocity next to ED. The vorticity rises to zero adjacent to D and becomes negative again adjacent to DC corresponding to reverse flow in the bubble. It rises to zero at reattachment and continues to

-----+-----+----+-

~-----+-----+----b---4

----+-

Fig. 6. Detailed

--+-+---+---+--+“-c----+

flow behaviour

near the sharp trailing

edge.

317

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited

S/h

Fig. 7. “Surface’ vorticity

distribution

for a rearward-facing

step.

rise downstream. Eventually, as C is approached, 6 begins to fall with increasing X; this corresponds to the growth of a new boundary layer. The pressure distribution at the solid surface (strictly at the edge of the computational domain) is shown in Fig. 8. This has been obtained using the procedure described in Section 2.7 ((46) and (47)). Like the vorticity distribution, the pressure distribution is qualitatively the same for both Reynolds numbers. A value of the pressure coefficient, Cp = 0, is expected for

-Reh = 126

gFE ‘-2.00

D T

a.00

2.00

4.00

6.00

0.00

to.00

12.00

1

14.00

l2.00

for a rearward-facing

step.

S/h

Fig. 8. ‘Surface’

pressure

distribution

I

I

12.00

I

1 20.00

31x

C.A.J.

Fletcher.

K. Srinivas,

Stream

function

vorticity

revisited

the flow past a flat plate (i.e., at F, Fig. 2(a)). The expansion of the flow around the corner (E in Fig. 2(a)) causes a slightly sub-freestream pressure which persists, within the separation bubble. until reattachment is approached (on CD in Fig. 2(a)). The pressure rises through reattachment to a small positive value and then subsides to zero as the boundary layer flow over a flat plate re-establishes itself (approaching C in Fig. 2(a)). The pressure rise during reattachment is a characteristic feature that has also been observed in experiments [32]. 3..?. Now ocer u forward-facing

step

The computational domain for this problem is shown in Fig. 2(b). As for the flow over a rearward-facing step, a uniform grid is used in the y-direction up to F’C’. Typically the boundary layer height only extends to half the height of the step DE. The grid grows geometrically in both x-directions from the line DA’. However, downstream of the line E’B’ the grid is uniform in the x-direction. The results presented in Figs. 9-l 1 arc based on case F (Table 1). In contrast to the flow over the rearward-facing step, the present flowfield is dominated by the sharp edge E (Fig. 2(b)). Even with a surface layer the vorticity adjacent to E reaches a very large value. The difference between the influence of a rearward and forward-facing step can be appreciated from the definition of vorticity.

The contribution au/+ is considerably larger close to and upstream of a forward-facing step than a rearward-facing step. The term c?ulJx is large and negative just upstream of a forward facing step. For a rearward-facing step au/ax is approximately zero. The value of [E (node adjacent to corner E, Fig. 2(b)). besides being large and positive for a forward-facing step, is also slow to converge. Consequently. the number of iterations to reach the steady state is controlled by the convergence of the solution adjacent to E. The variation of the u velocity component in the y-direction for various x-stations, is shown in Fig. 9. A small separation bubble occurs ahead of corner D (Fig. 2(b)), causing the u-component to be locally negative (x/h = -0.35). The flow does not separate behind the shoulder E at this low Reynolds number; instead there is a region of slow-moving fluid. At higher Reynolds numbers one would expect a well-defined bubble to develop immediately downstream of the shoulder (E in Fig. 2(b)). Behind the shoulder (E) a boundary layer profile adjacent to EC develops slowly with an increase in the fluid velocity above the undisturbed freestream value at about one step-height above EC. The variation of the ‘surface’ vorticity (strictly at the edge of the computational domain) is shown in Fig. 10 for two step-height Reynolds numbers (S is a surface co-ordinate). The qualitative behaviour is similar in both cases. Starting at F, the vorticity decreases rapidly in the direction of D (Fig. 2(b)), reaching a small negative value in the separation bubble ahead of D. As the face DE is traversed from D to E the vorticity initially takes a small negative value associated with the separation bubble. Thereafter, the vorticity grows rapidly reaching a maximum positive value adjacent to E. Immediately downstream of E the vorticity falls rapidly to a minimum and subsequently climbs slowly in the direction of C. The reason for this rise, which is more pronounced at Re, = 160. is not obvious.

C.A.J. Fletcher, K. Sriniuas, Stream function vorticity revisited

319

IL

8.

r; 8ri

X/h=-

X/h=04

X/h=-035

I.74

7

X/h=3.28

l

i

*

a

s

*

l

*

*

I

1

i

s

*

I *

l

i

a

X/h=12-35 l n

I

l

i

a ‘

.

8d

Fig. 9. u-velocity

profile

for a forward-facing

step, Reh = 160.

The pressure distribution over a forward-facing step is shown in Fig. 11 for two Reynolds numbers. In contrast to the situation for a rearward-facing step, the pressure distribution for a forward-facing step shows substantial variation. Starting at F (Fig. 2(b)) there is a steady increase in pressure as D is approached. The pressure continues to rise marginally from D to E. The acceleration of the flow around the shoulder at E causes the pressure coefficient to drop to a large negative value just downstream of E. From E to C a steady recovery in pressure takes place.

R53~84.82 ?I $

= 2 4.0

1

S/h

Fig. 10. ‘Surface’

vorticity

distribution

for a forward-facing

step,

320

C.A.J.

Fletcher,

K. Srinivas,

Stream function

vorticity

revisited

S/h

Fig. 11. ‘Surface”

pressure

distribution

for a forward-facing

step.

It is clear, particularly from Figs. 10 and 11, that to obtain an accurate solution adjacent to E would require a very refined grid locally, even when the solution is matched to the local analytic solution as in the present procedure. It is also clear that this problem gets worse as Re, increases. In contrast, the solution at a downstream-facing step (E in Fig, Z(a)) does not require a significantly refined grid as long as the solution is matched to the analytic solution at the corner across the surface layer.

4. Conclusion It is concluded that the present time-split group finite element formulation gives accurate and economical solutions for steady separated laminar flow for geometric configurations known to cause singular behaviour. The use of a surface layer to isolate the rearward-facing edge from the computational domain is effective without the need for additional local grid refinement. However, the singularity in vorticity at a forward-facing edge is considerably more severe and the surface-layer concept probably needs to be supplemented with a more refined grid locally. The fully implicit, second-order three-level marching algorithm is found to give faster convergence than a two-level crank-Nicolson scheme. The one-dimensional integration of the energy function and transverse extrapolation provides an economical and accurate method for computing the surface pressure distribution.

References (11 P.J. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, NM, 1972). [2] S.V. Patankar and D.B. Spalding, A calculation procedure for heat and mass transfer and momentum transfer in three dimensional parabolic flows, Internat. J. Heat Mass Transfer 15 (1972) 1787-1806. [3] W.R. Briley, Numerical method for predicting three dimensional flows in ducts, J. Comput. Phys. 14 (1974) S-28.

C.A.J. Fletcher, K. Srinivas, Stream function vorticity revisited [4] W.R. Briley and H. McDonald,

321

Analysis and Computation of Viscous Subsonic Primary and Secondary Flows, AIAA Paper 79-1453. [S] K.N. Ghia and J.S. Sokhey, Laminar incompressible viscous flow in curved ducts of rectangular cross-sections, J. Fluids Engrg. 99 (1977) 640-645. [6] H. Fasel, Recent developments in the numerical solution of the Navier-Stokes equations and hydrodynamic stability problems, VKI Lecture Series 78-4 (1978) l-90. [7] G.D. Mallinson and G. de Vahl Davis, The method of the false transient for the solution of coupled elliptic equations, J. Comput. Phys. 12 (1973) 43.5-461. [8] T. Davis and S.G. Rubin, Non-Navier-Stokes viscous flow computations, Comput. & Fluids 8 (1980) 101-131. [9] C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Comput. & Fluids 1 (1973) 73-100. [lo] M.S. Engelman, FIDAP: A fluid dynamics analysis program, Adv. Engrg. Software 4 (1982) 163-166. [ll] A. Campion-Renson and M.J. Crochet, On the stream function-vorticity finite element solutions of the Navier-Stokes equations, Internat. J. Numer. Meths. Engrg. 12 (1978) 1809-1818. [12] T.J.R. Hughes, W.K. Liu and A. Brooks, Finite element analysis of incompressible viscous flows by the penalty function formuIation, J. Comput. Phys. 30 (1979) l-60. [13] R.L. Sani, P.M. Gresho, R.L. Lee, D.F. Griffiths and M. Engleman, The cause and cure of the spurious pressures generated by certain FEM solutions of the incompressible Navier-Stokes equations: part 2, Internat. J. Numer. Meths. Fluids 1 (1981) 171-204. [14] A. Sugavanam and J.C. Wu, Numerical study of separated turbulent flow over aerofoils, AIAA J. 20 (1982) 464-470. [IS] C.A.J. Fletcher, On an alternating direction implicit finite element method for flow problems, Comput. Meths. Appl. Mech. Engrg. 30 (1982) 307-322. 1161 K. Srinivas and C.A.J. Fletcher, Finite element solutions for laminar and turbulent compressible flow, Internat. J. Numer. Meths. Fluids (1984), to appear. [17] A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1967) 12-16. [18] T.J.R. Hughes, T.E. Tezduyar and A.N. Brooks, Streamline upwind formuIations for advection-diffusion, Navier-Stokes and first-order hyperbolic equations, in: T. Kawai, ed., Finite Element Flow Analysis (Univ. of Tokyo Press 1982) 97-104. [19] H.J. Lugt and E.W. Schwiderski, Flows around dihedral angles, I. Eigenmotion analysis, Proc. Roy. Sot. Ser. A 285 (1965) 382-399. [20] C.A.J. Fletcher, The group finite element formulation, Comput. Meths. Appl. Mech. Engrg. 37 (1983) 225-244. [21] C.A.J. Fletcher, A comparison of finite element and finite difference solutions of the one and two-dimensional Burgers’ equations, J. Comput. Phys. 51 (1983) 159-188. [22] C.A.J. Fletcher, Computational Galerkin Methods (Springer-Verlag, Berlin, 1983). [23] C.A.J. Fletcher and K. Srinivas, Composite finite element/finite difference formulations for incompressible viscous flow (in preparation). [24] D.K. Gartling, Convective heat transfer analysis by the finite element method, Comput. Meths. Appl. Mech. Engrg. 12 (1977) 365-382. [25] C.A.J. Fletcher, A comparison of the finite element and finite difference methods for computationa fluid dynamics, in: T. Kawai, ed., Finite Element Flow Anaiysis (Univ. of Tokyo Press, Tokyo, 1982) 1003-1010. [26] A.J. Baker and M.O. Soliman, Utility of a finite element solution algorithm for initial-value problems, J. Comput. Phys. 32 (1979) 289-324. [27] P.N. Swartztrauber, A direct method for the discrete solution of separable elliptic equations, SIAM J. Numer. Anal. 11 (1974) 1136-1150. ]28] A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov-Ga~erkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Meths. Appl. Mech. Engrg. 30 (1982) 199-259. [29] C. Lanczos, Applied Analysis (Prentice-Hall, Englewood Cliffs, NJ, 1956). [30] L. Howarth, Modern Developments in Fluid Dynamics, High Speed Flow, Vol. 1 (Oxford Univ. Press, Oxford, 1953). 1311 R.M. Beam and R.F. Warming, An implicit factored scheme for the compressible Navier-Stokes equations, II: The numerical ODE connection, AIAA Paper 79-1446.

322

C.A.J.

Fletcher, K. Srinivas,

Stream function

vorticity revisited

[32] S.N. Sinha, A.K. Gupta and M.M. Oberai, Laminar Separating flow over backsteps and cavities, part I: Backsteps, AIAA J. 19 (1981) 1527-1530. [33] R.J. Goldstein, V.L. Ericksen, R.M. Olson and E.R.G. Eckert, Laminar separation, reattachment and transition of the flow over a downstream-facing step, ASME Trans Ser. D, J. of Basic Engrg. 92 (1970) 732-741. [34] P.J. Roache and T.J. Mueller, Numerical solutions of laminar separated flows, AIAA J. 8 (1970) 530-538. [35] A.J. Baker and M.O. Soliman, A finite element algorithm for computational fluid dynamics, AIAA J. 21 (1983) 816-827. [36] K. Srinivas and C.A.J. Fletcher, Time-split finite element method for compressible aerofoil trailing-edge flows, Proc. 3rd Internat. Conf. Numer. Meth. for Laminar and Turbulent Flow, Seattle, WA (3 Pineridge, Swansea, 1983).