- Email: [email protected]
Analysis of unsteady compressible boundary layer flow via finite elements
Computers and F/aids. 1976. Vol 4. pp 1-12 Pergamon Press. Printed in Grcal Rntam
ANALYSIS OF UNSTEADY COMPRESSIBLE BOUNDARY LAYER FLOW VIA FINITE ELEMENTS T. J. CHUNG'I"and J. N. Cmou:l: Department of Mechanical Engineering, The University of Alabama in Huntsville, Huntsville, AL 35807, U.S.A.
(Received 15 June 1975) Abstract--This paper is concerned with the discrete formulation and numerical solution of unsteady compressible boundary layer flows using the Galerkin-finite element method. Linear interpolation functions for the velocity, density, temperature and pressure are used in the momentum equation and equations of continuity, energy and state. The coupled nonlinear finite element equations are approximated by a third order Taylor series expansion as temporal operator to integrate in time with Newton-Raphson type iterations performed until convergence within each time step. As an example, a boundary layer problem of a perfect gas behind a normal shock wave is solved. A comparison of the results with those by other method indicates a favorable agreement.
NOMENCLATURE V, P T a~ tu,
F, t, n, E q, h d,j c~ c,, t
p aN, ,SN,3',~,I/',,, A, ~t K
velocity pressure temperature t~, = acceleration Cauchy stress tensor body force per unit mass surface pressure unit vector normal to surface internal energy density per unit mass heat flux per unit surface area heat supply per unit mass deformation rate tensor constant volume specific heat constant pressure specific heat time mass density interpolation functions for velocity, density, pressure, and temperature, respectively viscosity constants thermal conductivity
Indices--upper case and lower case represent element node numbers and spatial dimension numbers, respectively, unless otherwise noted. INTRODUCTION
Supersonic aerodynamics in space technology prompted an interest in compressible, nonsteady boundary layer flows. Boundary layers of this kind may arise behind shock waves and expansion trains in shock tubes, or in a fast-moving body subjected to deceleration or acceleration in its flight and to aerodynamic heating. In this connection, the analytical or numerical solution of unsteady compressible boundary layer flows still pose a formidable task, particularly for complicated geometries and boundary conditions. In recent years, considerable efforts have been made toward application of the finite element method to solve fluid mechanics problems, as evidenced by the recent international symposium on this subject[l]. In the present study, the nonlinear Navier-Stokes differential equations governing the motion are cast in a Galerkin integral approximation, leading to a coupled set of nonlinear integrodifferential finite element equations. The results are similar to the approximation for the NavierStokes equations as obtained by a direct method or the so-called virtual velocity formulation[2]. ?Professor of Mechanical Engineering. tGraduate Student. CAF Vo~ 4. No $--A
1
T.
J. ('l.tlN(i and J. N. CHIOI'
The energy equation written in terms of internal energy and equations of continuity and state are then combined together with the equation of motion for simultaneous solution, integrating across isoparametric finite elements[3] in spatial domain and using suitable temporal operators in time domain. A stable convergent solution is obtained by the temporal operator with Taylor series terms of the third degree in time and Newton-Raphson type iterations in spatial domain within each integration time step. To demonstrate the effectiveness of the present approach, a boundary layer problem of a perfect gas behind a shock wave is solved and compared with the results of others 14, 51.
INTERPOLATION FUNCTIONS AND WEIGHTING FUNCTIONS To derive the finite element analogues for equations of momentum, continuity, energy and state, it is necessary to introduce interpolation functions for the velocity V, density p, pressure P, and temperature T, in the form V,(X, t) = a,~(X)VY(t) p(X, t) =/3N(X)pN(t)
Ilb)
P(X, t) = yN(X)PN(t)
(lc)
T(X, t) = tON(X)TN(t)
(ld)
(la)
where V,N, pN, p~, and T N are the nodal values of velocity, density, pressure, and temperature, which are functions of time only. Note that the interpolation functions aN,/3~,, yN and qJN are dependent on the spatial coordinates X but free from dependence on time. Here i = 1, 2 for two-dimensional flow, N denotes the element node, and repeated indices imply summing. In the Galerkin's integral approximations, these interpolation functions are taken as weighting functions which require the basic convergence criteria of completeness and continuity[2, 3]. In the method of weighted residuals, we construct an inner product of the residual and the weighting functions. Often, the residual consists of more than one variable with each variable assumed to vary spatially with polynomials of different degrees. For example, in the momentum equation of a compressible viscous flow, the variables are the density, pressure and temperature, as well as the velocity. A question arises as to an appropriate weighting function for the momentum equation. The general rule is to seek one of the variables in the differential equation to construct an inner product of the residual of the differential equation and this variable so that the inner product will yield a spatial invariant. In this case, the obvious variable is the velocity, and we obtain an invariant I, = ( ~ , Vk) = (~k, aN(X)V~N(t)). Therefore, the weighting function is the velocity interpolation function aN(X). Note that the density p, pressure P, and temperature T are disqualified because they are scalars and do not possess a free index k required to yield an invariant as a result of an inner product with the residual ~k of the momentum equation. Likewise, for the continuity equation, we form an invariant as an inner product of the residual and the density. Thus I2 = (~, p) = (~, flN(X)pN(t)). Obviously the density interpolation function/3~(X) is the appropriate weighting function. For incompressible fluids, however, the density is constant and v,., = 4 .
Here the only variable involved is the velocity. However, we cannot choose the velocity as a subspace since the inner product (~, V~) does not yield an invariant. In this case, incompressibil-
Analysisof unsteadycompressibleboundarylayerflowviafiniteelements
3
ity requires that Is = (~, P) = (~, Tr~(X)PN). Consequently, the pressure interpolation function 3,N(X) is the chosen weighting function. For the energy equation, we find that temperature, pressure, density, and velocity are all involved as variables. The residual is denoted ~ and all variables except the velocity appear to provide an invariant from the inner product. In the absence of dissipative energy due to velocity gradients, the spatial distribution of energy is contributed by the gradient of heat flux or V~T. Thus it is apparent that an invariant is constructed as I3 = (~, T) = (~, #~(X)T N) and the temperature interpolation 0N(X) is the weighting function. Lastly, the equation of state for a perfect gas assumes an invariant from the inner product with the pressure. Therefore, the pressure interpolation function 7N(X) is taken as the weighting function. In general, it may be argued that, for a very large velocity gradient, higher order approximations for the velocity should be used. In the present study, however, a linear spatial variation of all the variables is assumed for simplicity. In particular, two-dimensional four node linear isoparametric elements are used for the velocity, density, pressure, and temperature. Consequently, the weighting functions for all four equations of momentum, continuity, energy and state are of the polynomial of the first degree. MOMENTUM EQUATION A linear momentum equation is given by pa~
-
pF~ - tl~., = 0
(2)
where ah = (DVk/Dt) = (:~ + V, Vk., is the acceleration, F~ is body force per unit mass, t,k is the Cauchy stress tensor, and the comma denotes a partial derivative. Let ~ be the residual of (2). We then have
3k = pa~ - pFk
-
t,k., = p(f'k + V, Vk.,)- pFk
-
t,~.,.
Using Galerkin's approximation with a~, as weighting functions, we require
f ~lkaN dv = O. This leads to f [p(~Yk+ VrVk.,)-pF~ - ta.,las dv = O.
(3)
For Stokesian fluid, the stress tensor assumes the form tt~ = - P & k
+ Ad,,&~
+ 21~d,k
where, for Stokes hypothesis, the viscosity constants/z and A are related by X = -2/3~., and the deformation rate tensor da is d,k =~(I V ~,k + Vk.,)
Noting that
f t,,.,a,,, dv = f t,~,,,,n, dA - f t,,.a,,,., dv = f t~,,,, dA - f t.,a,,,., dv
T.J.
4
CHt?NG a n d J. N . CHIOU
in which dA is the boundary surface, tk = t,~nk represents the components of boundary stresses and 8,~ is the interpolation function for the boundary surface pressure, N being only the boundary nodes involved instead of all nixies of an element. Assuming the viscosity constants/x and A to be a function of temperature such that #, = aT and A = bT where a and b are constants, we can rewrite (3) in the form, A~,N(/k M + BbN~ V, ~' + Z~tNkP ea = XN~
l~l, 1) ~ XNk
~
(4)
I~(2) ~_ L"(3)
~ N k "T- r" N k T
& " Nk
and A~,N = f fll.aMaN dvp L B ~Nk = C ~Nk + D'MNk C',aNk = f bez.OtM.,CtN.k dv T t" D ~ k = ~ aqSLA'Maa~.~ dv T ~ ZMNk = -- f 7~,aN.k d v
(1)-- f ~.F~a~ dvp R FNk-(2) F Nk -
t~tN dA
1
FN - --f fl~a~c~Lao., dvpMV, LVk ° (3)k -with A ~,~k= a~.k&" + aA,.~&" Here, AMN is the mass matrix, B ~,M~is the viscosity matrix, Z~,N~ is the dilation matrix, F "~ is the body force vector, F ~2~is the surface pressure vector, and F~k is the convective term hereafter called the inertia pseudo force vector.
EQUATION OF CONTINUITY
The equation of continuity is written in the form, + (p Vk ).~ =0. Again, using Galerkin's approximation residual + (pV~).k =
with weighting functions, tiM, we obtain
f
~#M
dv = 0
or
f[,o
"~-(p Vk ),k ]flM
dv =0.
Substituting the interpolation functions and performing the differentiations, we obtain GMNp N + N~Np N = 0
(5)
Analysisof unsteady compressibleboundary layer flowvia finite elements
5
where f GMN= J ~M~N dv HM~ = f ~SM(~NaL.~+ ~8~.~aL)dv VkL. ENERGY EQUATION The first law of thermodynamics requires the balance of energy as well as the conservation of linear momentum. The energy equation is given by De
p - ~ - t,d, - q,., - oh = 0
(6)
where • is internal energy per unit mass, q~ is the heat flux per unit surface area, h is the heat supply per unit mass. The residual function of the energy equation is D~
= p - ~ - tljd~j - qu - ph = p~ + pV~.l - t,do - ql,~ - ph.
Once again, using the Galerkin's approximation with weighting functions ~ , , we have f ~o'~ dv
0
or
f(p~ +
-
- q,., - ph )~bN dv = 0.
(7)
Here, for a perfect gas, = c.T - Pip
(8)
where cp is the specific heat at constant pressure. Making use of the Fourier law, qj = KT.~ where ~ is the thermal conductivity, and in view of the continuity equation, we may rewrite (7) in the form f (pcpi" + p % V~T., - P - A Vk.~ - 21~d, do - V~P., - q,., - ph )~bN dv = O.
(9)
From the Green-Gauss theorem it is possible to write
0o) where ~ = q,n~ is the input heat flux normal to the surface. Substituting (1) and (10) into (9) and simplifying yield CMNi "M + DMNT M + KMNP M + SMNP M = YN
with
and CMn = f Cp[3LOM~bn d v p L DMn = f KOM.I~bN.I dv
(11)
T. J. CHUNG and J. N. CHIOt"
KuN = - f
yuON dv
S~N = - f yM.iaLtPN dv
Vil"
QN"'= .//3,~h~Ndrp M Q ,2,= f (1~. dA QNTM = f b~uOtL,kfltQdd.INdv TMvkLvi Q =
2a~bu~NAL~Ap.dvVk Vt T
Q~'" = f - C,JSLapd/U.,OZ~dvp L V/.T u. For a perfect gas the internal energy is defined alternatively as • = c~T
(12)
where c~ is the specific heat at constant volume. Using (12) instead of (8), one obtains the finite element energy equation of the form (11) without the term containing the pressure rate. In this case, eqn (ll) is modified to CMt~7"M + DMNT u + LuNP u = YN
(13)
where
CuN = f C.[3L~bM$Ndvp L Lu~ = f yuaL.i~b~ dv V,L QN'" = f - CdJLaotPu.,ON dvp L VIOTu with other quantities remaining the same as before. E Q U A T I O N OF S T A T E
The equation of state which will now complete the description of compressible thermal flow is given by
P = pRT
(14)
where R is the gas constant. Proceeding as before, but this time with weighting functions y~,, we obtain
f (P - pRT)yN dv = 0 or
Eu~P u = fN where EuN = f y u y ~ dr;
fN = R f flL~Uy~ d v p L T u.
(15)
Analysisof unsteadycompressibleboundarylayerflowviafiniteelements
7
SOLUTION PROCEDURE The local finite element equations for momentum, continuity, energy, and state assume the global forms upon assembly[2, 3], respectively, A~,~'~ ~ + B ' ~ V 7
+ Z.,~,P ° = X~,
(16) (17)
God~" +H°op ° = 0 CoaT ° + D o l T ~" + K,,¢~ ° + Jo~P*" = Yo E,,j~P*=fa
(18)
(19)
with a,/3 = 1, 2. . . . . & ~ being the total number of nodes in the entire domain under investigation. Note that integration of all element coefficient matrices must be carried out before assembly. A computerized Gaussian quadrature may be easily applied for the integration. The time derivatives may be approximated by various difference operators. In the present study, the 3rd order Taylor series expansion for the temporal operator was used. The first time derivatives for velocity, density, temperature, and pressure may be expanded in terms of histories of the past time steps. For example, if g represents any one of the unknown vectors, then we have[6] 1
lg, = ~-~-(1 lg, - 18g,_, + 9g,-2 - 2g,_,)
(20)
where At is the time increment at any time step i. Thus, the incremental equations in recursive matrix form for momentum, continuity, and energy become 1 IA + B'~V ] i+ZP,=X.
IlG+
x
+ OAt .A (18V,_1 - 9V,_2 + 2V,-,)
G
H)p, =-~(18p.-9p,_:
+ 2p.)
(~-~t + D)T, + (~-~t + S)P, = Y,-I + £(18"1',-i - 9T,-2 + 2T,_,) + ~At ( 1 8 P . - 9Pi-2 + 2P~-s). To these the equation of state must be added. Writing all governing equations in a matrix form, we obtain N 0 0
0 R 0
0p= ST
Ii ° ° !llilIll
where IIA M=6--~'+B, c = X.
I1G_ N=~--~-+H,
_ 11C IlK R=6--~-+D , S = 6 - ~ - + J
A
+ ~-~'(18V,- t - 9V,-2 + 2V,-3)
d = £(18p.
- 9p,-2 + 2p,-3) C
e = Y,-I + ~ - ~ ( 1 8 T . - 9T~-2 + 2T~-3) K
+ 6-'~ (18P,-I - 9P,-2 + 2P,-3). If the energy equation is given in the form of (13), then K = 0 and S are replaced by L = L.~.
(21)
T. J.
CHun~ and
J. N. Cmot:
Initial and boundary conditions may be applied using the Lagrange multipliers or simply rearranging the matrix arrays of the global form of (21) upon assembly. The assembled form of (21) may be written as fL..X. = A.,
t22)
where m , n = 1,2 . . . . . r with r = total number of equations, and X. represents the solution vector consisting of all unknowns of velocity, density, temperature, and pressure. In general, homogeneous boundary conditions may be written in the form, bq,.X,,, = 0
(23)
where bqm is the boundary condition matrix with q = !, 2 . . . . . s and s = total number of boundary conditions. We now introduce the Lagrange multipliers A~ such that Aqb..X.=
(24)
O.
Consider, in view of (22) and (23), a quadratic functional ~ in the form = ½~...x..x.
- A , . X , . + Xqbq..X...
(25)
Minimizing the functional (25) gives 0---~-=0 and ~ - ~ = 0 which may be combined into a form 11,.. b,° ~ ; ' ] [~X~]--?O~J •
(26)
For example, if a boundary condition is described as X~ = X2, then Xt - X: = 0 and the boundary condition matrix is bq.=[b,,
b,2]=[l
-1].
A short form of (26) may be expressed as 1~,...,~. = A .
(27)
This is the modified form of (21) or (22) with appropriate boundary conditions. Here it should be noted that rearrangement of the matrix is necessary to avoid zeros on the diagonal of lar~,. We observe that the simultaneous equations shown in (21) are highly nonlinear. In the present study, the initial and boundary values specified at the first time increment allow the assembled form of (21) to be rearranged such that all coefficients and right-hand-side vector values are known and solution is obtained at the end of the first time increment for all unknowns. Specifically, the Newton-Raphson method may be used to solve these nonlinear equations within each time increment. Dropping * in (27) for simplicity, we may write in the form F(X) = f t ( X ) X - A(X)= 0.
(28)
It is possible to expand (28) in Taylor series and we obtain the Newton-Raphson formula of the form X = Xo - Jo-I F(Xo)
(29)
Analysis of unsteady compressible boundary layer flow via finite elements
in which the Jacobian Jo is given by Jo = aF,(Xo)laXj. For m 'h iteration, (29) becomes (30)
X , =Xm_t-J~LtF(X._t) with
J,.-I
=
OF1 OX, aF~ #X,
OFt OX~ OF2 OX2
aFro OX. OF OX,
OF, aF, aXl OX2
aF, aXn
(31)
'n-I.
With an initial guess for Xo together with boundary conditions, we initiate computation in the first cycle of the first time increment. Upon convergence of solution of (30), we move to the next time increment and the same procedure is repeated until enough time steps are advanced. An approximate Jacobian matrix may be obtained by abandoning the dependency of 11 and A in (28) on X temporarily but recalculating J in terms of the updated values of X at each iteration cycle. This will lead to J = II. Further simplification results if one uses a relationship of the form, xm = t'V'(X,,, _,)A(X,,,_,).
(32)
The process in (30) or (32) may be repeated until we obtain
(33)
Ux~ - x ~ _,ll/llx~,-,ll < e
in which e represents an allowable error. APPLICATIONS
The discrete formulation and numerical procedures in the foregoing sections are applicable for any types of flows since no simplifications or approximations are made specifically for boundary layer flows. However, smaller elements are required in the vicinity of the boundary layer whereas very large elements may be used in other areas. In this manner, singularity of the leading edge as would occur in an analytical solution can be adequately accounted for in the finite element numerical solutions. For a specific application, the boundary layer flow behind the moving normal shock as shown in Fig. 1 is solved. The normal shock wave moves at a constant speed U, into a fluid at rest, as depicted in Fig. l(a). The finite element discretization domain constitutes the stationary shock system in which the shock is stationary and the wall moves at a speed U, is the negative x-direction so that the fluid behind the shock wave also moves in the negative x-direction, as shown in Fig. l(b). To compare the results of the finite element method with those of other methods, a review of
Moving
y
Shock i -
Stationary
y
Shock t__
~Us N
--~ Tw ~-" (a) Movin 9 shock system
~X
tI g7
"-
F
(b) Stationary shock system
Fig. I. Boundary layer flow behind a moving shock wave.
"'-.
N.
---X
T. J. ('tit ~t; and J. N. (.'Htot
the literature concerning the same problem appears to be useful. Let us cite the approximate analytical solution by Mirels[2,31. To solve analytically, many simplifying assumptions were required. Mirels assumes that the flow properties of the external flow are independent of x and t or that the effect of a growing boundary layer on the external flow is negligible. This implies that the Iogitudinal pressure gradient is zero. With these assumptions. Mirels obtained the solutions for the boundary layer velocity and temperature in the forms, 134)
u = UJ'(~)
and T,
-
~ M~"rOT)-
M=2r(O)+ 1-~-~ s(~).
(35)
Here, U~ = velocity outside the boundary layer, [ ( ~ ) = dimensionless streamline function, T, = temperature outside the boundary layer, ~, = cp Ic,,, M~ =Mach number outside the boundary layer, T,~ and T, are wall and recovery temperatures, respectively, the functions, r(~) and s (r/) are solutions satisfying the equations for linear momentum, continuity, energy, and state. A total of 6 isoparametric elements[3] are used in the present study (Fig. 2). The y-dimension is estimated by the so-called Rayleigh's problem. The fluid properties, initial conditions, and boundary conditions for the stationary shock system are summarized in Table i. The finite element solution of the velocity distribution in a laminar boundary layer behind a normal shock is shown in Fig. 3. It is seen that steady state velocities are reached at approximately 3 × 10-6 sec. The steady state velocity profiles at distances of 2 ft and 5 ft behind the shock wave are shown in Fig. 4. The results of the present study and those of Mirels appear to be close, although exact agreement cannot be expected because of the basic difference between the rough approximations in the Mirels work and the discrete solution of the complete governing differential equations. It should be noted that these plots represent velocities in the direction of the shock relative to the wall. The temperature vs time plots as shown in Fig. 5 indicate a slower approach of temperature to a steady state than in the case of velocity. The temperature distribution which has reached a steady state at approximately 12 x 10-6 sec is presented in Fig. 6. The comparison between the finite element solution and the analytical solution is again favorable. The results presented above are based on iterative cyclic solutions of (30) with J = 11 within each time increment. The integration time increment in (20) was 1 ~sec. A proper assessment of the accuracy of the present method is involved in rigorous studies on convergence for solution of nonlinear equations and stability of temporal operators. These tasks are under progress and the preliminary results are reported elsewhere[6].
Shock T i . . . . . . . . T~
®
05
® "2
.025 ~_
® u ~2
......... 5
4;---. 2
4
.01
~,
Q 0
~X
Distance from the leading edge fit) Fig. 2. Finite element model for the boundary layer behind a normal shock wave.
Analysis of unsteady compressible boundary layer flow via finite elements
II
Table 1. Input data L
Fluid Properties u.
external
T o = 200°R
flow temperature
external flow pressure
P, = I0 Ibf/ft 2 p.
external
~s - 1.424 x I0 -q Ibf/see/ft =
fluid viscosity
external flow Reynold's number
= 1.536 x l0 "~ ibf/sec/ft 2
R,
= 4.125 x 105 @ L = 5 ft L Pr = 0.72
Prandtl No. Wall temperature
T~ = 180°R
Constant pressure specific heat
cp = 13.73 Btu/Slug-°R
Constant volume specific heat
e, ~ 5.46 Btu/Slug-°R
Free stream Mach No.
H=
Constants
a = .712 x I0 -~ Ibf/ft s OR, b = -(2/3)a
Initial
Conditions
u = -695 f t / s e c u =-347.5
IIL
@y = 0
ft/sec
Boundary C o n d i t i o n s
elsewhere
~ y = 0
u =-347.5 ft/sec v-0@y=0
T = 180°R @ y = 0
T -
200°R e l s e w h e r e
p = I0 psf
- 2.65
u =-695 f t / s e c
V = 0 everywhere
T "
= 3.38 x 10 -6 slug/ft z
external fluid density
external fluid thermal conductivity ~
II.
" 347.5 f t / s e c
external flow velocity
180°R @ y -
@ y = 0.05 ft
0
T = 200°R @ y - 0 . 0 5 f t
everywhere
P -
0 " 3 . 7 6 x 10 -s S l u g / f t s @
y = 0
10 psf
and x - 0
@ y -
0.05
and x = 0
ft
and x -
0
0 = 3 . 3 8 x I0 "~ s l u g / f t s @ y - 0.05 f t and x = 0
p " 3.38 X 10 -s s l u g / f t 3 e l s e w h e r e
350
3~
310
Node 6
== 290 u
o
"~ >
27C
~Node 5 25C Time
FiB. 3. - -
o
Mirel's
Finite
(tx IOSsec)
Velocityvs time at nodes 5 and 6. Shock
Solution
Element
D4
Solution
boundary I
~
.o3
.°.I OI
~ Distance
~ from
2 the
leading
I
o°
edge ( f t )
Fig. 4. V e l o c i t y distribution in the l a m i n a r y b o u n d a r y layer behind a normal shock wave,
12
1. J. (HI~4, and J. N. ('HIO~
ooii
.'-Node
f95i
.............
6
zf.. /
tgo I.
/
P 0
ZNode
5
(3.
E i
--~
180
4
6
8
~0
r'2
---
Time (txlO6sec) Fig. 5. Temperature vs time at nodes 5 and 6. ---c>
Mirel's Finite
Solutron Element
Solution
Shock
.05
L. . . .
;i .02 0c
_ .-
,of ~ u
--s
i
.. "
4
J
0o
>
Distance from t h e l e a d i n g edge(ft) Fig. 6. Temperature distribution in the laminar boundary layer behind a normal shock wave. CONCI.USIONS
example, the power of the finite element method has been amply demonstrated. The formulation retained all viscous and dissipation terms in the momentum and energy equations. Although the usual simplified boundary layer equations may be used in the finite element method, the present study has demonstrated that the computer program, once coded for general purposes, can handle all possible types of flows and geometries. No special difficulties were encountered in the iterative solution of the nonlinear equations. Reasonably small sizes of elements along the boundary layer and computationally feasible small time increments are the only requirements for an acceptable accuracy. In t h e f o r e g o i n g
Acknowledgements--The support of this study by the U.S. Air Force Office of Scientific Research under Contract No. F4462-69-C-0124 is gratefully acknowledged. Discussions with Drs. S. T. Wu and J. J. Brainerd are also appreciated. REFERENCES I. Proc. 1st Intern. Syrup. Finite Element Applications to Flow Problems (Edited by J. T. Oden, O. C. Zienkiewicz, R. H. Gallagher and C. Taylor), January 1974, The University of Swansea, England. 2. J. T. Oden and C. Wellford, The analysis of flow of viscous fluids by the finite element method. A.I.A.A.J. 1~ 12) (1972). 3. O. C. Zienkiewicz, The Finite Element Method in Engineering Sciem'e. McGraw-Hill, New York (1972). 4. H. Schlichting, Boundary Layer Theory. McGraw-Hill, New York (1968). 5. H. Mirels, Boundary layer behind shock or thin expansion wave moving into stationary fluid. NACA TN 3712 (1956). 6. T. J. Chung, Convergence and stability of nonlinear finite element equations. A.I.A.A J. 1~7) (19751.
ANALYSIS OF UNSTEADY COMPRESSIBLE BOUNDARY LAYER FLOW VIA FINITE ELEMENTS T. J. CHUNG'I"and J. N. Cmou:l: Department of Mechanical Engineering, The University of Alabama in Huntsville, Huntsville, AL 35807, U.S.A.
(Received 15 June 1975) Abstract--This paper is concerned with the discrete formulation and numerical solution of unsteady compressible boundary layer flows using the Galerkin-finite element method. Linear interpolation functions for the velocity, density, temperature and pressure are used in the momentum equation and equations of continuity, energy and state. The coupled nonlinear finite element equations are approximated by a third order Taylor series expansion as temporal operator to integrate in time with Newton-Raphson type iterations performed until convergence within each time step. As an example, a boundary layer problem of a perfect gas behind a normal shock wave is solved. A comparison of the results with those by other method indicates a favorable agreement.
NOMENCLATURE V, P T a~ tu,
F, t, n, E q, h d,j c~ c,, t
p aN, ,SN,3',~,I/',,, A, ~t K
velocity pressure temperature t~, = acceleration Cauchy stress tensor body force per unit mass surface pressure unit vector normal to surface internal energy density per unit mass heat flux per unit surface area heat supply per unit mass deformation rate tensor constant volume specific heat constant pressure specific heat time mass density interpolation functions for velocity, density, pressure, and temperature, respectively viscosity constants thermal conductivity
Indices--upper case and lower case represent element node numbers and spatial dimension numbers, respectively, unless otherwise noted. INTRODUCTION
Supersonic aerodynamics in space technology prompted an interest in compressible, nonsteady boundary layer flows. Boundary layers of this kind may arise behind shock waves and expansion trains in shock tubes, or in a fast-moving body subjected to deceleration or acceleration in its flight and to aerodynamic heating. In this connection, the analytical or numerical solution of unsteady compressible boundary layer flows still pose a formidable task, particularly for complicated geometries and boundary conditions. In recent years, considerable efforts have been made toward application of the finite element method to solve fluid mechanics problems, as evidenced by the recent international symposium on this subject[l]. In the present study, the nonlinear Navier-Stokes differential equations governing the motion are cast in a Galerkin integral approximation, leading to a coupled set of nonlinear integrodifferential finite element equations. The results are similar to the approximation for the NavierStokes equations as obtained by a direct method or the so-called virtual velocity formulation[2]. ?Professor of Mechanical Engineering. tGraduate Student. CAF Vo~ 4. No $--A
1
T.
J. ('l.tlN(i and J. N. CHIOI'
The energy equation written in terms of internal energy and equations of continuity and state are then combined together with the equation of motion for simultaneous solution, integrating across isoparametric finite elements[3] in spatial domain and using suitable temporal operators in time domain. A stable convergent solution is obtained by the temporal operator with Taylor series terms of the third degree in time and Newton-Raphson type iterations in spatial domain within each integration time step. To demonstrate the effectiveness of the present approach, a boundary layer problem of a perfect gas behind a shock wave is solved and compared with the results of others 14, 51.
INTERPOLATION FUNCTIONS AND WEIGHTING FUNCTIONS To derive the finite element analogues for equations of momentum, continuity, energy and state, it is necessary to introduce interpolation functions for the velocity V, density p, pressure P, and temperature T, in the form V,(X, t) = a,~(X)VY(t) p(X, t) =/3N(X)pN(t)
Ilb)
P(X, t) = yN(X)PN(t)
(lc)
T(X, t) = tON(X)TN(t)
(ld)
(la)
where V,N, pN, p~, and T N are the nodal values of velocity, density, pressure, and temperature, which are functions of time only. Note that the interpolation functions aN,/3~,, yN and qJN are dependent on the spatial coordinates X but free from dependence on time. Here i = 1, 2 for two-dimensional flow, N denotes the element node, and repeated indices imply summing. In the Galerkin's integral approximations, these interpolation functions are taken as weighting functions which require the basic convergence criteria of completeness and continuity[2, 3]. In the method of weighted residuals, we construct an inner product of the residual and the weighting functions. Often, the residual consists of more than one variable with each variable assumed to vary spatially with polynomials of different degrees. For example, in the momentum equation of a compressible viscous flow, the variables are the density, pressure and temperature, as well as the velocity. A question arises as to an appropriate weighting function for the momentum equation. The general rule is to seek one of the variables in the differential equation to construct an inner product of the residual of the differential equation and this variable so that the inner product will yield a spatial invariant. In this case, the obvious variable is the velocity, and we obtain an invariant I, = ( ~ , Vk) = (~k, aN(X)V~N(t)). Therefore, the weighting function is the velocity interpolation function aN(X). Note that the density p, pressure P, and temperature T are disqualified because they are scalars and do not possess a free index k required to yield an invariant as a result of an inner product with the residual ~k of the momentum equation. Likewise, for the continuity equation, we form an invariant as an inner product of the residual and the density. Thus I2 = (~, p) = (~, flN(X)pN(t)). Obviously the density interpolation function/3~(X) is the appropriate weighting function. For incompressible fluids, however, the density is constant and v,., = 4 .
Here the only variable involved is the velocity. However, we cannot choose the velocity as a subspace since the inner product (~, V~) does not yield an invariant. In this case, incompressibil-
Analysisof unsteadycompressibleboundarylayerflowviafiniteelements
3
ity requires that Is = (~, P) = (~, Tr~(X)PN). Consequently, the pressure interpolation function 3,N(X) is the chosen weighting function. For the energy equation, we find that temperature, pressure, density, and velocity are all involved as variables. The residual is denoted ~ and all variables except the velocity appear to provide an invariant from the inner product. In the absence of dissipative energy due to velocity gradients, the spatial distribution of energy is contributed by the gradient of heat flux or V~T. Thus it is apparent that an invariant is constructed as I3 = (~, T) = (~, #~(X)T N) and the temperature interpolation 0N(X) is the weighting function. Lastly, the equation of state for a perfect gas assumes an invariant from the inner product with the pressure. Therefore, the pressure interpolation function 7N(X) is taken as the weighting function. In general, it may be argued that, for a very large velocity gradient, higher order approximations for the velocity should be used. In the present study, however, a linear spatial variation of all the variables is assumed for simplicity. In particular, two-dimensional four node linear isoparametric elements are used for the velocity, density, pressure, and temperature. Consequently, the weighting functions for all four equations of momentum, continuity, energy and state are of the polynomial of the first degree. MOMENTUM EQUATION A linear momentum equation is given by pa~
-
pF~ - tl~., = 0
(2)
where ah = (DVk/Dt) = (:~ + V, Vk., is the acceleration, F~ is body force per unit mass, t,k is the Cauchy stress tensor, and the comma denotes a partial derivative. Let ~ be the residual of (2). We then have
3k = pa~ - pFk
-
t,k., = p(f'k + V, Vk.,)- pFk
-
t,~.,.
Using Galerkin's approximation with a~, as weighting functions, we require
f ~lkaN dv = O. This leads to f [p(~Yk+ VrVk.,)-pF~ - ta.,las dv = O.
(3)
For Stokesian fluid, the stress tensor assumes the form tt~ = - P & k
+ Ad,,&~
+ 21~d,k
where, for Stokes hypothesis, the viscosity constants/z and A are related by X = -2/3~., and the deformation rate tensor da is d,k =~(I V ~,k + Vk.,)
Noting that
f t,,.,a,,, dv = f t,~,,,,n, dA - f t,,.a,,,., dv = f t~,,,, dA - f t.,a,,,., dv
T.J.
4
CHt?NG a n d J. N . CHIOU
in which dA is the boundary surface, tk = t,~nk represents the components of boundary stresses and 8,~ is the interpolation function for the boundary surface pressure, N being only the boundary nodes involved instead of all nixies of an element. Assuming the viscosity constants/x and A to be a function of temperature such that #, = aT and A = bT where a and b are constants, we can rewrite (3) in the form, A~,N(/k M + BbN~ V, ~' + Z~tNkP ea = XN~
l~l, 1) ~ XNk
~
(4)
I~(2) ~_ L"(3)
~ N k "T- r" N k T
& " Nk
and A~,N = f fll.aMaN dvp L B ~Nk = C ~Nk + D'MNk C',aNk = f bez.OtM.,CtN.k dv T t" D ~ k = ~ aqSLA'Maa~.~ dv T ~ ZMNk = -- f 7~,aN.k d v
(1)-- f ~.F~a~ dvp R FNk-(2) F Nk -
t~tN dA
1
FN - --f fl~a~c~Lao., dvpMV, LVk ° (3)k -with A ~,~k= a~.k&" + aA,.~&" Here, AMN is the mass matrix, B ~,M~is the viscosity matrix, Z~,N~ is the dilation matrix, F "~ is the body force vector, F ~2~is the surface pressure vector, and F~k is the convective term hereafter called the inertia pseudo force vector.
EQUATION OF CONTINUITY
The equation of continuity is written in the form, + (p Vk ).~ =0. Again, using Galerkin's approximation residual + (pV~).k =
with weighting functions, tiM, we obtain
f
~#M
dv = 0
or
f[,o
"~-(p Vk ),k ]flM
dv =0.
Substituting the interpolation functions and performing the differentiations, we obtain GMNp N + N~Np N = 0
(5)
Analysisof unsteady compressibleboundary layer flowvia finite elements
5
where f GMN= J ~M~N dv HM~ = f ~SM(~NaL.~+ ~8~.~aL)dv VkL. ENERGY EQUATION The first law of thermodynamics requires the balance of energy as well as the conservation of linear momentum. The energy equation is given by De
p - ~ - t,d, - q,., - oh = 0
(6)
where • is internal energy per unit mass, q~ is the heat flux per unit surface area, h is the heat supply per unit mass. The residual function of the energy equation is D~
= p - ~ - tljd~j - qu - ph = p~ + pV~.l - t,do - ql,~ - ph.
Once again, using the Galerkin's approximation with weighting functions ~ , , we have f ~o'~ dv
0
or
f(p~ +
-
- q,., - ph )~bN dv = 0.
(7)
Here, for a perfect gas, = c.T - Pip
(8)
where cp is the specific heat at constant pressure. Making use of the Fourier law, qj = KT.~ where ~ is the thermal conductivity, and in view of the continuity equation, we may rewrite (7) in the form f (pcpi" + p % V~T., - P - A Vk.~ - 21~d, do - V~P., - q,., - ph )~bN dv = O.
(9)
From the Green-Gauss theorem it is possible to write
0o) where ~ = q,n~ is the input heat flux normal to the surface. Substituting (1) and (10) into (9) and simplifying yield CMNi "M + DMNT M + KMNP M + SMNP M = YN
with
and CMn = f Cp[3LOM~bn d v p L DMn = f KOM.I~bN.I dv
(11)
T. J. CHUNG and J. N. CHIOt"
KuN = - f
yuON dv
S~N = - f yM.iaLtPN dv
Vil"
QN"'= .//3,~h~Ndrp M Q ,2,= f (1~. dA QNTM = f b~uOtL,kfltQdd.INdv TMvkLvi Q =
2a~bu~NAL~Ap.dvVk Vt T
Q~'" = f - C,JSLapd/U.,OZ~dvp L V/.T u. For a perfect gas the internal energy is defined alternatively as • = c~T
(12)
where c~ is the specific heat at constant volume. Using (12) instead of (8), one obtains the finite element energy equation of the form (11) without the term containing the pressure rate. In this case, eqn (ll) is modified to CMt~7"M + DMNT u + LuNP u = YN
(13)
where
CuN = f C.[3L~bM$Ndvp L Lu~ = f yuaL.i~b~ dv V,L QN'" = f - CdJLaotPu.,ON dvp L VIOTu with other quantities remaining the same as before. E Q U A T I O N OF S T A T E
The equation of state which will now complete the description of compressible thermal flow is given by
P = pRT
(14)
where R is the gas constant. Proceeding as before, but this time with weighting functions y~,, we obtain
f (P - pRT)yN dv = 0 or
Eu~P u = fN where EuN = f y u y ~ dr;
fN = R f flL~Uy~ d v p L T u.
(15)
Analysisof unsteadycompressibleboundarylayerflowviafiniteelements
7
SOLUTION PROCEDURE The local finite element equations for momentum, continuity, energy, and state assume the global forms upon assembly[2, 3], respectively, A~,~'~ ~ + B ' ~ V 7
+ Z.,~,P ° = X~,
(16) (17)
God~" +H°op ° = 0 CoaT ° + D o l T ~" + K,,¢~ ° + Jo~P*" = Yo E,,j~P*=fa
(18)
(19)
with a,/3 = 1, 2. . . . . & ~ being the total number of nodes in the entire domain under investigation. Note that integration of all element coefficient matrices must be carried out before assembly. A computerized Gaussian quadrature may be easily applied for the integration. The time derivatives may be approximated by various difference operators. In the present study, the 3rd order Taylor series expansion for the temporal operator was used. The first time derivatives for velocity, density, temperature, and pressure may be expanded in terms of histories of the past time steps. For example, if g represents any one of the unknown vectors, then we have[6] 1
lg, = ~-~-(1 lg, - 18g,_, + 9g,-2 - 2g,_,)
(20)
where At is the time increment at any time step i. Thus, the incremental equations in recursive matrix form for momentum, continuity, and energy become 1 IA + B'~V ] i+ZP,=X.
IlG+
x
+ OAt .A (18V,_1 - 9V,_2 + 2V,-,)
G
H)p, =-~(18p.-9p,_:
+ 2p.)
(~-~t + D)T, + (~-~t + S)P, = Y,-I + £(18"1',-i - 9T,-2 + 2T,_,) + ~At ( 1 8 P . - 9Pi-2 + 2P~-s). To these the equation of state must be added. Writing all governing equations in a matrix form, we obtain N 0 0
0 R 0
0p= ST
Ii ° ° !llilIll
where IIA M=6--~'+B, c = X.
I1G_ N=~--~-+H,
_ 11C IlK R=6--~-+D , S = 6 - ~ - + J
A
+ ~-~'(18V,- t - 9V,-2 + 2V,-3)
d = £(18p.
- 9p,-2 + 2p,-3) C
e = Y,-I + ~ - ~ ( 1 8 T . - 9T~-2 + 2T~-3) K
+ 6-'~ (18P,-I - 9P,-2 + 2P,-3). If the energy equation is given in the form of (13), then K = 0 and S are replaced by L = L.~.
(21)
T. J.
CHun~ and
J. N. Cmot:
Initial and boundary conditions may be applied using the Lagrange multipliers or simply rearranging the matrix arrays of the global form of (21) upon assembly. The assembled form of (21) may be written as fL..X. = A.,
t22)
where m , n = 1,2 . . . . . r with r = total number of equations, and X. represents the solution vector consisting of all unknowns of velocity, density, temperature, and pressure. In general, homogeneous boundary conditions may be written in the form, bq,.X,,, = 0
(23)
where bqm is the boundary condition matrix with q = !, 2 . . . . . s and s = total number of boundary conditions. We now introduce the Lagrange multipliers A~ such that Aqb..X.=
(24)
O.
Consider, in view of (22) and (23), a quadratic functional ~ in the form = ½~...x..x.
- A , . X , . + Xqbq..X...
(25)
Minimizing the functional (25) gives 0---~-=0 and ~ - ~ = 0 which may be combined into a form 11,.. b,° ~ ; ' ] [~X~]--?O~J •
(26)
For example, if a boundary condition is described as X~ = X2, then Xt - X: = 0 and the boundary condition matrix is bq.=[b,,
b,2]=[l
-1].
A short form of (26) may be expressed as 1~,...,~. = A .
(27)
This is the modified form of (21) or (22) with appropriate boundary conditions. Here it should be noted that rearrangement of the matrix is necessary to avoid zeros on the diagonal of lar~,. We observe that the simultaneous equations shown in (21) are highly nonlinear. In the present study, the initial and boundary values specified at the first time increment allow the assembled form of (21) to be rearranged such that all coefficients and right-hand-side vector values are known and solution is obtained at the end of the first time increment for all unknowns. Specifically, the Newton-Raphson method may be used to solve these nonlinear equations within each time increment. Dropping * in (27) for simplicity, we may write in the form F(X) = f t ( X ) X - A(X)= 0.
(28)
It is possible to expand (28) in Taylor series and we obtain the Newton-Raphson formula of the form X = Xo - Jo-I F(Xo)
(29)
Analysis of unsteady compressible boundary layer flow via finite elements
in which the Jacobian Jo is given by Jo = aF,(Xo)laXj. For m 'h iteration, (29) becomes (30)
X , =Xm_t-J~LtF(X._t) with
J,.-I
=
OF1 OX, aF~ #X,
OFt OX~ OF2 OX2
aFro OX. OF OX,
OF, aF, aXl OX2
aF, aXn
(31)
'n-I.
With an initial guess for Xo together with boundary conditions, we initiate computation in the first cycle of the first time increment. Upon convergence of solution of (30), we move to the next time increment and the same procedure is repeated until enough time steps are advanced. An approximate Jacobian matrix may be obtained by abandoning the dependency of 11 and A in (28) on X temporarily but recalculating J in terms of the updated values of X at each iteration cycle. This will lead to J = II. Further simplification results if one uses a relationship of the form, xm = t'V'(X,,, _,)A(X,,,_,).
(32)
The process in (30) or (32) may be repeated until we obtain
(33)
Ux~ - x ~ _,ll/llx~,-,ll < e
in which e represents an allowable error. APPLICATIONS
The discrete formulation and numerical procedures in the foregoing sections are applicable for any types of flows since no simplifications or approximations are made specifically for boundary layer flows. However, smaller elements are required in the vicinity of the boundary layer whereas very large elements may be used in other areas. In this manner, singularity of the leading edge as would occur in an analytical solution can be adequately accounted for in the finite element numerical solutions. For a specific application, the boundary layer flow behind the moving normal shock as shown in Fig. 1 is solved. The normal shock wave moves at a constant speed U, into a fluid at rest, as depicted in Fig. l(a). The finite element discretization domain constitutes the stationary shock system in which the shock is stationary and the wall moves at a speed U, is the negative x-direction so that the fluid behind the shock wave also moves in the negative x-direction, as shown in Fig. l(b). To compare the results of the finite element method with those of other methods, a review of
Moving
y
Shock i -
Stationary
y
Shock t__
~Us N
--~ Tw ~-" (a) Movin 9 shock system
~X
tI g7
"-
F
(b) Stationary shock system
Fig. I. Boundary layer flow behind a moving shock wave.
"'-.
N.
---X
T. J. ('tit ~t; and J. N. (.'Htot
the literature concerning the same problem appears to be useful. Let us cite the approximate analytical solution by Mirels[2,31. To solve analytically, many simplifying assumptions were required. Mirels assumes that the flow properties of the external flow are independent of x and t or that the effect of a growing boundary layer on the external flow is negligible. This implies that the Iogitudinal pressure gradient is zero. With these assumptions. Mirels obtained the solutions for the boundary layer velocity and temperature in the forms, 134)
u = UJ'(~)
and T,
-
~ M~"rOT)-
M=2r(O)+ 1-~-~ s(~).
(35)
Here, U~ = velocity outside the boundary layer, [ ( ~ ) = dimensionless streamline function, T, = temperature outside the boundary layer, ~, = cp Ic,,, M~ =Mach number outside the boundary layer, T,~ and T, are wall and recovery temperatures, respectively, the functions, r(~) and s (r/) are solutions satisfying the equations for linear momentum, continuity, energy, and state. A total of 6 isoparametric elements[3] are used in the present study (Fig. 2). The y-dimension is estimated by the so-called Rayleigh's problem. The fluid properties, initial conditions, and boundary conditions for the stationary shock system are summarized in Table i. The finite element solution of the velocity distribution in a laminar boundary layer behind a normal shock is shown in Fig. 3. It is seen that steady state velocities are reached at approximately 3 × 10-6 sec. The steady state velocity profiles at distances of 2 ft and 5 ft behind the shock wave are shown in Fig. 4. The results of the present study and those of Mirels appear to be close, although exact agreement cannot be expected because of the basic difference between the rough approximations in the Mirels work and the discrete solution of the complete governing differential equations. It should be noted that these plots represent velocities in the direction of the shock relative to the wall. The temperature vs time plots as shown in Fig. 5 indicate a slower approach of temperature to a steady state than in the case of velocity. The temperature distribution which has reached a steady state at approximately 12 x 10-6 sec is presented in Fig. 6. The comparison between the finite element solution and the analytical solution is again favorable. The results presented above are based on iterative cyclic solutions of (30) with J = 11 within each time increment. The integration time increment in (20) was 1 ~sec. A proper assessment of the accuracy of the present method is involved in rigorous studies on convergence for solution of nonlinear equations and stability of temporal operators. These tasks are under progress and the preliminary results are reported elsewhere[6].
Shock T i . . . . . . . . T~
®
05
® "2
.025 ~_
® u ~2
......... 5
4;---. 2
4
.01
~,
Q 0
~X
Distance from the leading edge fit) Fig. 2. Finite element model for the boundary layer behind a normal shock wave.
Analysis of unsteady compressible boundary layer flow via finite elements
II
Table 1. Input data L
Fluid Properties u.
external
T o = 200°R
flow temperature
external flow pressure
P, = I0 Ibf/ft 2 p.
external
~s - 1.424 x I0 -q Ibf/see/ft =
fluid viscosity
external flow Reynold's number
= 1.536 x l0 "~ ibf/sec/ft 2
R,
= 4.125 x 105 @ L = 5 ft L Pr = 0.72
Prandtl No. Wall temperature
T~ = 180°R
Constant pressure specific heat
cp = 13.73 Btu/Slug-°R
Constant volume specific heat
e, ~ 5.46 Btu/Slug-°R
Free stream Mach No.
H=
Constants
a = .712 x I0 -~ Ibf/ft s OR, b = -(2/3)a
Initial
Conditions
u = -695 f t / s e c u =-347.5
IIL
@y = 0
ft/sec
Boundary C o n d i t i o n s
elsewhere
~ y = 0
u =-347.5 ft/sec v-0@y=0
T = 180°R @ y = 0
T -
200°R e l s e w h e r e
p = I0 psf
- 2.65
u =-695 f t / s e c
V = 0 everywhere
T "
= 3.38 x 10 -6 slug/ft z
external fluid density
external fluid thermal conductivity ~
II.
" 347.5 f t / s e c
external flow velocity
180°R @ y -
@ y = 0.05 ft
0
T = 200°R @ y - 0 . 0 5 f t
everywhere
P -
0 " 3 . 7 6 x 10 -s S l u g / f t s @
y = 0
10 psf
and x - 0
@ y -
0.05
and x = 0
ft
and x -
0
0 = 3 . 3 8 x I0 "~ s l u g / f t s @ y - 0.05 f t and x = 0
p " 3.38 X 10 -s s l u g / f t 3 e l s e w h e r e
350
3~
310
Node 6
== 290 u
o
"~ >
27C
~Node 5 25C Time
FiB. 3. - -
o
Mirel's
Finite
(tx IOSsec)
Velocityvs time at nodes 5 and 6. Shock
Solution
Element
D4
Solution
boundary I
~
.o3
.°.I OI
~ Distance
~ from
2 the
leading
I
o°
edge ( f t )
Fig. 4. V e l o c i t y distribution in the l a m i n a r y b o u n d a r y layer behind a normal shock wave,
12
1. J. (HI~4, and J. N. ('HIO~
ooii
.'-Node
f95i
.............
6
zf.. /
tgo I.
/
P 0
ZNode
5
(3.
E i
--~
180
4
6
8
~0
r'2
---
Time (txlO6sec) Fig. 5. Temperature vs time at nodes 5 and 6. ---c>
Mirel's Finite
Solutron Element
Solution
Shock
.05
L. . . .
;i .02 0c
_ .-
,of ~ u
--s
i
.. "
4
J
0o
>
Distance from t h e l e a d i n g edge(ft) Fig. 6. Temperature distribution in the laminar boundary layer behind a normal shock wave. CONCI.USIONS
example, the power of the finite element method has been amply demonstrated. The formulation retained all viscous and dissipation terms in the momentum and energy equations. Although the usual simplified boundary layer equations may be used in the finite element method, the present study has demonstrated that the computer program, once coded for general purposes, can handle all possible types of flows and geometries. No special difficulties were encountered in the iterative solution of the nonlinear equations. Reasonably small sizes of elements along the boundary layer and computationally feasible small time increments are the only requirements for an acceptable accuracy. In t h e f o r e g o i n g
Acknowledgements--The support of this study by the U.S. Air Force Office of Scientific Research under Contract No. F4462-69-C-0124 is gratefully acknowledged. Discussions with Drs. S. T. Wu and J. J. Brainerd are also appreciated. REFERENCES I. Proc. 1st Intern. Syrup. Finite Element Applications to Flow Problems (Edited by J. T. Oden, O. C. Zienkiewicz, R. H. Gallagher and C. Taylor), January 1974, The University of Swansea, England. 2. J. T. Oden and C. Wellford, The analysis of flow of viscous fluids by the finite element method. A.I.A.A.J. 1~ 12) (1972). 3. O. C. Zienkiewicz, The Finite Element Method in Engineering Sciem'e. McGraw-Hill, New York (1972). 4. H. Schlichting, Boundary Layer Theory. McGraw-Hill, New York (1968). 5. H. Mirels, Boundary layer behind shock or thin expansion wave moving into stationary fluid. NACA TN 3712 (1956). 6. T. J. Chung, Convergence and stability of nonlinear finite element equations. A.I.A.A J. 1~7) (19751.