Theory and implementation of a semi-implicit finite element method for viscous incompressible flow

Computers FluidsVol. 22, No. 6, pp. 725-747, 1993 Printed in Great Britain. All rights reserved

0045-7930/93 $6.00+ 0.00 Copyright © 1993 Pergamon Press Ltd

THEORY AND IMPLEMENTATION OF A SEMI-IMPLICIT FINITE ELEMENT METHOD FOR VISCOUS INCOMPRESSIBLE FLOW B. RAMASWAMY Department of Mechanical Engineering and Materials Science, Rice University, Houston, TX 77251-1892, U.S.A.

(Received 13 June 1991; in revised form 16 April 1992) Abstract--A numerical procedure for solving the time-dependent, incompressible Navier-Stokes equations is derived based on the operator-splitting technique. This operator split allows separate operations on each of the variable fields to enable pressure-velocity coupling. Discretizations of the equations are formed on a nonstaggered finite element mesh and the solutions are obtained in a time-marching fashion. Several benchmark problems, including a standing vortex problem, a lid-driven cavity and a flow around a rectangular cylinder, are studied to demonstrate the robustness and accuracy of the present algorithm.

1. I N T R O D U C T I O N

The solution of practical fluid flow problems by the finite element method has received increasing attention during recent years. In particular, the solution of incompressible fluid flow phenomena can now be performed efficiently in many practical cases. This potential originates in its ease in handling very complex geometries and the ability to naturally incorporate differential-type boundary conditions. Among the available combinations of the solution variables, the methods using primitive variables are considered important--primarily because these variables are more physical, have lower order equations and this form of equations provides a relatively straightforward extension to three dimensions. However, the primitive variable formulations have difficulties with the determination of pressure, which has no time term and is coupled implicitly with the divergence-free constraint on the velocity. This constraint, which is the continuity constraint equation, prohibits time integration of the incompressible flow equations in a straightforward manner. This is one of the major differences from compressible flow calculations. The most primitive formulation for the incompressible Navier-Stokes equations is constructed by the mixed formulation [1], where the velocity components and the pressure are used as variables. The coupled equation system of the incompressibility equation and the Navier-Stokes equations can be solved simultaneously [2] by a direct solver, such as LU decomposition or Gauss elimination. The use of the same interpolation for the velocity and pressure generates parasitic pressure. These elements do not satisfy the inf-sup condition of Babuska-Brezzi, as discussed in Ref. [3]. Therefore, various mixed finite element approximations have been developed recently. Usually, the order of the interpolation function for the velocity is chosen higher than that of the pressure to obtain physically reasonable solutions. For example, the multilinear interpolation function for the velocity and piecewise-constant one for the pressure are commonly used in the mixed formulation. However, this element occasionally yields spatial oscillation under particular boundary conditions. The penalty method [4-7] is a variation of the mixed formulation. In this method it is assumed that the pressure equals the negative of a penalty parameter times the continuity equation. This expression is then substituted into the pressure gradient terms in the Navier-Stokes equation and pressure is eliminated. The resulting equations are solved for the velocity components using the finite element method. The penalty parameter must be made large enough to enforce continuity, but not so large that the other terms in the Navier-Stokes equations are calculated inaccurately. The pressure is then recovered from the initial relation. In this approach, the velocity and pressure should also be approximated by the interpolation functions which are consistent with the Babuska-Brezzi condition. This problem can be resolved using a lower order integration for the terms associated with the penalty function [8]. It should be noted that the reduced integration 725

726

B. RAMASWAMY

technique for the penalty term leads to inconsistent imposition of the incompressibility constraint, since the incompressibility cannot be retained in the discretized system, as explained by Engelman et al. [9] On the other hand, low- and same-order elements are very attractive for their programming simplicity and the purpose of this paper is to develop a stable and efficient method using bilinear quadrilateral elements. The present effort is restricted to presenting a modification to the existing velocity correction algorithms [10-14] that reduces the restrictions on the time step with a simultaneous reduction in computing effort at each time step. In the following sections of this paper more details about the theory and numerical methods used in the program for laminar flow are given. We consider three numerical tests: the standing vortex problem [15]; the lid-driven cavity flow with different Reynolds number ranging from Re = 1 to 10,000; and flow past a rectangular cylinder at Re = 250. The main purpose of testing the standing vortex test problem is to see how much artificial diffusion will be introduced by the present scheme. The second test problem is interesting to the numerical analyst because of the non-linearity of the governing equations and the presence of high-pressure singularities near the upper corners. The third test example [16, 17] is to verify the features of the outflow boundary procedure. 2. GOVERNING EQUATIONS Consider a spatial domain f~ and a temporal domain (0, T) with x and t representing the coordinates associated with fl and (0, T). The governing equations consist of the following velocity-pressure formulation of the incompressible Navier-Stokes equations: du a--t"+ (u' V)u = V" a, on f~ x (0, T)

(I)

and V'u=0,

o n f ~ × ( 0 , T)

(2)

where u is the velocity and a is the total stress tensor (per unit density of the fluid), given as tr = - p I + 2vE(u)

(3)

E(U) = ½[Vu + (Vu)T].

(4)

with Here p is the pressure (per unit density), v is the kinematic viscosity and I is the identity tensor. For the purpose of describing the numerical formulation used, we assume that the boundary F of the domain f~ is composed of two parts: on F, the velocity vector is specified; and on F2 the stress "traction" vector is specified. To specify the velocity boundary conditions, we use the following equations along the boundary F, of the computational domain f~: u=b

onFl.

(5)

The normal and tangential components of the velocity boundary conditions can be written as n . U[r, = n . b

(6)

and n

x ulr, =

n

x b.

(7)

Natural boundary conditions are specified on the boundary/'2 as n.a=e

onF2,

(8)

where Ft and F2 are mutually complementary subsets of F, b represents the velocity vector prescribed on F,, e designates the traction vector prescribed on/'2 and n is the outward unit vector normal to F. The initial condition for the Navier-Stokes problem is u(x, 0) = u0(x) in ~.

(9)

A solution for viscous incompressible flow

727

The boundary and initial data must satisfy certain compatibility conditions if equations (1)-(9) are to represent a well-posed problem. In particular, we must have

V. Uo = 0

(lO)

nu0. n dF = 0.

(11)

and

When the velocity is prescribed on the entire boundary, i.e. when Ft = F and F2 = 0, the velocity boundary data b must satisfy

f

n . b d F = 0.

(12)

3. WEAK F O R M U L A T I O N To recast equations (1)-(9) into a functional setting, we introduce the following notations. Let Ht(t)) be the standard Sobolev space of scalar functions defined on f~,

Hl(~) = {u]u ~ L2(~), u' ~ L2(f~)},

(13)

where the prime denotes the first derivative with respect to the spatial variables. Similarly, let H~(f~)d be the Sobolev space of vector functions, Hl(~')) d = {Bin ~ L2(~"))d, u' u L2(~)d}.

(14)

In the spaces L 2(fl) and L2(f~) a of scalar and vector square-integrable functions, the scalar product is defined, in the normal way, by # (u, v) = |dlauv (15) 3

and

f (u, v) =

d

(u,,

dt)uv = ~

i=1

v,).

(16)

A weak formulation of problem (1)-(12) is obtained by multiplying equation (1) by a test function v ~ HI(r) d, VJr,= 0; and equation (2) by a test function q ~ L2(~). After integration by parts and noting that E(u) is symmetric and that equation (8) holds, the stress term gives

where we have defined

(V" a, v) = (p, V" v) - 2v((u, v)) + fdF2c • v,

(17)

d ((u, v)) = X (E.(u), c.(v)).

(18)

i,j=l

Thus, problem (1)-(12) can now be restated as follows: Given b, c and Uo, find u(x, t) ~ Ht(fl) u and p(x, t) ~ L2(£~), such that

- ~ , v + ((u. V)u, v) = ~ , V. v) - 2v((u, v)) +

f

drhc- v,

gqeL2(l)),

(V-u, q) = 0,

U]rl = b

V v ~ H ' ( t 2 y , Vlr~ = 0,

(19) (20) (21)

and

u(x, 0)

=

uo(x).

(22)

Under suitable smoothness conditions for the boundary data b and e, the Navier-Stokes problem (19)-(22) admits at least one solution [e.g. 18].

728

B. RAMASWAMY

4. O P E R A T O R - S P L I T T I N G SCHEME The time integration schemes can be classified into the explicit method and implicit method. The former has the advantage in computational storage, while the time step must be small enough to satisfy the stability condition. On the other hand, the latter is not restricted by the stability condition usually, but the inversions and storage of large matrices are required in the time integration process. Therefore, the most eligible time integration should be adopted for each problem, taking the characteristic feature of the problem into account. The ultimate goal of the present study is to construct an efficient and robust finite element scheme for incompressible fluid flow problems. Such a finite element scheme is assumed to be competitive with the finite difference method, with the two methods being more or less complementary to each other; the ideal numerical scheme would be a hybrid containing the advantageous properties of both methods. Therefore, the present study takes the approach of extracting the advanced numerical techniques from the finite difference method [19] and assimilating them with the finite element method. Emphasis is placed on the mesh subdivision, solver and algorithms used in the finite element computations. In this paper, the incompressible Navier-Stokes equations are decomposed into the hyperbolic, parabolic and elliptic equations and each equation is integrated fractionally by different suitable time integration schemes. The hyperbolic convective equation is temporally integrated using the two-step Adams-Bashforth scheme to save on computational effort. The fully implicit time integration is used for the incompressibility equation, because the fluid is always incompressible. The same time integration scheme is applied to the time integration of the pressure gradient term, since the pressure plays the role of the Lagrange multiplier for the incompressibility constraint from the viewpoint of the variational principle. Such a time integration yields the elliptic pressure Poisson equation, which is solved by direct solvers. The diffusion equations are integrated temporally via the implicit backward Euler scheme. The underlying philosophy of the present study is to reduce the computational costs as much as possible, daring to sacrifice the accuracy. The author believes the operator-splitting method is the most suitable for this purpose, not only because the algorithm can be constructed on semi-implicit and iterative bases, which naturally require less storage resources than fully implicit and direct bases, but also because the CPU time would be reduced by utilization of an appropriate time integration scheme for each equation. The resulting algorithm is given below. We will consider an operator-splitting scheme which is based on an orthogonal decomposition theorem due to Ladyzhenskaya [20]. This theorem states that any vector field which is the gradient of a scalar function is orthogonal to any solenoidal vector field with zero normal component on the boundary. By virtue of this theorem, the momentum equation and the incompressibility condition in the Navier-Stokes problem can be treated according to a procedure consisting of two distinct steps. Let us denote by u'~ and pm the velocity and pressure solutions of the discretized equations at time level t '~, where t m = t m - ~ + A t m and t o = 0. From u mand the boundary conditions the fields um+t and pm+~ at time t m+t are calculated as follows: Step

1

The calculation starts from u° = Uo, with V • Uo= 0. Then from um and boundary data b m+t =b(x, t "+t) and em+~= e(x, tm+~), find the intermediate velocity field U* ~ Hl(~'))d such that U*

--

Um )

.St-~.~~

,v

= - [~ ((u ~. V ) u ~, v) - ~((u~ -'. V ) u ~ - ~, v)]

-

2v ((u*,v)) + fdF2(c m + I+ p m + ~n)" v, V v e Hi(f~)a,Vlr,= 0, (23) (urn+l,v) = (P.u*, v), V v ~ H'(f~)d

(24)

U*lr, = IIe'+ i.

(25)

and

A solution for viscous incompressibleflow

729

Step 2 Finally the solution is advanced from u* to u"+~ via: find Hm+IEHI(~'~) d and p'+~ ~L2(f~) such that - ~-,-Zi

,v

)

VveH~(f~)U,n.vlc~=O,

V'v)=0, (V. um+~,q) = 0 ,

Vq ELZ(O)

(26) (27)

and n. u "+ lift= n- b "+ i.

(28)

In this step, only the normal component n •b" + ~ is imposed on the final velocity field ur,+ I. As a result, the boundary conditions for the tangential components of the velocity will be satisfied by u" +~ in a weak sense only. It can also be said that, in this formulation, the boundary conditions for the normal component are of an

essential type, whereas those for the tangential components are of a natural-like character. 5. F I N I T E E L E M E N T M E T H O D This section introduces a few of the basic terms used in finite element methodology. The finite element approach to discretizing a problem divides the domain of interest into elements determined by some number of nodes. In general, these dements can be of any shape. In this paper, all elements are quadrilaterals with four nodes. Quadrilaterals were chosen over triangles since quadrilateral grids typically have about as many element as nodes, while triangular grids typically have twice as many element as nodes. Since solution accuracy is primarily a function of the number of nodes, for a given accuracy a quadrilateral mesh will have fewer elements and will therefore use less memory than a comparable triangular mesh. When the finite element approximation is applied to equations (23)-(28) by the Galerkin method, it yields non-linear simultaneous ordinary differential equations of the following forms.

For Step 1 M ( u * - u")/At"

+1 =

_ [3K(um)um

__

½K(u"- J)um- 1] _ v S ( u * ) +

F m +'

(29)

and u*lrl = bin+ i,

(30)

where M is the consistent mass matrix, K is the advection matrix, S is the diffusion matrix resulting from the scalar product (E;j(u), Ei:(v)) and F results from the imposition of natural boundary conditions and the exertion of a body force.

For Step 2 The spatially discrete form of problem (26)-(27) is M ( u " +' -

u*)/At

m +1 _

Cpra +1 =

DU''+ ~= 0

0,

(31) (32)

and n . u "+' = n .

b "+l,

(33)

where C and D are pressure gradient and divergence matrices, respectively. Equations (31) and (32) represent a linear system of coupled equations for the unknowns um+~ and pro+~. The dimension of the system is equal to the number of degrees of freedom of the velocity and pressure altogether. The orthogonal decomposition theorem, which is the fundamental for the operator-splitting method, can be stated as follows: L 2(f~)a= H ~ H ±, (34) where H is the space H = {ulu ~ L2(f~)a, 17- u = 0, n. ulr = 0} (35)

730

B. RAMASWAMV

and H z is the orthogonal complement of H (in L2(~"~)d) and can be characterized as follows: H ~ = {uIue L2(f~) a, u = Vp, p ~ H~(f~)}. (36) in equation (34), the symbol • denotes an orthogonal sum. Relation (34) amounts to saying that um+l is the orthogonal projection of u* on the space H in L2(fl) d. Thus, the element u~+~ is well-defined by relation (24); we write U m + l = Pnu*, (37) where Pn denotes the orthogonal projection operator in L2(f~)a on the space H. The construction of this operator also suggests a convenient elimination scheme for solving the linear system (31)-(32) without storing its exceedingly large matrix of coefficients. Equation (31), solved with respect to u" + 1, gives um+~ =tl* + A t m + l M Icpm+l. (38) Then substituting u" + ~from equation (38) into equation (32), we obtain the following linear system which governs the pressure field: 1 ( D M - I C ) p m+l At.,+l Du* (39) The solution p " + t of equation (39) can be expressed as 1 p'~+ ~ = -Atm+ l ( D M - I C ) - t Du *.

Once p~+~ is obtained, the velocities um+* are obtained from equation (38).

0

: pm+l (unknown)

Q : P= (known) Fig. 1. Coordinatesystem for outflowboundary conditions.

(40)

A solution for viscous incompressible flow

6. OUTFLOW

BOUNDARY

CONDITION

731

FOR PRESSURE

In general for time-dependent problems, neither the pressure gradient nor the pressure itself on the open boundary problems can be prescribed because urn+’ on rJ (downstream boundary) is unknown. Here we follow the procedure successively adopted by Shimura and Kawahara [21] for the explicit velocity correction method. Here we calculate the pressure on rj by using a similar procedure to that for the inner domain of R. Consider the normal and tangential projections of the Navier-Stokes equations on r, where the coordinate system is shown in Fig. 1: WI+1 8P =vV'u;-(um.V)u::-(II;+'-u;)/At an

(41)

m+l ap -=vV2ur:-(u".V)u:-(u;+'-u;)/At a7

(42)

and

clef’ -+ an

ac+l --La7

=o.

(43)

Taking the normal and tangential derivatives of equations (38) and (39) and using equation (40), the boundary pressure Poisson equation is derived as

a2 a2

(

@+a22

)

p

m+,

2ww a7 an

,2auzw

an

a7

’

(44)

An additional Dirichlet boundary condition is required to solve equation (41) using the Galerkin method at each time step. Then, the nodes on which the Dirichlet boundary conditions were given will not vary as we move further in time. To avoid this inconvenience, the adjacent nodes to the open boundary are considered in an approximate sense. The integration of equation (41) can be considered as those over the adjacent to rj and can be written as follows: [Llp;:‘]=[;]?

(45)

where pm+’ on rs are unknown and pm on the adjacent nodes are known, L is the Laplace matrix and f mis the velocity gradient vector on the RHS of equation (44). In this approach, only the initial value of the pressure is necessary to solve the boundary pressure Poisson equation (45). 7. STABILITY

ANALYSIS

In this section, the stability criteria are analyzed for the two-dimensional advection-diffusion equation. The conditions for numerical stability of the equations resulting from the linear spatial discretization combined with the first-order Euler scheme as an explicit time integration technique are provided. It is generally known that the explicit method is conditionally stable and the time increment At has to satisfy the following condition: 0
au au dt+,u,dx+,ul%

g+“” ay

. w >

(

(46)

Using an explicit-Euler discretization for the time derivative term, forward difference for the advection terms and central difference for the diffusion terms, we obtain the difference equation: It+1 uij

-

u; = _

>

At

”

-v

Ui+lj_2U",+Ur_,j+U~+,-2U~+U~_,

AX2

(47)

B. RAMASWAMY

132

Assume an initial condition in the form of a spatial Fourier component as given below: uyj = a;, exp( - i2rrkxi)exp( - i2nlyj).

(48)

To determine the analytical relationship between the spatial frequency and the frequency of the nth component, equation (48) is introduced into equation (47) to give ai: ’ exp( - i2nkx,)exp( - i27tlyi) - at, exp( - i2nkxi)exp( - i2111yj) =At

--(u(

a$ (exp[ - i27ck(x, + Ax)] - exp[ - i2rrkxi])exp( - i27c1y,) Ax

[ _ ,v, a”,,exp( - i2nkx,){exp[ - i2rrf(yj + Ay)] - exp( - i27c1yj)) AY

~“,,{expt-i2~k(xi+Ax)]-2exp(-i2nkxi)+exp[-i2~k(xi-Ax)]}exp(-i2~~yj) Ax2 + a& exp( - i2nkxi){exp[ - i2d(yj + Ay)] - 2 exp( - i27c1yj)+ exp[- i2d(jj - Ay)]} )].

Ay2 In dividing both sides of equation (49) by a”,,exp( -i2nkx,)exp(-i2nlyj), rewritten in the following form: II+1 OkI -= a%

l+At

-

(u(

N

(49)

equation (49) can be

cos(27ck Ax) - 1 + i sin(2rrk Ax) Ax

cos(2~1 Ay) - 1 + i sin(2nl Ay) + 101 AY > _ 2v 1 - cos(2nk Ax) Ax’

(50)

The amplification factor for this problem is n+l akl

l-l

<

1.

(51)

a;!

For the stability, it is necessary to have - 1 < 1 - A$4 + 2vBI < 1,

(52)

where A=lul

cos(2nk Ax) - 1 + i sin(2nk Ax) cos(2nl Ay) - 1 + i sin(2nl Ay) + loI Ax AY

and B = 1 - cos(2nk Ax) + 1 - cos(27rl Ay) Ay’ ’ Ax2 The inequality equation (52) can be rearranged as below: 2 >/ AtJA + 2vBI > 0.

(53)

The RHS is always satisfied, so we just need to check the LHS value. For high Re flow, the advection terms dominate. The inequality equation can be written as -

2

IAI

> At.

(54)

Take Ax = Ay = h, k = 1, then ,A,

=

(,u,

+

,v,)

P - 2c0@nkhnkhl~.~ h

*

(55)

A solution for viscous incompressibleflow

733

Then taking the minimum value of the RHS of equation (55), we can get h (lul + Ivl) ~> At.

(56)

On the other hand, for low Re flow, the diffusion terms will dominate. The inequality equation can be written as 1 --/>

vlBI

At.

(57)

For simplicity, we take Ax = Ay = h, k = / , then Inl = 211 - cos(2rrkh)]h 2 .

(58)

Taking the minimum value of the RHS value of equation (58), we can obtain the criteria for stability as h2 4-~ ~>At. (59) Taking the characteristic velocity as (lul + Iv I), to be consistent with equation (56), then the stability condition becomes h 2 Re I> At, (60) 4h(lul + Ivl) where Re =

v

(lul + Ivl)h"

In our numerical calculations, based on the semi-implicit finite element method, we have chosen At 10 and 4 times bigger at the low- and high-Re ranges, respectively, when compared to the explicit scheme [12]. The semi-implicit scheme was found to greatly enhance stability and, in addition, provided the same steady-state solution at less computational cost than the equivalent explicit scheme. 8. N U M E R I C A L EXAMPLES

8.1. Standing vortex problem This problem was initiated by Gresho [15] at Lawrence Livermore National Laboratory and was extensively studied by Tezduyar and his research group [23] at the University of Minnesota using various velocity-pressure and vorticity-stream function formulations. The problem is governed by the assumption of an inviscid fluid flow in a square cavity. The initial condition (see Fig. 2) for this problem consists of zero radial velocity and with circumferential velocity is given by Uo= {5r for r < 0.2, 2 - 5r for 0.2 < r < 0.4, 0 for r > 0.4}. This initial condition was assumed by the above-mentioned investigators [15, 23] for this problem, The mesh system yields 961 nodes and 900 elements. The computed results are shown in Figs 3 and 4. At t = 3, the present scheme introduces an artificial diffusion by an amount of 2.5%. Table 1 illustrates, how much energy is preserved in the domain at different time level. 8.2. The shear-driven cavity The driven cavity flow is taken as the second application of the proposed scheme. Problem definition is given in Fig. 5. For the numerical analyst, it is a challenging numerical problem, because of the nonlinearity of the Navier-Stokes equations of motion, and because of the presence of regions in which the solution changes rapidly. From a purely computational viewpoint, the cavity flow is an ideal prototype nonlinear problem which is readily posed for numerical solution. Its geometric simplicity, comparatively minor singularities, well-defined boundary conditions and no preferred flow directions makes it very attractive as a test case for new numerical techniques. The reported calculations are for the asymptotic steady-state results in a square driven cavity at Re = 1,400, 5000 and 10,000. Two grid systems, 41 x 41 and 65 x 65 nonuniform grids with finer

734

B.R.AMASWAMY

..

.

.

.

.

.

.

.

g

IJ

~1 I oL~ I tN W[D

o?8

o

!

E

I",

<

I:: 0

0

l:l,

o

i"-

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

:,,,', ~~"~,,"x"~.,

.....

.,

:.,

~',,~~~.-..?;

,t' N ..... , ~! ~ , % , ~ ~ / / :

....'1 . . . . ,~II'fI

....

. . . . .

iiiiii!!!)i!!!!iiiiiiiiiii!i!

'..... ....

......

I,.::;

Z.~,....

glt."'

,

.......

Jolol^ Xl!OOlO^ (p)

lii'l~tlll'id I°U°lllUlUl!O £ (t)

O

1.0-3#!;£'0=dals Jnoluoo 'llnoluo 0 6 O0+3ttgL'O-=U!R 'i~O'-3t:O~"O=XOFt llnoluoo lunsuJd (:~)

Z 0 - ] 0 / - t ' 0 - d e l s Jno~uo3 'tlrlO)eu°~) 6 Z 0 - 3 9 9 g ' 0 - = u ! f l ' L0-':1~'1. L"0=XD~II ilnoluoo Xi! .O!lJOA (q)

tO-]t,~-t'o=dels Jnoluo~ 'lunoluo"J ~ i ~0-3~ |~'O-=U!M '00+3L6 I.'O--XDM UU!lUa~qS (o)

tll!o!po ^ lOUO!liulul!o £ (o)

'(~0'0 = IV) poqlatu qlJojtls~] stuvpv ~qi ~u!sn ~ = ~ li~ tu~lqoad XOI.tOA~'UtpUVIS ~ql JO suo!lnlo S "~ lt!.i

736

B. RAMASWAMY 1.2

1.01 O 0.8

== "-

0.6

._~

.E

0.4

Second order Adams-Bashforth 0.2

0.0

,

I

0

,

1

I

,

2

3

Time

Fig, 4. Time history of the vortex kinetic energy for the standing vortex problem. Table 1. Time history of the vortex kineticenergy Time K.E. (%)

Adams-Bashforth

0.5 0.996

1.0 0.991

1.5 0.987

2.0 0.983

2.5 0.979

3.0 0.975

grids clustered in the near-wall regions, were used. The cavity lid is always the upper surface, the lid is always impulsively started, and the lid always moves from the left to the right. The transient calculation is started from initial conditions of zero velocity and continued until a steady-state regime is reached. The criterion for convergence to steady-state was generally taken to be Ilu m÷~ - umll~ ~< 1.0 × 10 -5 , Ilum+'ll,

(61)

where Ilul[, is the L, norm obtained as the sum of the absolute value of all velocity components on the interior grid multiplied by the cell size. The reported computational results for the u=l,v=O

u=O,

~

v=O

•

u=0,

v=O

u=O,v=0 Fig. 5. Sketch of a shear-driven cavity and boundary conditions.

A solution for viscous incompressible flow

ttl

~ & t t { l I !!!!.,.~

[t11111

"" \ ' \ ' \ ~ \ \ ' % ' "

II!I1

[tlfllt llllkV liult~

~xxx.x.-..-..... _~

"-'~.~

\

,

k_ "~ II

o~

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

\\\,,~,

.......

ttI!I1tti f,',': -,~~ ~V~ ~, ,,,,, lHlll~',,,~,,__~-"-"/1

1 1 / / ; ;

~

l

Z

7

-'----.

I

1 i

b

/

,

...........

t l t l l lll{tt~ " . ( / ( ! I i I ~, ~,, i

~

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

filltlllll I ~ I . - _ ~ - - . ~ . ~ x . . . ~ \

i

.

, , ............

t11tlltlIt; ~ ~ - - - % . R ~

t

~

t ~ t i i i ff.q,~ ~ ~ ~ ~. ~. -:~,~ ~11

tiltlllt~~"'-~"-~.-. i•__k• tll

e~

( 0

~ " 1 1 Z / ! / II l ~.'":] i / i Z Z Z Z II ~ '""\:I
It.

l

II

...... I

/ I (£ I'F~'lllll]:::[[[ I

i~\%'<\H-~

~

~ / / ~

" X V \'\'\'\'\\~'"'""1 ~l p I'pl'ytiiil.''"l

ttlll~\~ \\\\\-~-_

737

..........

g

"~

........ ~

o , , , - "

...... .~"

s~

ltlllllTI;l~g~_ P " \ \ \ \ \ \ ~ "'' t1111tlII;(;~7~_-'.".""- ~ \ \ \ \ ~ ~ ' ' I1111111t1/I, - " \ \ \ \ \ \ \ ~~~"

............ ............. ............

g

Itl,li(~_'"~ IIk%2-"//i I\"%2_-~<
l~ ~ , , , , / / ~ ~''"

2 ," / / ~ , , , ,

.~

........... ..........

...........

._~ 0"

0

II

' t¢,,~¢¢¢

•

• ,~ /

.... m / , ~ i f ; ".,,it}. kit////

',.at/l/i,

kDIt

.,,~I~tI7777

,~tlll """;fJ,t,f,f,f

~.~-~

/ / //,.

,,~lllltA&k\\ .,,~\\\\\\~ x;.,\ X\\ x\ , . ~_ . ,

" > ,

\ \\\\\,~tt~l ,' , ,*, i

I

~l ttlli~tlll

......~\ \ \ \ % ;

~\"-t---t-~-~-m~Dl

........, ~ , ~ I~IIIttl

#

I

I !

',S,Iiltlltl

!

tiiltl

t>"~

l,J

' , ~

zT,

"0

,-i

o

O

;u'~vA~svvtv~I "1]

~t7/

A

solution for viscous incompressible flow

739

Table 2. Values of the pressure contours for the shear-driven cavity flow Re Contour No. 1 100 400 1000 5000 10000 1 --9.6000 -0.0800 --0.1000 -0.1000 -0.1000 -0.0800 2 3 4 5 6 7 8 9

-6.4000 -4.0000 - 1.5000 0.0000 1.5000

-0.0600 -0.0400 --0.0200 0.0000 0.0200

-0.0800 -0.0600 -0.0400

4.0000

--

--

--

6.4000 9.6000

---

---

---

-0.0200 0.0000

-0.0800 -0.0600 -0.0400 -0.0200 0.0000

-0.0800 -0.0600 -0.0400 -0.0200 0.0000 0.0200 0.0400 0.0600

-0.0600 -0.0400 -0.0200 0.0000 0.0200 0.0400 0.0600 --

asymptotic steady-state cases are compared with the published data in Ref. [22]. The asymptotic results for the square driven cavity are given in Figs 6--11 for Re = 1,400, I000, 5000 and 10,000. The 65 x 65 grid resolution is too coarse for the resolution of any additional vortices in the corner vortex cascades. Pressure contour values are given in Table 2. The profiles of the dimensionless horizontal velocity along the vertical centerline of the driven cavity (x = 0.5) and the dimensionless vertical velocity along the horizontal centerline of the driven cavity (y = 0.5) are illustrated for various Re in Fig. 12. The results from Ghia et al. [22] for a grid with 128 x 128 cells are overlayed on these plots. Very good agreement can be found for both the maximum velocity position and value for Re values up to 10,000. In particular, it must be pointed out that with the present method there is no need for upwinding at high Re. Finally, in Table 3 we compare the properties of the primary and secondary vortices with those given in Ref. [22]; the agreement should be considered satisfactory. In this table, we also provide the time increment At and CPU time required for each Re. LO

/

06

g >

0.2

---e---- S e m i - i m p l i c i t

~

oih G /

I

(41"41)

et0~(129"129) •

{

-0.2

8 -0.6

-I 0 10

i •

n

,

I

"

'

.

'

-

~

Re = 10,000 I Semi- I mpliclt (65" 65) • Ghia

~11 J r

etoZ(129*1291~l~

06

0.2

~

-0.2

8 -0.6

-1.0

n

L

-7.0

-0.6

-02

0.2

06

1.0 -1.0

-0.6

x Coordinote - u Velocity

Fig. 12. Velocity profiles for the shear-driven cavity: (a) R e = 4 0 0 ; (d) R e = 10,000.

-0.2

0.2

06

x Coordinote - u Velocity

(b) R e = 1000; (c) R e = 5 0 0 0 ;

1.0

B. RAMASWAMY

740

Table 3. Stream function extrema, C P U and At for the shear-driven cavity flow

Re

Bottom left 1 vortex (large)

Bottom left 2 vortex (small)

Bottom right 1 vortex (large)

Bottom right 2 vortex (small)

Method

Primary vortex

T o p left vortex

Semi-implicit

-0.103

--

7.15 x 10 -6 1.17 x 10 -5

-__

Ghia et al. [22]

-0.103

--

--

Semi-implicit

-0.114

--

Ghia et al. [22]

-0.114

--

Semi-implicit

-0.118

--

1.75 X 10 6 1.25 x 10 5 1.63 x 10 5 6.77 x 10 4 1.42 x 10 ~ 6.42 x l0 4 2.25 x 10 -4

100

400

1 . 7 7 × 10

1000

Ghia et al. [22]

-0.118

--

Semi-implicit

-0.117

1.24 x 10 -3

5000 Ghia et al. [22]

-0.119

1.46 x 10 -3

2.31 1.75 1.32 2.99 1.36 3.08

x × x x x x

CPU' 276.6

At 0.040

--

-__ --

553.1

0.040

829.7

0.040

2765.6

0.020

--

3

--4.80 x 10 -7

10 4 10 -3 10 3 10 3 10 3 10 -3

__ -9.32 x --3.80 X -7.09 x -1.43 x

10 s 10 -7

10 -s 10 -6

__

I

=CPU time in seconds using a C R A Y X-MP.

8.3. Vortex shedding behind a rectangular cylinder In this subsection numerical results for ftow past a rectangular cylinder at Re = 250 are presented in order to verify the features of the proposed procedure. The calculated flow fields are compared with experimentally observed ones, which are reported by Davis and Moore [16]. The finite element mesh and the boundary conditions are shown in Fig. 13. The mesh system yields 3010 nodes and 2880 elements. At the inlet boundary, a uniform velocity distribution is assumed. To solve Poisson's equation for pressure, the outflow boundary condition as described in Section 6 is used successfully. (a)

u : l , v:O o ~ !

u=O, v=O on cylinder surface

r~

(j C>

b

~3

:E

C>

"5 C>

u=l, v=O

,_ _I I. 4,5H

(b)

H

i | l l l l l l ilinIHili iUimUnll

ii..

14,5H

~ I ~ ~ llqlllllll|linUillilllmUllilllUlilalll !ll1111111illlilDlllmilliliOlilliml

iimRiimI

i

Il'~:',:ll-'-'ll=--'||'--=lll-'-=-'-=l-'==-=:=-'-=-lliil

.....

i i g i n g | i

"'"'iii

ii!iii!ili!iiF!'irrEriiiiiiiiiiii|

i i l D i i l l l RBHHBlil

t hllillilllliBninlmiliill|llRililHilI i I mlti||iniiiiI iIilniIiiIiiI

lilHImi nil | | ililliHlll f a l l i B l e

~ Illlt I I I I I l l g i I i IllnmllIlII oI IlllliUllililillilmilmilnllilIRmHiH ;~t l I l l ~ U l l l l l i l l U i l l i i l l l l l l l l l i l i l i i n l

"" -----.,, ......,, ...... l i i m i I l l l

Iil I r l l l l l i l l l i l l l i l i l l I i l i l l i l i I I i l i I i m l l i

if

llillIIlIli

',IIIpIII:IIIIIIUIIIlIIIlIIIIIlIIIlI,

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

i H*IIIlIII!ilI!'II'-'-II--'!I!I|'!|'I!'III

........--.....,-----.,,--' ..... "'"" .l l. i . l .l .i l. l . l . .

IIlilIiIinii

:;! .......................................

I

!1 I i

!iii ii!iiiliiiililiiiililiiliililii

i;II'~IIIilIIIIIIIIIIIIIIIIIIIUlIIUnUU

Fig. 13. F l o w p a s t a r e c t a n g u l a r c y l i n d e r : (a) p r o b l e m definition; (b) finite e l e m e n t m e s h .

A solution for viscous incompressible flow

741

.__L___ _

_

8 e~

O

.g

II

b ~0

._=

O

E

e~ O

11111111~11111~I111111 I 11111111~IIIII~ITITIII I I v

~r

.Q

<~>

(e)

.

.

.

.

.

.

.

.

~. . . . . . .

~::

-

.

i

.

.

.

.

.

.

:

.

.

.

.

.

.

.

.

.

.

.

.

~-

.

_~.-~..-.._.,,...~.-.._....

k,,.__

~

C L =

_

--~.,,r-,,,---

max.

:

~

~._...__.~...-~..--..--.

Fig. 15. Velocity vectors and dimensionless pressure contours during regularly periodic sheddings at Re = 250, At =0.015,

'[/jt~~

.

.

........ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

(b)

(8)

_

_

.

_

_

I

.--.---.~__..:...._...~....~_....~.

Fig. 16. Velocity vectors (a) and dimensionless pressure contours (b) during regularly periodic sheddings at Re = 250, At = 0.015, CL = 0.

_

0

ct~

8

0

¢~

>

(b)

.

.

.

.

.

.

.

.

-------I

~

~

~

- -

-

~

.

.

.

.

.

.

.

.

.

~

.

i~,--

i

.

-

-

.

.

.

.

-:

.

.

_

.

.

.

.

.

.

.

.

.

.

~

~-~.~,,,',i'~'i.,,.,,,......a,~.~-.~

~

.

.

-

~

-

,

.

_

.

-

.

.

_ :

.

-

.

.

.

..........

~

_..

=

•

,

_

_-......

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

~..-.

~

~

...:..~1..i3~

~

~

. ' , , .'i". . ' ;

~ .-...

: ~ . - . ~ ' . . . - - ~ ~ . ~ . . . . _ ~

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

.... ~

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

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

~-'.,.-..-...',,.--.._

~

Fig. 17. Velocity vectors (a) and dimensionless pressure contours (b) during regularly periodic sheddings at Re = 250, A~ = 0.015, CL = rain.

-

. . . . . . . . . . . . . . . . .

.

.

Z,

>

A s o l u t i o n f o r viscous i n c o m p r e s s i b l e flow

745

Table 4. Values of the pressure contours for vortex shedding behind a rectangular cylinder Contour N

Symmetric vortex (Fig. 14)

C t = max (Fig. 15)

CL = 0 (Fig. 16)

C L = rain (Fig. 17)

I 2 3 4 5 6 7 8 9 10

79.4500 79,6000 79.7500 79.8250 79.9000 79.9750 80.0500 80.1250 80.2000 --

113.375 113.600 113.675 113.750 113.825 113.900 113.975 114.050 114.125 114.200

111.950 112.025 112.100 112.175 112.250 112.325 112.400 112.475 112.625 --

111.450 111.525 111.675 111.750 111.900 111.975 112.050 112.125 112.200 112.275

Along the other two boundaries, tangential tractions and the normal velocity components are assumed to be zero. The nonslip condition is prescribed at the cylinder surfaces. Note that the velocity vector and pressure contour plots given in Figs 14-17 do not display the whole flow field domain defined by the mesh system. The partial domain plotted is from 1H to 19H in the x-direction (instead of 0 to 20H for the full domain) and 3.7H to 11.3H in the y-direction (instead of 0 to 15H for the full domain). The nondimensional freestream velocity and Re chosen for the case investigated in this study are 1 and 250, respectively. The time step is At = 0.015 s. Figure 14 shows the velocity vector and pressure contour plots for the flow field at 30 s. The flow has not developed fully yet, and the flow appears symmetric, with twin vortices forming behind the cylinder. As time progresses, the flow becomes asymmetric and periodic vortex shedding over the cylinder is attained. The velocity vectors and pressure contours during a half cycle of shedding are shown in Figs 15-17. Pressure contour 0.6

1.74

0.4

1.72

1.70 0.2

1.68

0.0

~

-0.2

1.66 1.64 1.62

-0.4

1.60

-0.6

1.58 162

166

170

174

178

' 162

Time

' 166

170

174

178

Time

(4) 0.6 0.4

0.2 0.0 -0.2 -0.4

s-

•

s -

-0,6

X

×

/,

-0.8

i

162

i

166

i

¢-,

j,-~

/, i

170 Time

i

...~ i

174

i

178 (b)

Fig. 18. T i m e histories o f (a) the d r a g a n d the lift a n d (b) the s u r f a c e pressures.

746

B. RAMASWAMY Table 5. Some global quantities associated with vortex shedding Source

Shedding period

Wavelength

Drag coefficient

Lift coefficient

Davis and Moore [16] Semi-implicit

5.9-6.3 ~ 5.73

5.9 5.7

1.73-1.81 1.59-1.73

-0.5~-0.50 -0.51 0.51

~Numerical (several meshes) and experimental,

values are given in Table 4. Figure 17 shows the results at 187.5 s for a fully developed flow where the period o f vortex shedding is constant. The lift coefficient, drag coefficient and the nondimensional pressure values at the center of each cylinder side are plotted in Fig. 18 for t ranging from 162 to 178 s. The CL is the difference between the pressures on sides $2 and $4, and CD is correlated with the variation in the pressure on the back side $3. N o t e that the oscillation frequency o f CD is twice that o f Q . At CL = max or CL = rain, Figs 15 and 17 show that the flow separates at the leading edge, and the vortex is detached from the cylinder. At CL = 0, Fig. 16 shows that the growing vortex is ready to depart from the cylinder. The resulting shedding period is 5.73 s, and the corresponding Strouhal frequency is 0.174. This is agreeable with the experimental results obtained by Davis and M o o r e [16]. Table 5 shows the comparison o f the c o m p u t e d results with the experimental observations made by Davis and M o o r e [16], 9. C O N C L U D I N G

REMARKS

The present paper is devoted to describing a solution procedure o f the time-dependent, incompressible finite element Navier-Stokes equations o f the primitive variable formulation. The "segregated" a p p r o a c h is followed and the equations o f m o t i o n are considered sequentially. The fundamental concepts and characteristics o f the formulations and the solution m e t h o d o l o g y used are described in technical detail, The semi-implicit scheme proposed in this paper uses equal orders o f interpolation for the velocity and pressure and this eases the implementation effort. The algorithm has been shown to have g o o d convergence properties in attaining steady-state solutions at an acceptable practical level o f time step. Results obtained are in g o o d qualitative and quantitative agreement with those published previously. Our future efforts will be directed at employing the present semi-implicit scheme to transient three-dimensional problems in conjunction with iterative schemes (such as the element-by-element method) and its parallel implementation. Acknowledgements---This research was carried out with the support of the NSF, under Grant DMS-9112847. The

computation were performed on a CRAY Y-MP at the University of Illinois, Urbana-Champaign, under Grant NSF TRA910175. The author would also like to thank Dr T. C. Jue for his many helpful discussions. REFERENCES I. R. L. Sani, P. M. Gresho, R. L. Lee and D. F. Griffiths, The cause and cure (?) of the spurious pressure generated by certain FEM solutions of the incompressible Navier-Stokes equations: Part 1. Int, J. Numer. Meth. Fluids 1, 17 (1981). 2. C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using FEM technique. Computers Fluids 1, 73 (1973). 3. F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multiplier. R.A.LR.O. 8, 129 (1974). 4. M. Bercovier and M. Engelman, A finite element for the numerical solution of viscous incompressible flows. J. Comput. Phys. 30, 181 (1979). 5. T. J. R. Hughes, W. K. Lui and A. Brooks, Finite element analysis of incompressible viscous flows by the penalty function formulation. J. Comput. Phys. 30, 1 (1979). 6. J. N. Reddy, On penalty function methods in the finite element analysis of flow problems. Int. J. Numer. Meth. Fluids 2, 152 (1982). 7. J. N. Reddy, Penalty-finite-element analysis of 3-D Navier-Stokes equation. Comput. Meth. Appl. Mech. Engng 35, 87 (1982). 8. D. S. Malkus and T. J. R. Hughes, Mixed finite element methods-reduced and selective integration technique: a unification of concept. Comput. Meth. AppL Mech. Engng 15, 63 (1978). 9. M. S. Engelman, R. L. Sani, P. M. Gresho and M. Bercovier, Consistent vs reduced integration penalty methods for incompressible media using several old and new elements. Int. J. Numer. Meth. Fluids 2, 25 (1982). 10. J. Donea, S. Giuliani, H. Laval and L. Quartapelle, Finite element solution of the unsteady Navier-Stokes equations by a fractional step method. Comput. Meth. Appl, Mech. Engng 30, 53 (1982). 11. P. M. Gresho, S. Chan, C. Upson and R. L. Lee, A modified finite element method for solving the time-dependent, incompressible Navier-Stokes equation. Int. J. Numer. Meth. Fluids 14, 557 (1984). 12. B. Ramaswamy, Finite element solution for advection and natural convection flows. Computers Fluids 16, 349 (1988).

A solution for viscous incompressible flow

747

13. B. Ramaswamy, Efficient finite element method for two-dimensional fluid flow and heat transfer problems. Numer. Heat Transfer 17, 123 (1990). 14. B. Ramaswamy, T. C. Jue and J. E. Akin, Finite element analysis of unsteady two-dimensional Navier-Stokes equations. Numer. Heat Transfer 21, 147 (1992). 15. P. M. Gresho and S. T. Chan, Semi-consistent mass matrix techniques for solving the incompressible Navier-Stokes equations. Lawrence Livermore National Lab. Preprint UCRL-99503 (1988). 16. R. W. Davis and E. F. Moore, A numerical study of vortex shedding from rectangles. J. Fluid Mech. 116, 475 (1982). 17. Y. Yoshida and T. Nomura, A transient solution method for the finite element incompressible Navier-Stokes equations. Int. J. Numer. Meth. Fluids 5, 873 (1985). 18. R. Temam, On the Theory and Numerical Analysis of the Navier-Stokes Equations. North-Holland, Amsterdam (1977). 19. A. J. Chorin, Numerical solution of the Navier-Stokes equations. Maths Comput. 22, 745 (1968). 20. O. A. Ladyzhenskaya, Mathematical Problems in the Dynamics of an Incompressible Flow. Fizmatgiz, Moscow (1961). [English Translation, Gordon & Breach, New York, MR27 # 50349b (1963).] 21. M. Shimura and M. Kawahara, A new approach for pressure boundary conditions using the explicit velocity correction method. In Computational Methods in Flow Analysis (edited by H. Niki and M. Kawahara), Vol. I, pp. 446-453 (1988). 22. U. Ghia, K. N. Ghia and C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and the multigrid method. J. Comput. Phys. 48, 387 (1982). 23. T. E. Tezduyar, J. Liou and D. K. Ganjoo, Incompressible flow computations based on the vorticity-stream function and velocity-pressure formulations. Computers Structures 35, 445 (1990).

C A F 22/~--E