# Numerical analysis of unsteady compressible laminar boundary layer flow

## Numerical analysis of unsteady compressible laminar boundary layer flow

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 0 NORTH-HOLLAND PUBLISHING COMPANY 19 (1979) 10-204 NUMERICAL ANALYSIS OF UNSTEADY ~O~PR~SIBLE...

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 0 NORTH-HOLLAND PUBLISHING COMPANY

19 (1979) 10-204

NUMERICAL ANALYSIS OF UNSTEADY ~O~PR~SIBLE LAMINAR BOUNDARY LAYER FLOW NGUYEN J.R.C. ISPRA,

System Analysis

Hung

Diuision,

i-21020 ISPRA

(VA), Italy

Received 3 October 1977 The unsteady compressible laminar boundary layer flow on an arbitrary cylinder due to an incident stream whose velocity varies arbitrarily with time is considered. The method presented is based on the separation between the convective and diffusive quantities. By defining some new variables, the splitting appears rather naturally, and the initialisation problem can be solved without difficulty. The transformed equations are solved with the help of a semi-implicit finite difference scheme which is unconditionally linearly stable. The computations have been applied to flows past a cylinder with constant and fluctuating free-stream velocities.

1. Introduction The subject of unsteady boundary layers - that is the addition of unsteadiness to viscous flows near boundaries has again attracted considerable attention in recent years. The proceedings of the IUTAM symposium on unsteady boundary layers [l] trace out an impressive and vivid picture of the various activities in this particular field. Unsteady boundary layers play an important role for the understanding of particular flow problems in many technical applications such as turbomachines, helicopter rotors, shock tubes, propellers, unsteady flight problems, duct and nozzle flow, combustion chambers, charging and discharging motions etc. In order to determine the friction drag and the rate of heat transfer through the surface for such unsteady motions, the nature of time-dependent laminar compressible boundary layer flow of a fluid past an obstacle must be investigated. Attempts to obtain practical solutions, exact or approximate, to the complete set of boundary layer equations lead to very great difficulties. Three basic approaches have been successfully developed. The first one, usually refered to as the integral method, is a natural extension of the well-known von Karman-Pohlhausen procedure to unsteady boundary layers. Mathematically speaking, these methods belong to a wider class of approximate methods based upon weighted residuals. Using this method, Lighthill [2] has considered the case of an arbitrary cylinder in a stream which fluctuates in magnitude but not in direction. Using the same procedure, Yang [3], [4] has studied the unsteady two-dimensional incompressible laminar boundary layer on an arbitrary cylinder subjected to a stream whose velocity varies arbitrarily with time. Unfortunately, the accuracy of the results obtained from these methods depends strongly on the choice of the single parameter family of profiles for velocity and temperature. Moreover, integral methods are only approximate, not exact. The second way uses series expansion methods (Gribben [s], Zeytounian [6]), but they are somewhat laborious to apply and require many terms for large times.

188

H. Nguyen, Numerical unalysis of unsteady compressible lominur boundury luyerflorv

The third method is the finite difference numerical technique. For compressible flows, results have been presented by Piquet [7] and Vimala and Nath [8]. Using Crocco’s transformation, Piquet has calculated some velocity and temperature boundary layers on a flat plate, especially behind a shock. However, this transformation has some disadvantages, the most severe of which is the problem of uniqueness of the solutions in cases where the velocity profile exhibits an overshoot. Vimala and Nath have considered the unsteady laminar compressible boundary layer flow in the immediate vicinity of a two-dimensional stagnation point due to an incident stream with a velocity varying arbitrarily with the time. The governing partial differential equations, involving both time and the independent similarity variable, are transformed into new coordinates with finite range by means of a transformation which maps an infinite interval into a finite one. The solutions given by this method are interesting. Unfortunately, problems of complex nature, such as nonsimilar and semisimilar flows, have not been investigated by the authors. The present paper deals with the unsteady laminar compressible boundary layer flow on an arbitrary cylinder, where unsteadiness arises from arbitrary time variations of the velocity of the incident steam. The method bescribed here is based on the separation between the convective and the diffusive quantities. By defining some new variables, we can obtain a system of four simultaneous equations in which this separation is achieved. The computations have been carried out for the flow past an elliptic cylinder with the following free-stream velocity distributions: (i) a stream moving with constant velocity (elliptic cylinder started impulsively), (ii) a stream fluctuating about a steady mean. As for the unsteady boundary layer over an impulsively started cylinder, using our transformation and the Blasius hypothesis, the initialisation problem can be solved without difficulty. The fluctuating flow problem is solved as an initial value problem starting with a steady solution. When the ratio of the semiaxes of the elliptic cylinder becomes infinity, our method gives the solution of the stagnation flow.

2. Equations of motion The equations of motion for the unsteady, two-dimensional, laminar boundary layer flow of a compressible fluid past an obstacle can be written in nondimensional form as

S&r’)+ \$clr”u) +J\$lr”u) =0, aY

P

(

au

au

Sz+“ax+“ay

ah

ah

s,+uax+“aJ:

au ap a =-ax+ay > ah

>

au

(> Pay’

=(P-_I)M:(S~+~~)+a~(,~)+(~-~)M~~(~~,

(2)

(3)

H. Nguyen, Numerical unulysis

of unsteody compressible luminur boundury layerflow

189

where t is the time; x and y are’ curvilinear coordinates along and ~rpe~dicular to the boundary; u and u are the corresponding velocity components; and p, p, p, h, Pr, S and M, are the pressure, density, viscosity, enthalpy, Prandtl number, Stourhal number and Mach number, respectively. The obstacle is defined by r(x); E = 1 corresponds to a body of revolution and 8 = 0 to a flat plate (fig. 1). As in Gribben [5] we assume that (a) the gas is perfect, (b) the viscosity is proportional to the absolute temperature T, and (c) the compressibility effects within the boundary layer arise from an arbitrary (constant) wall temperature while the flow external to the boundary layer is incompressible. These assumptions are expressed by: pp = p,p,

=

T, = const,

const,

Pe = const,

where the suffixes w and e refer to the wall condition and to main-stream Terms of the order of the Mach number M, are also neglected. The boundary conditions are

u=u=o,

h = h, = const

at y = 0,

u = %&, 0,

h = h, = const

at y-‘co.

values, respectively.

(4)

However, in order that the problem be well-posed, it is necessary to give the form of the velocity and the enthalpy at t = t, and at x = .x0. The outer velocity and the enthalpy are connected with the pressure by the relations

he - l)M2,

1

+ ‘il2 2 e *

In order to avoid the explicit calculation of the density in the basic set of equations, the transformation of Dorodnitsyn will be employed :

s Y

5 = x,

vl =

pr”dy,

-C=

t.

Fig. I. Boundary layer near a body of revolution (system of coordinates)

(6)

190

H. Nguyen, Numerical unulysis oJ’unsteudy compressible lantinar boundary luyer,flow

For small Mach numbers - in which case the dissipation term may be neglected - eq. (l)-(3) are now written as

ah

ah

,ah

s,+.,,+,,=pr,

+ a

ah

( > PP,,

3

where

Putting N = ppr2” and changing the notation thus: 4: + x,

‘I -+ Y,

r--t 4

6 --f u,

we obtain

aU+%() ax

ay

7

(11)

(12)

(13) As we have seen, in the present problem with a slow flow, compressibility effects do not extend to the main stream. We may therefore consider the unsteadiness of the stream of the type discussed by Lighthill [2] in his work on incompressible boundary layers, and we write

Consider an elliptic cylinder with semiaxes a and b, and let

H. Nguyen, Numerical analysis

be their ratio. The velocity distribution

191

laminar boundary layerflow

V(x) is then given by (see [9])

(1 + 1)sin 8 V(x) = (sin20 + l2 cos2e)1/2 ’

(14)

with (1 + 1)P cos 8 dV dx = (sin20 + I2 cos28)2 ’

(15)

where

x =

(sin28 + l2 c0s2ep2 de,

with 0 E [0,17].

(16)

s The function 4(t) reflects the nature of the unsteadiness in the external stream; it is defined by 0

1

if

40) = 1 if

t
(17)

t>O

for an elliptic cylinder started impulsively and by 4(t) = 1 + Esin t

(E is

a small constant)

(18)

3. Proposed transformation By defining some new variables, we can split the equations of motion to yield separate equations for the convected and diffused quantities. Let

(19)

H. Nguyen, Numerical analysis of unsteady compressible laminar boundary layerflow

192

Introducing G=;,

and

e

we obtain, instead of (1 I), (12), (13),

Sag = -;

'SF

+ ueq5F2 + P(x,

t)G

+‘\$,

(20)

e

1 a2G Pr ay2’

S;=u,+FG+--

(21) (22)

W=

-(g+c;),

(23)

=;s-Z+\$

(24)

with P(x,t)

e

This system may be written in the matrix form

a32 ay27

(25)

ay> (Lm+g+cm,

(26)

S+2+Bn=-

where the vectors and the matrices are given by

;2 + A=

[

e

0

u&F

P(x, 0

1, B=[; i],

uell/F

H. Nguyen,

Numerical

analysis

compressible

laminar

boundary

layerflow

193

The boundary conditions are y = 0:c = 0,

h with a = 2, he

(27)

y+o:*=[-yq; *=[;].

The above transformation has been applied successfully to the problem of unsteady incompressible laminar boundary layer [lo]. The advantages are obvious: (i) The main equation (25) is linear. (ii) The tirst spatial derivatives are absent in the parabolic equation (25) (thus the knowledge of the phenomenon results only from that of the flow characteristics at the initial time). (iii) The profiles of the reduced velocity and the reduced enthalpy can be obtained by solving eq. (25) and (26) directly.

4. Numerical analysis 4. I. Difference equations

Let x,, = nAx = nk,

n=0,1,2

,... N,

y, = mAy = mh,

m=0,1,2

,... M,

ti = t, +

iAt = t, + iz,

i = 0, 1,2, . . . ,

where Ax = k, Ay = h, At = z, are the space and time increments, as shown schematically in fig. 2. The approximate solution is denoted by

Fig. 2. Finite difference

boundary

layers

194

H. Nguyen, Numerical analysis

n~i

~

s2(X,,

luminar boundary layer Jlo,\

ym, ti).

In addition, we introduce a step parameter At s=(4?,)2>0.

The time derivative is written as the difference m

n,i

C-J =

at

Lp;i+l -

qy) +

O(z),

(28)

and the partial derivative with respect to x is expressed as m

n,i

1

(6)

2k(fi;+l,i

=

i3X

-

a;- lJ) + O(k2).

(29)

The following relations are used to eliminate the partial derivatives with respect to y (see [l 11,

WI)

(f\$)“:l+lO(\$)“i

+ (\$)Y_,=

;([email protected]\$,

- ut\$

+ a::,)

+ O(P)

(30)

and

(g)“:l +4(!?)y+(g)y_,=;(f2::I-n:‘J

+

ow.

(31)

With the above expressions eq. (25) can be expressed by means of two semi-implicit schemes: X?RL X#‘\$;’

1

+

Y,,\$2;* + ZJ2;: + Y#;‘+

r = Rh,

1 + z;i_+l1 = R,,

II= O,l,...N,

(32)

n=O,l,...N,

(33)

where the vectors and coefficient matrices are given by

(34)

H. Nguyen, Numerical analysis of unsteady compressible laminar boundary layerflow

195

and H, and 6; are the operators defined by

q+m1 = AIf1 -

(35)

w, + dL-1.

From the known values of the dependent variables a:’ at i corresponding to time ti = t, + iz, the dependent variables &I:* can be computed by solving the tridiagonal system (32) for 1,. . . A4 - 1 by a very efficient method called the method of factorization, n = o,... Nandm= or the double-sweep method, which is a particular form of the Gauss elimination procedure. Having determined the values F:* and G:* for all values m and n, we may proceed to compfute the values of cz* from eq. (22) written in the difference form F&y 1 _ p* c h* m m+1

1

F”,*_ m

I-

_

_h 1 au, n,i +l UC: 3 ( u, ax )

I +

4F;* + F;:,) (36)

with

The values of &* and \$:* are given by (37) with

and by

H. Nguyen,

196

Numricul

unulysis

ofunstrudy cornprrssible luminur boundury layer /h

Then, R,, is known, and the values of ft\$+l can be obtained from the eq. (33). Replacing the asterisk by i + 1 in the eqs. (36), (37) and (38). we may proceed to compute the values of czi’i, &i+ 1 and +;i+ 1. 4.2. Boundary and initial conditions In terms of the new notation the wall condition at y = 0 is written (39) The boundary condition at y -+ + can then be approximated

by (40)

(41) where E is a small quantity which, for the present calculations, was given the value 10W4. The determination of the initial condition depends on the nature of the unsteadiness in the external stream : (i) Consider firstly the elliptic cylinder started impulsiveIy from rest and suppose that the Blasius hypothesis is valid. Then, for short times, the convective terms of the momentum and energy equations are much smaller than the unsteady and viscous terms [9], and hence the zeroth-order problem of the asymptotic expansion becomes

(42)

&?f= 1 a2h at y = 0:

y--+ m:

(431

Pr ay2’ u = 0, i.l=Li 0

h = h,, h = h,.

(44)

The familiar solutions to the above problems are

(45) h = (h,

H. Nguyen, Numerical unulysis

where er f( ) is the error function values

lurtrimv boundary

luyerflorc

197

The above solutions are evaluated at t = f0 = 0.005, and the

(46)

serve as initial conditions. (ii) Suppose now that the stream fluctuates about a steady mean; the initial conditions are the steady solutions of the previous ease. The boundary conditions at x = 0 may be calculated from eqs, (321,(33), (361, (37) and (38) for Iz = 0. 4.3. 7”Cteskin friction nnd the rate of heat transfer The quantities of physical interest are the skin friction and the rate of heat transfer through the surface. The equation defining the wall skin friction is (47) The rate of heat transfer from the wall to the fluid is

Thus @Flay), and @G/ay),, the velocity and enthaipy gradients at the surface, are the critical skin friction and heat transfer parameters of the flow.

The method of solution described in the previous section has been programmed in PL/I, and solutions have been obtained on an IBM 3705158 computer for two different unsteady free stream velocity distributions characterized by eqs. (1’7)and (18). Computations have been carried out with tl =

0.5,

Ay = 0.15

and

At = 0.01.

198

compressible

laminur

boundary

layer.florc

5.1. An elliptic cylinder started impulsively from rest X1.1. A circular cylinder (1 = 1)

This case is calculated with Prandtl number Pr = 0.72 corresponding

s=

to the case of air and Strouhal number

1.

The computer program requires approximately 650 seconds of execution time. The time required for zero skin-friction to appear at the rear stagnation point is t, =

0.515.

This value is larger than the time obtained in the incompressible case (t, = 0.325). It can be explained by the interaction between the velocity and the thermal boundary layers. Fig. 3 depicts the evolution of the velocity and enthalpy profiles with time at x = n. The evolution of the skin-friction and heat-transfer parameters with time are shown in fig. 4. Our solution indicates that the point of zero skin-friction tends asymptotically to the steady position of separation. However, since the flow presents a point of separation in the steady case, it is impossible to compute beyond a certain time t*. This has been pointed out by several authors for the incompressible case 1131, [14-J, [15]. C onsequently, it is interesting to propose a criterion for

0

I

I

I,

3

6

9

Y

Fig. 3. Velocity and enthalpy profiles at x = x on an impulsively moved circular cylinder

H. Nguyen,

Numerical

analysis

Fig. 4. Evolution

199

of 7,. and qu.

accepting the numerical solution obtained from our numerical procedure and to tell us when to stop the calculations. Our criterion is based on the quantity Q(x, t) which is defined from the wall friction by the relation

Qk t) =

‘2 +uR(x,t) ‘2,

(49)

where uR(x, t) = min 24(x,y, t).

(50)

Y

This quantity Q(x, t) possesses some interesting properties. Indeed, if the flow becomes stationnary, we have

3-0 at

-

and

uR- -0

(thus Q(x, t) = 0).

In the unsteady boundary layers before the time t, is required for zero skin-friction to appear at the rear stagnation point the situation is the following:

& at

2

and

uR- -0

(thus Q(x, t) < 0).

After the time t, there exist three possibilities: a) 0 < x < x*, with r,(x*, t) = 0:

H. Nguyen, Numerical analysis of unsteady compressible laminar boundary layerflow

200

az”
and

b) x** < x < I7,

JS

at

c)

< 0

with

IJR- -0

(thus Q(x, t) < 0);

- 0: X=X** -

%v ->o ax

and

UR < 0

(thus Q(x, t) < 0);

and

UR < 0

(thus Q(z, t)< 0

x*
!s at <03

a7

2ax d 0

or

Q(x, t) 2 0).

It is interesting to observe that Q(x, t) is always negative or zero except perhaps in the domain after the time t,. Therefore, the inequality

can be used as a criterion to accept a numerical solution. Fig. 5 shows the evolution of the quantity Q(x, t) with time. We note that, for t < t* = 0.980 and for all values of x E [0, II], the variation of Q(x, t) is continuous. The discontinuity appears at t = t*, and for x, = 7ZI/lO the value Q becomes positive. The criterion proposed is thus violated. The time t* is the limitation time for our calculation for x E [0, II]. However, if we restrict the value of the coordinate x in the domain [0, xJ, we can continue the calculations until we obtain the steady solutions. The evolution of the velocity and enthalpy profiles with time for: x=--

3l-l 5

are shown in fig. 6. 1

6.31

cY_

t = 0.980 n

.-.-A.

p\ v _ ’\

-10

-

R I

-2 I

1

0

-.

,

1.0.815 ‘\

\

\

1.0.515

\

\

\I \

\

\

\

,’

/

/

A_& \

.// ‘I

‘-’

1 :

t=o.205

,I /I

\

t.o.105 \

/ /I

\$ ‘\

.d

1’

Fig. 5. Evolution

of Q(x, t)

* x

H. Nguyen, Numerical

analysis

laminar boundary layerflow

201

Pr I 0.72 s

Fig. 6. Velocity

and enthalpy

-1

profiles at x = 345 on an inpulsively

moved circular

cylinder

5.1.2. The stagnationflow (1 = CO) For 1+ co,

V(x) = 2x

the problem studied is the stagnation flow. In this case the steady flow has no separation point. Then the criterion

is always satisfied, and we can continue the computation as far as possible. So the steady solutions are obtained for t = 2. In fig. 7 we have plotted, the velocity and the enthalpy profiles for different times. The steady solutions obtained here can be used as initial conditions for the unsteady boundary layer problem where a stream fluctuates about a mean steady value.

Pr = 0.7; s zo.5

-

Fig. 7. Velocity

and enthalpy

proliles

for stagnation

flow with

d(t) =

0 if t
202

H. Nguyen, Nurnericul anulysis of unsteady compressible lantinar boundary layer flow

5.2. A streamfluctuating

This problem is treated for the stagnation flow only (I + co). The initial conditions here are perfectly determined by the previous example. The computations are achieved with Prandtl numbers Pr = 0.72,

Pr = 7.2

and

Pr = 10,

S=l

and

s=5

Strouhal numbers

s = 0.5,

and E = 0.1. The evolution of the velocity and the enthalpy profiles is shown in fig. 8. The variation of the critical wall parameters with time is presented in fig. 9 and 10. From fig. 9 we can see that the skin-friction parameter reacts less to the fluctuations of the free-stream velocity than the heat-transfer parameter. Our results are the same as these of Vimala and Nath [S]. It is clear from fig. 9 that the maximum amplitude of the skin-friction and the heat transfer parameters change only slightly with Strouhal numbers. Fig. 10 shows that the maximum amplitude of the heat-transfer parameter is an increasing function of Prandtl numbers.

6. Conclusions The flow in unsteady laminar compressible boundary layers past cylinders and in the immediate neighbourhood of a two-dimensional stagnation point due to an incident stream with unsteady velocity has been investigated by applying an accurate numerical method for solving nonsimilar boundary layer equations. The equations are solved after having introduced a transformation with the help of a semi-implicit finite difference scheme which is unconditionally linearly stable.

0.5

0(t)

E Pr s

Fig. 8. Velocity

and enthalpy

profiles for stagnation

=I+Esint =O.l ~0.72 =l

flow with

4(t) = 1 + E sin t

H. Nguyen,

Numerical

unnIysis oftinstwdg

flow

compressible Iaminar bowndary Layer

203

When a stream moves with constant velocity and the flow presents a separation point in the steady state (elliptic cylinder), a criterion is proposed to stop the computation. 3t may be noted that the stagnation flow for the same free stream velocity distribution has no separation point; we can then continue the computation as far as possible, and a steady solution is obtained with excellent accuracy.

Acknowledgment The author wishes to thank Dr. H.J. Wirz for valuable remarks on this paper,

0 466

Fig. 10. Variation of[r\$).

- (F),;,1.

for stagnation flow with \$~(t)= t. + csin t; E = 0.1; s = !

H. Nguyen, Numerical analysis

204