Finite-different solutions for compressible laminar boundary-layer flows of dusty gas over a semi-infinite flat plate

FINITE-DIFFERENCE SOLUTIONS

FOR COMPRESSIBLE LAMINAR BOUNDARY-LAYER FLOWS OF A DUSTY GAS

OVER A SEMI-INFINITE FLAT PLATE

August 1986

5

J

11987

by

B. Y. Wang and I. I. Glass

UTIAS Report No. 311

CN ISSN 0082-5255

~ ~

PLATE1: ILLUSTRATION OF A DUSTY-GAS BOUNDARY LAYER. THE PHOTOS SHOW THE

AIR-SAND BOUNDARY-LAYER DEVELOPMENT STARTING FROM THE SHORE OF THE GULF OF AQUABA, JORDAN, AS VIEWED FROM EILAT, ISRAEL (PHOTOS BY

1.

I. GLASS).

FINITE-DIFFERENCE SOLUTIONS

FOR COMPRESSIBLE LAMINAR BOUNDARY-LAYER FLOWS OF A DUSTY GAS

OVER A SEMI-INFINITE FLAT PLATE

by

B. Y. Wang and I. I. Glass

Submitted February, 1986

UTIAS Report No. 311

CN ISSN 0082-5255

Acknowledgements

We wish to thank Dr. W. S. Liu and Prof. J. P. Sislian for their valuable discussions.

One of us (B. Y. Wang) is grateful to the Institute of Mechanics, Academi a Si n i ca, Bei j i ng, Chi na, and to UTIAS for the opportun ity to study and to do research during the past two years.

The financial assistance received from the Natural Science and Engi neeri ng Research Counci 1 of Canada under grant No. A1647, the U.S. Ai r Force under grant AF-AFOSR-82-0096, from the U.S. Defence Nuclear Agency under DNA Contract 001-85-0368, and the Defence Research Establishment Suffield (DRES), is acknowledged with thanks.

...

Summary

A finite-difference method is used to investigate compressible, laminar boundary-layer flows of a dilute dusty gas over a semi-infinite flat plate. Detai 1 s are gi ven of the impl i cit fi nite-difference schemes as well as the boundary conditions, initial conditions and compatibility conditions for solving the gas-particle boundary-layer eqlJations. The flow profiles for both the gas and particle phases were obtained numerically along the whole 1 ength of the pl ate from the 1 eadi ng edge to far downstream of it. The finite-difference solutions in the large-slip region and the small-slip region are compared with the asymptotic solutions and good agreement is achieved. The boundary-layer characteristics of interest, including the wall shear stress, the wall heat-transfer rate and the displacement thickness, are calculated. The alteration of the flow properties owing to the presence of particles is discussed in detail. It was found that the boundary-layer flow of a dusty gas can be ·divided into three distinct flow regi mes whi ch are characteri zed by quas i -frozen, nonequil i bri um and quasi-equilibrium flows and that at a critical distance from the leading edge the particle velocity at the wall decelerates to zero and near-equilibrium is achieved between the gas and particle flows. For the 1 ami nar boundary 1 ayer of a dusty gas, the shear stress and the heat-transfer at the wall are i ncreased and the di spl acement thi ckness is decreased compared with the pure-gas case alone •

'

.

Contents

Acknowledgements

i i

Summary

i i i

Contents

;v

Notation

v

1 • I NTROOUCT

1

ON

1

2. MATHEMATICAL DESCRIPTION OF COMPRESSIBLE, LAMINAR, OUSTY-GAS

BOUNOARY-LAYER FLOWS

5

3. FINITE-OIFFERENCE SCHEMES AND RESULTING FINITE-DIFFERENCE

EQUATIONS

9

4. METHOOS OF SOLUTION OF THE FINITE-OIFFERENCE EQUATIONS

18

5. RELATIONS FOR SHEAR STRESS, HEAT TRANSFER, AND DISPLACEMENT

THICKNESS

23

6. COMPUTER PROGRAM

25

7. NUMERICAL RESULTS AND DISCUSSIONS

27

8. CONCLUD I NG REMARKS

32

REFERENCES

34

FIGURES

APPENDIX A:

OERIVATION OF

THE

FINITE-DIFFERENCE EQUATIONS WITH A

SIX-POINT SCHEME

APPENDIX B:

DERIVATION OF

THE

FINITE-OIFFERENCE EQUATIONS WITH A

FOUR-POINT SCHEME

APPENDIX C:

DERIVATION

OF

THE

RELATIONS

FOR

SHEAR

STRESS,

HEAT

TRANSFER ANO DISPLACEMENT THICKNESS

APPENDIX 0: COMPUTER PROGRAM FDBLEP

APPENDIX E: AN ADDITIONAL DISCUSSION REGARDING THE ASSUMPTION OF THE

PARTICLE-DENSITY PROFILE AFTER THE CRITICAL POINT

D

Notation

coefficients of the finite-difference equations, Eqs.

(3.12)-(3.20) and (3.26)-(3.27)

9

rid parameters for the six-point difference scheme, Eq.

(3.14)

coefficient for fitting a polynomial to the gas temperature

near the wall, Eq. (5.11)

coefficient for fitting a polynomial to the gas velocity

near the wall, Eq. (5.10)

coefficients of the finite-difference equations, Eqs.

(3.12)-(3.20) and (3.26)-(3.27)

9

_{rid parameters for the six-point difference scheme, Eq.}

(3.14)

coefficient for fitting a polynomial to the gas temperature

near the wall, Eq. (5.11)

coefficient for fitting a polynomial to the gas velocity

near the wall, Eq. (5.10)

general drag coefficient for a sphere in a viscous fluid

Stokesian drag coefficient for a sphere in a viscous fluid

coefficients of the finite-difference equations, Eqs.

(3.12)-(3.20) and (3.26)-(3.27)

grid parameters for the six-point difference scheme, Eq.

(3.14)

specific heat at constant pressure for the gas phase

specific heat for the partiele phase

coefficient for fitting a polynomial to the gas temperature

near the wall, Eq. (5.11)

coefficient for fitting a polynomial to the gas velocity

near the wall, Eq. (5.10)

normalized drag coefficient for a sphere in a viscous

fluid

d

Ec

k m N

Nu

n

Pr

p • q R

coefficients of the finite-difference equations, Eqs.

(3.12)-(3.20)

and

(3.26)-(3.27)

diameter of the particles

9

rid parameters for the six-point difference scheme, Eq.

(3.14)

coefficient for fitting a polynomial to the gas temperature

near the wall, Eq. (5.11)

coefficient for fittin9 a polynomial to the gas velocity

near the wall, Eq. (5.10)

gas Eckert number based on the freestream temperature,

Ec

=

u~2/cpT~

recurrence coefficients in the Thomas algorithm, Eq.

(4.1)

integrated function of the nondimensional

displacement

thickness, Eq. (5.14)

recurrence coefficients in the Thomas algorithm, Eq.

(4.1)

transformation function for the asymptotic solution

recurrence coefficients in the Thomas algorithm, Eq. (4.1)

ratio of consecutive step sizes in the y-direction,

K =

I1Yn/I1Yn-1

heat conductivity for the gas phase

grid line in the y-direction

grid point at the outer edge of the boundary layer

Nusselt number based on the particle diameter

grid line in the x-diretion

gas Prandtl number, Pr

=

c*~*/k* _p

gas pressure

heat-transfer rate

gas constant

Re

oo

T

u

v

W iI X Y

Greek Symbol s

Cl

f3

ö e: II

e

À*

00

IJ.

IJ. I

P

Ps

'Ç w

Subscripts

asy

flow Reynolds number based on the particle equilibrium

1

ength, Re

oo

=

p~~À!/

IJ.!

temperature

tangential velocity in the x-direction

normal velocity in the y-direction

function representing any flow property, u, v, P or T

coordinate along the wall

coordinate normal to the wall

ratio of the specific heats for the two phases,

Cl

=

cp/c~

mass loading ratio of the particles in the freestream,

~

=

pp)p!

displacement thickness of the boundary layer

sma

11

quant i ty used i n test i ng for the outer edge of the

boundary layer, Eq.

(4.6)

similarity variable for the asymptotic solution

weighting factor for the finite-difference schemes

particle velocity-equilibrium length, À!

=

p~d*2u~181J.!

dynamic viscosity for the gas phase

deri vat i ve of the gas vi scos i ty with respect to the gas

temperature, IJ.I

=

dlJ./dT

density

density of the particle material

shear stress

power index for the gas viscosity

cri m n p

s

w

o

Superscripts

*

critical point

grid line in the y-direction grid line in the x-direction particle

s 1 i P quant i ty wall conditi ons

freestream conditions initial conditions

di mens i onal

index for dependent variables: i

=

1 , 2 , 3 , 4 , 5 , 6

..

1. INTRODUCTION

Boundary-layer flows of a dusty gas have been investigated using several analytical methods: a series method [1-7], an integral method [8-11], and a fi nite-difference method [12-15]. All the work menti oned above, however, dealt with incompressible-flow cases for the gas phase. Very few authors [16-18] considered the problems of compressible boundary-layer flows where the density of the gas phase can be changed due to compressibility. As pointed out by Singleton [17], Chiu [16] employed incorrect boundary-layer equat i ons and assumed that the parti cl e density i s constant. Si ngl eton extended Marble's analysis [1] to compressible boundary-layer flows. He applied the coordinate perturbation method and obtained asymptotic solutions for two limiting regions (see Fig. 1): for the large-slip (or quasi-frozen) region near the leading edge (I) and for the small-slip (or quasi-equilibrium) region far downstream of the leading edge (111). Zhao [18] used a similar series-expansion method and improved Singleton's analysis. However, these series solutions in the form of asymptotic expansions could provide only one term in addition to the frozen or equilibrium-flow values, owing to the complexity of the problem. Moreover, this solution does not provide any information on the boundary-layer development in the nonequilibrium transition region where the slip is moderate. A thorough understandi ng of the compress i b 1 e dusty-gas boundary-layer over the entire-flow region is important, since these flows have practical applications in many scientific and technical fields such as solid rocket exhaust nozzles, nuclear reactors with gas-solid feeds, ablation cooling, blast waves moving over the Earth's surface, conveying of powdered materials, fluidized bed and environmental pollution, as mentioned in Refs. [12, 19].

In the present paper, the behaviour of compressible, laminar boundary-layer flows of a dusty gas over a semi-infinite flat plate along the whole length of the plate is studied using a finite-difference methode The problem of two-phase suspension flows is solved in the framework of a mode 1 of two i nterpenetrat i ng and i nteract i ng cont i nuous medi a, whi ch is called a two-way coupling model or a two-fluid approach [20,21]. The following assumptions are made in this analysis: (l) The gas-particle mixture is a dilute system where the volume fraction of the particle phase is neglected. (2) The gas phase is a perfect gas. (3) The particles are spheres of uniform size without random kinetic motion. There are no mutual collisions or other interactions among the particles. (4) Only the drag and heat-transfer processes couple the particles to the gas. The momentum and energy exchange between the two phases can be cal cul ated from avail ab 1 e analytical solutions for the viscous flow field around a single sphere.

Finite-difference methods of solution of single-phase boundary-layer equations have been studied for many years. A review of this work is given in Ref. [22]. Flü"gge-Lotz and Blottner [23] developed an implicit difference techni que. They used a six-poi nt scheme for the momentum and energy equations and a four-point scheme for the continuity equation. This fi nite-difference procedure was appl i ed success fully to vari ous studi es of

pure-gas boundary-layer flows. However, in the dusty-gas case, the nature of the governing equations requi res some changes which result in considerable complexity. First, in addition to some new interaction terms in the conservation equations for the gas phase, there is an extra set of conservation equations for the particle phase. The partial differential equati ons for the gas phase are of second order, whil e those for the particle phase are of first order. Secondly, there is no correponding state equation for the particle phase, since the particle phase has no analog of flow pressure. In order to close the system of basic equations, the y-momentum equation for the partiele phase cannot be omitted as for the gas phase. Finally, the flow properties of the particles present quite different features i n di fferent flow regi ons. In the near 1 eadi ng-edge region, very large velocity slip and temperature defect between the two phases appear, whereas a quasi-equilibrium state can be reached in the far-downstream region where the flow profiles for the two phases are almost the same. Between these two regi ons, there i s a trans it i on regi on whi eh is characterized by a nonequilibrium flow. In general, the two phases in this regi on have moderate differences in velocity and temperature across the boundary layer. It is interesting to note that in the transition region, there exists a special position along the flat plate, which is defined as the critical point in this analysis. At the critical point, the tangential velocity of the particlesat the wall vanishes, that is, there is no slip between the part iel es and the gas. Th i sis due to the fact that in the two-phase boundary 1 ayer, the gas decel erates from its freestream velocity at the outer edge to zero at the wall and then the particles are retarded by the gas. The velocity of the particles at the wall may be reduced to zero provided the distance is long enough for the particles to adjust to the gas. Of course, equalization of the gas and particle velocities at the wall does not mean that the di sparity of the two phases has di ed out because across the boundary layer, equilibrium between the particles and the gas is still not attained. Nevertheless, it can be said that at the critical point, the two-phase system completes essentially the transition from a nonequilibrium flow to an equilibrium flow, since the equilibrium state is reached first on the surface of the plate and this process is continued until the two phases are in equilibrium across the whole boundary layer far downstream. As the particles are slowed down, the density of the particle phase near the wa" increases. When the particle velocity becomes zero, the particles tend to accumulate at the wa". In other words, deposition of the particles at the wa" may occur if there is no diffusion. Therefore, as discussed by Soo [24,25], there are two possible situations when the particle velocity decelerates to near zero:

(1) For large particles, their Brownian motion is neglected, the particles slowed down at the wall will deposit and form a sliding layer (or bed of particles). This compacted layer may build up or erode away, and even a steady equilibrium condition may be achieved when the shear stresses in this dense layer of the particulate matter and in the suspension mixture are equalized. The velocity at which such a layer moves depends on the materials and surfaces of the particles and the wa". Because of deposition of the particle phase, the density of the partieles at the

wa11 becomes very large. However, if the partic1e density is too high,

the present ana1ysis wil 1 fail since the assumptions concerning the

interactions between the two phases or among the particles in this paper

can be considered correct just for a di1ute gas-partic1e system.

(2) For sma11 particles, the Brownian diffusion is significant in the region

near the wa11, a1though the intensity of Brownian motion of the

particu1ate c10ud is usually small across the boundary 1ayer.

The

den s i ty of the pa rt i cl e phase at the wa 11 is then cont ro 11 ed by the

Brownian diffusion process.

The 1ayer of deposited particles may not

exist at all because the diffusion due to the Brownian motion prevents

the formation of a dense bed of particles.

It is shown in Soo's

analyses that, if the Schmidt number of the partic1e Brownian diffusion

is of order unity or 1ess, the who1e two-phase system behaves 1ike a

gaseous mixture and the density profiles for the partic1e phase reduces

to its origina1 one as in the freestream.

Soo studied on1y the

incompressib1e

boundary-1ayer

case.

Marb1e

[1,26J

treated

the

case of compressib1e f10ws and obtained a similar result.

Some other

studies [17,18J on the cornpressib1e boundary-1ayer flow of a

gas-partic1e mixture came to the same conc1usion.

Frorn the asymptotic

solution for the small-s1ip limit, it is found that the zeroth-order

approximation of the dimension1ess density for the partic1e phase is the

same as that for the gas phase.

This means that the 10ading ratio of

the part i cl es i s constant across the boundary 1 ayer and equa1 to its

origina1 va1ue in the freestream. Physica11y, in this quasi-equilibrium

flow region, the particles a1ways rernain attached or fixed to their

origina1 gas mass and move together with the gas. Then the gas-partic1e

mixture behaves 1ike a perfect gas with the modified properties. This

implies that the flow process in the sma11-s1ip region is main1y

diffusion-contro11ed for both the gas and particles. Therefore, in this

paper, it is assumed that, af ter the critica1 point, the partic1e

density is determined from the gas density and the loading ratio. Using

these cons i derat i ons, the fi nite-difference schemes for the dusty-gas

boundary-1ayer f10ws can be constructed.

In this ana1ysis, the

fi ni t e-d i ffe rence scheme deve10ped by F1tJgge-Lotz and B10ttner [22,23J

were emp10yed for the gas phase.

For the partic1e phase, a four-point

scheme was used. For comparison, the six-point scheme was used to solve

the x-momentum

and energy equations of the particles, emp10ying

additiona1 boundary conditions obtained from the compatibi1ity

conditions.

Af ter the critica1 point, very simp1e compatibi1ity

conditions for the tangentia1 velocity and temperature of the particles

can be derived:

at the wa11, the particles have the same velocity and

temperature as the gas.

With this finite-difference scheme, the flow properties of the dusty-gas

boundary 1ayer over the entire 1ength of a semi-infinite f1at-p1ate were

ca1cu1ated numerica11y.

The flow profiles of u, v and T for the two phases

are presented at different distances from the 1eading edge.

From these

results , it is shown that the boundary-1 ayer f10ws of a dusty gas have

different characteristics in the three distinct regions.

In the large-slip

region, the particles have a little deviation from their freestream uniform motion and then the differences in the flow quantities of the two phases are quite large. While slipping through the gas downstream, interaction between the two phases increases the gas velocity and temperature but decreases the particle velocity and temperature as well. Thus, in the transition region, the differences in the flow properties of the two phases are significantly reduced. Of course, the particles and the gas are still in nonequilibrium. In this region, the velocity slip and temperature defect are moderate compared with those in the other two limiting regions. In the small-slip region far downstream, the flow profiles for the particle phase become a 1 mos t i dent i ca 1 with those for the gas phase, that is, the two phases approach nearly equilibrium and the slip quantities are very smalle In fact, the only reason the part i cl es do not actually attai n the local gas velocity and temperature is that slip is induced along the gas streamlines by the gas retardation associated with thickening of the gas boundary layer. In addition to the flow profiles of u, v and T, some boundary-layer characteristic quantities of interest, i.e., the shear stress and heat-transfer rate at the wall and the displacement thickness, are calculated in this analysis. It is noted that owing to the presence of the particles, the shear stress and heat transfer increase while the displacement thickness decreases in the case of laminar boundary-layer flows, since the interaction between the gas and particles causes an increase in the gas velocity and temperature.

In this paper, the quasi-frozen flow properties in the near leading-edge region and the quasi-equilibrium flow properties in the far-downstream regi on were compared wi th the correspondi ng asymptot i c val ues [27]. The agreement was very good. For the finite-difference sol ut ion in the nonequilibrium transition region, it is found that the results are phys i ca lly reasonab 1 e • Al though it i snot poss i b 1 e at present to make any direct comparison between our finite-difference solution and other relevant results, since there are no experimental or other analytical data available for the nonequilibrium-flow region. Nevertheless, the fact that the fi nite-difference sol ut i on in the far-downstream regi on agrees quite well with the asymptotic small-slip solution provides confidence in the difference solution for this transition region, since the boundary-layer equations are parabolic, which is classified as a marching problem [30J. Thus, the solution procedure of finite difference begins with certain initial profiles at or near the leading edge, then through the large-slip region, the transition region, and finally ends in the small-slip region downstream. It is clear that the finite-difference solution for the small-slip region would not be correct if there were some mistakes in the difference solution for the transition region.

The numerical study of boundary-layer flows in dusty gases provides a good introduction to the dynamics of a two-phase system. The quasi-frozen flow, nonequilibrium flow and quasi-equilibrium flow are all encountered and analysed using the finite-difference methode The difference solution gives the complete and exact information about modifications of the boundary-layer

flows due to gas-particle interaction. Moreover, it provides a basis for the experimental investigation of dusty-gas boundary-layer flows.

2. MATHEMATICAL DESCR!PTION OF COMPRESS!BLE, LAM! NAR, DUSTY-GAS BOUNDARY-LAYER FLOWS

The basic boundary-layer equations for steady, two-dimensional,

I

compressible, laminar, dusty-gas flows over a flat plate are given by [27J Conti nuity:

~ p*u* + 0 p*v* = 0

~x* ~y*

(2.1)

Momenturn:

p* (u*

~u*

+ v*

~*)

=

0 (IJ.*

~*)

+ P! (u* - u*) u!, IJ.* 0

~x* ~y* ~y* ~y* . p p ~! IJ.!,

(2.2)

Energy:

p*c*(u* ~T* + v* ~T*)

P ~x* ~y*

=

~

(k* oT*) + IJ.*

(~*

) 2 + p* [( u* _ u*) 2

ê:tj* ~y* '0::1* p p

+ (v*p _ v*)

2]

u!. IJ.* 0 + _1_ P*pc*p(T*p _ T*) u!, IJ.* Nu

~!. IJ.! 3P r ~!, IJ.!

State:

p*

=

p*R*T* for the gas phase, and

Cont i nu ity :

x-momentum:

y-momentum:

~ p* u* + _0_ p* v*

=

0

~x* p p ~y* p p

ou* ~* u!. IJ.*

p* (u*

..::.:E.

+ v*

..::E.)

= -

P*p (u*p - u*) - - 0

p p ~x* p ~y* ~!. IJ.!. ~v* 'àtI* p*(u* ~ + v* ~) p p ~x* p ~y* u!, IJ.*

=

-pp*(v*p - v*) - - 0 ~!. IJ.!.

(2.3)

(2.4)

(2.5)

(2.6)

(2.7)

Energy:

aT*

oT*

u*

P*p

c~ (U*p ~

+

v*p

..::..E.)

=

-

_1_

P*p c*p (T*p - T*)

~ IJ.* Nu

(2 .8)

ax*

oy*

3P r

À.! IJ.!

for the partic1e phase.

In the equations above, the partic1e

velocity-equilibrium 1ength À.! is

p* d*

2

À*

=

s

u*

<XI

181J.:

<XI

The starred quantities in Eqs.

(2.1)-(2.8)

have dimensions.

The

independent variables are the space coordinates x* and y* which are parallel

and perpendicu1ar to the wa 11 , respective1y.

The dependent variables are

the density p*, the velocity components u* and v* and the temperature T* for

the gas phase as we11 as the corresponding quantities Pp' up' vp and T

p

for

the partic1e phase. For the f1at-plate boundary layer, the gas pressure p*

is constant and equa1 to its freestream value.

Hence, in the dusty-gas

boundary-1 ayer prob 1 em, there are ei ght simultaneous equat i ons with ei ght

unknowns so that this system is c10sed.

Of course, it is required that the

other physical quantities appearing in Eqs.

(2.1)-(2.8)

are known functions

of the flow variables.

The normalized drag coefficient D and the Nusselt

number Nu can be expressed in terms of the sl i p Reyno1 ds number and the

Prandtl number. Here, the normalized drag coefficient D is defined as

D

=

CD

CD

o

(2.9)

where CD and CD are the real drag coefficient for the flow situation under

consideration anB the drag coefficient from the Stokes re1ation.

In this

ana1ysis, on1y Stokes' relation is used and consequent1y D

=

1.0

and Nu

=

2.0.

Regarding the thermodynamic properties, the fo110wing assumptions are

made:

(1)

The specific heats for the gas and partic1e phases (cp and

c~)

are constant; (2) the Prandt1 number of the gas (Pr) is constant; (3) the

viscosity coefficient of the gas (IJ.*) has a power-law form with temperature.

Consequent1y, the expression for the gas viscosity is given by

(2.10)

where

w

is the power index for the viscosity coefficient.

It is advantageous to write the basic equations and relative expressions

in nondimensiona1 form before the numeri cal computations are performed. For

the i nvest i gat i on of two-phase boundary-1 ayer fl ows, it is conveni ent to

choose the ve10city-equilibrium 1ength

À.~

as the characteristic length, and

x

=

x* >,,*' 00 u

=

u* u*' 00

*

-V =

'i-

{Re 00' u* 00 v* -v

=..=E.

{Re P U

*

00' 00

*

P

=

..e.::..

_*'

Poo p* 0...

=

~ 'IJ _Poo

*'

T

=

T* T* 00 (2.11) T* Tp

=..:2-T*' 00

where the flow Reynolds number based on the partiele velocity-equilibrium 1 en gt h Re 00 is

Now, a nondimensional form of the boundary-layer equations results in

(2.12)

p ( u ou + v

~)

= loL 02u + (d loL • oT)

ru

+

Pp (

u p - u)!JI)

ox oy oy 2 dT öy 'óy

(2.13)

(2.14)

P

=

.!.

T

(2.15)

ru

~

Up

-=-e.

+

v

-=.2.

= -(

Up - u)

!JIJ

ax

p 'àj

(2.17)

av

ê:N

Up

~ +

v

-=-e.

= -(

v - v)

!JIJ

ax

p

ay

p

(2.18)

aTp aT n _ a

up

+

_vp

--L.. - - -

(Tp -

T) ~u

ax

'àj

3Pr

(2.19)

where the gas Eekert number Ee, the gas Prandtl number Pr and the ratio of

specifie heats for the two phases

a

are respeetively defined as

c* 1-1* P r =

---'Pl:.---k* '

c* a

=

.:.e.

c*

_s

The viseosity relation in nondimensional form can be written as

(2.20)

In order to obtain a unique solution to the problem, it is necessary to

satisfy the boundary conditions.

Inspecting the basic equations

(2.12)-(2.19), there are seven partial differential equations and two of

them are of seeond order.

Therefore, nine boundary eonditions should be

specified. If the partcle phase is in equilibrium with the gas ph ase in the

external flow, the nondimensional boundary eonditions are given by:

(1) At the wallof the fl at plate

u(x, 0)

=

0,

v(x, 0)

=

0,

T(x, 0)

=

_Tw'

(2.21)

vp(x,O)

=

0

(2) At the outer edge of the boundary 1 ayer

u(x,

CD)

=

1,

T( x,

CD)

=

_1,

(2.22)

up(x,

CD)

=

1,

Tp(x,

CD)

=

1,

Pp (x,

CD)

=

~

where the mass 1 oadi ng rat i

0

of the partieles

~

is

p*

~

=

PCD

..

Besides, owing to the parabolic character of boundary-layer equations, the initial profiles of the dependent variables are required across the boundary layer at some point xo:

u (x 0' y) = u o(y), v (x 0' y) = v o(y)

T (x 0' y) = T o(Y) , p(x 0' y) = po(y)

(2.23) up (x 0, y) = up 0 (y) , vp(xo, y) = _{vp 0 (y)}

T P (x 0' y) = T Po (y) , P_{p (x} 0' y) = P_{p o(y)}

At the initial position xo' the finite-difference solution procedure starts and then proceeds downstream.

3. FINITE-DIFFERENCE SCHEMES AND RESULTING FINITE-DIFFERENCE EQUATIONS The basic boundary-layer equations (2.12)-(2.19) with the boundary conditions (2.21)-(2.22) and the initial conditions (2.23) can be solved numerically using a finite-difference method. In this way, the parti al different i al equat i ons are approximated by fi nite-difference equations and the flow field is divided into a rectangular grid or mesh. Generally, either equal or unequal intervals can be used. In this report, equal intervals in the x-direction and unequal intervals in the y-direction were used in order to reduce the computation time (see Fig. 2). The step size in the y-direction was increased in a geometrie progression as

6.Yn

6.Yn-1

=

K

where K is a constant and it is set with a value slightly greater than unity. When

K

=

1.0, the unequal-interval mesh reduces to an equal-interval mesh. In the difference procedure, it is assumed that the flow quantities are known at the grid points in the column (m) and unknown at the grid points in the column (m+1). The computation starts stepwise downstream with the initial profiles.

When the finite-difference scheme is employed, the derivatives are repl aced by difference quotients. There are numerous ways of construct i ng difference quotients. For the sake of stability, implicit schemes, which can be six-point or four-point, are used in this analysis [28J.

For the momentum and energy equat i ons of the gas phase, a si x-poi nt difference scheme was used • . With this scheme, six grid points (m, n-1),

(m, n), (m, n+1), (m+1, n-1), (m+1, n) and (m+1, n+1) are involved. Any function w(x, y) and its derivatives are evaluated at a mid-point (m+e, n):

aw

=

L

_{Wm+_{1 n - Wm n}} ~x!J.x ' , 2a< _{[Wm+1,n+1 - {K+1}Wm+1,n + KWm+1,n-1]} {K+ I} ÓYn 2 + 2{1-e}K _{[Wm,n+1 - {K+1}Wm,n + KWm,n-1]} {K+1} ÓYn 2

{3.2}

{3.3}

{3.4}

where the wei ghti ng factor a can be chosen as any val ue between zero and unity. When a = 0.5, it reduces to the six-point Crank-Nicolson scheme where the truncation error is of order {ó><2} [29J. When a takes the value of zero or unity, it gives the full explicit or implicit scheme, respectively. The last two schemes involve only four grid points and have a truncation error of order {ó><}. But with respect to the variable y,all three schemes above are of second order {óy2}, since the cent ral difference formula is used for derivatives in the y-direction.

A four-point difference scheme is applied to the gas continuity equation. In this scheme, four grid points {m, n-1}, {m, n}, (m+1, n-1) and

{m+1, n} are included and all the values of the function w{x, y} and its derivatives are calculated at a mid-point (m+a, n - 1/2):

W

=

t

_{[a{Wm+1,n + Wm+1,n-1) + {l-a}(Wm,n + Wm,n-1}]} {3.5}

{3.6}

~W _

a (

)

1-

a (

)

-_~y - - - W_{óYn-1 m+ 1 n - W'} m+ 1 n -1 + _{, óYn-1'} Wm n - Wm n-1 _, {3.7} When a

=

0 or a = 1 in the above formulae, the truncation error in the x-direction is of order {ó><}. When e

=

0.5, the scheme is known as the four-point Wendroff scheme and the truncation error in the x-direction is of second order {ó><2} [29J. However, from experience in the present analysis,

..

gas phase. It was found that this osci1lation problem can be avoided by using e

=

0.75, which produces a discretization error of order (fl)(1.5). In the y-direction, this four-point scheme has a truncation error of order (lly 2) as in the si x-poi nt scheme.

For the part i cl e phase, due to the stabi 1 ity requi rement , the y-deri vati yes are approximated by backward difference quotients i nstead of the central difference quotients which are used in the above schemes for the gas. Then, another four-point difference scheme for the particles is constructed as follows. The function w{x, y) and its derivatives are estimated at the point (m+e, n). The derivatives, both in the x- and y-directions, are replaced by backward quotients.

Forthe x-momentum, energy and continuity equations of the particles, the grid points (m, n), (m, n+l), (m+l, n) and (m+l, n+l) are involved and the difference scheme is

w

=

_{eWm+l,n + {l- e)Wm,n} (3.8)

(3.9)

(3.10)

For the y-momentum equation, the grid points (m, n), (m, n-1), (m+1, n) and (m+1, n-1) are involved and the difference scheme takes another form for the y-derivative instead of Eq. (3.10)

(3.11) aw _

-ay

The function wand the x-derivative aw/ax have the same forms as Eqs. (3.8) and (3.9).

Similarly, when e is equal to zero or unity, which represents, respectively, the explicit or implicit scheme, the above four-point schemes (3.8)-(3.11) reduce to three-point schemes. These schemes have a truncation error of first order (lly). For stabi1ity consideration, as mentioned before, the value of e is chosen to be 0.75.

With the above formulae, the finite-difference equations for compressible, laminar boundary-layer flows of a dusty gas over a semi-infinite flat plate are given by:

(I)

The momentum equation of gas phase,

where

Al

₌

a n (pv - J.l1 T Y ) m+ S , n - C n ~ S, n

n

BI

₌

(PU)m+S,n _{+ an (K2-1)(pv - J.l'TY)m+S,n + cn(K+1)IJm+S,n}

n + St.x( ~~)m+1,n

..

Cl

₌

-anK 2( pv - J.l'TY)m+S,n - cnKl1n+s,n n

(2) The energy equation of gas phase,

2 2 2 2 An Tm+1,n+1 + Bn Tm+1,n + Cn Tm+1,n-1

=

Dn (n

=

2, 3, ••• N-1)

(3.13)

where + S t.x (-.L Ap!JNU)m+ 1 n

3Pr

'

D~

=

[(pu)m+S n - (l-S)t.x(-.L P...!JNulm n]Tm n

_,

_3Pr

_.~

_,

_,

2

For the coefficient Dn' the term

(Vp. -

v}2/Re"" is a small quantity and can

be neglected. In all the above coefficient expressions, some parameters are

gi ven by

an

=

9óx

_bn

=

(1-9}6X

(K+ 1) flYn

,

(K+ 1) flYn

cn

=

2 9K flx

,

_dn

=

2{1- 9}Kflx

(K+ 1) flYn 2

_{(K+ 1) flY n 2}

flW

_{m ,n}

=

W

_{m,n+1 + (K 2-1}W m,n - K 2W}

_m,n-

1

(3.14)

where the function W represents the flow properties such as velocity,

temperature, density, etc.

(3) The state equation of the gas phase:

1

P_m

+1,n

=

_.-:...-Tm

+

1,n

(n

=

1, 2, ••• ,N)

(4)

The continuity equation of the gas phase:

( pv ) m+ 1 ,n

=

(pv ) m+ 1 , n -1 - 1

~

9

[( pv) m , n - (pv ) m, n -1 ]

flYn-1 [

]

29f1x

(pu}m+1,n - (pu}m,n + (PU}m+1,n-1 - (pu}m,n-1

(5)

The x-momentum equation of the particle phase:

A3 u +B3 u

n Pm+1,n+1

n Pm+1,n

(n

=

1, 2, ••• N-1)

(3.15)

(n = 2, 3, ••• N)

(3.16)

(3.17)

where A_n 3

=

_-a~ v _p l!.Yn m+ a,n u -

a~

v +

a~(~)m+l,n

Pm+9,n _{l!.Yn Pm+9,n} _[(l-9)t.x v Ju +

[U

+ (1-9)t.x v

l!.Y_n Pm+a,n Pm,n+l Pm+9,n _{l!.Y n} Pm+a,n

(6) The y-momentum equation of the partiele phase,

where A4

v

+B4

v

n Pm+l,n n Pm+l,n-l v Pm+a,n (n =

2, 3, ••• , N)

(l-9)t.x v_p - (1-9)l!.x(IJI))m nJvp l!.Yn_l m+a,n ' m,n + [(1- 9) l!.x v Jv + t.x ( v

~

) m+ a , n l!.Yn-l Pm+9,n Pm,n-l

(7) The energy equation of the partiele phase,

A5 T + B5 T

n Pm+l,n+l n Pm+l,n (n

=

1, 2, ••• , N-l)

(3.18)

where

Up _{- at,x vp} + 9t,x (~ llNu ~+ 1 n m+a,n _{b.Yn m+ a,n} 3Pr '

_ [(l-

9) t,x V ]T + [U + (1-a) ÄX V

b.Y n Pm+ 9,n Pm,n+l Pm+a,n b.Yn Pm+ a,n

- (l-a)b.x(~ ~u'. _{n]T p} + b.x(~ TI-LNu' +a n

3Pr Jm, m,n 3Pr Jm,

(8) The continuity equation of the particle phase,

6 6 6 A p + B p

=

en

n Pm+l,n+l n Pm+l,n (n

=

1, 2, ••• , N-l) (3.20) where 29 _{up +} (1 29) _{- up} + 92t,x (v - 2v ) m+l,n m,n _{b.Y n} Pm+l,n+l Pm+l,n + a{l-9)t,x (v - 2v ) b.Y n Pm,n+l Pm,n

[{l-a)t,x v

]p

_{+ [(2a-l)u p + 2(l-9)u p}

b.Yn Pm+a,n Pm,n+l m+l,n m,n

_ a(I-9)t,x (v _ 2v ) - {l-9)2b.x (v - 2v )]p b.Y_n Pm+l,n+l Pm+l,n _{b.Y n} Pm,n+l Pm,n Pm,n

Using the finite-difference equations in the form (3.12)-(3.20), a stable and convergent numeri cal solution to the dusty-gas boundary-layer equations was obtained when x is smaller than xcri. After the critical point (x ) xcri)' quite simple compatibility conditions were derived for the part i cl e velocity and temperature. These condit i ons provided suppl ementa 1 boundary conditions at the wall so that the six-point scheme could be used for the x-momentum and energy equations of the partieles when x ) xcri. At

the wall (y

=

0), with the boundary conditions u

=

_{0 and vp}

=

0, Eqs. (2.17) and (2.19) become

(3.21)

(3.22)

These two equations, Eqs. (3.21) and (3.22), are termed as compatibility equations from which cornpatibility conditions can be derived. From Eq.

(3.21), it is known that, as x increases, the particle velocity at the wall decreases until it becornes zero. The position of the critical point is determi ned by

(3.23)

Af ter the critical point (X) xcri), the two phases have zero velocity at

the wall and then the drag vani shes (Ow

=

0). Thus, for x ) xcri' Eq. ( 3 .21) 1 eads to

u

=

0

Pw

(3.24)

Substituting Eq. (3.24) into Eq. (3.22) yields T

pW

=

Tw (3.25)

Equation (3.25) is valid for X ) xcri' too. Equations (3.24) and (3.25) mean that af ter the critical point the particles and gas are in equilibrium at the wall. Now, concerning the tangential velocity and ternperature of the particles, there exist two boundary conditions: one is at the wall and the other is at the outer edge of the boundary layer, as in the case of the gas phase. For the normal velocity, however, no such sirnple compatibility conditions, such as Eqs. (3.24) and (3.25), can be derived.

With the six-point scheme, the x-rnomentum and energy equations of the particle phase are replaced by the following difference equations:

( 1) Moment urn : A7 u +B7 u +C7 u n Pm+1,n+l n Pm+1,n n Pm+l,n-1 (n

=

2, 3, ••• ,

N-1)

_(3.26) where 7 A

=

a v n n

Pm+

e,n

-a

K2

v

n Pm+9,n

[ u p - (l-9) Ax ( IJI)) m n Ju p

m+S,n ' m,n

(2) Energy:

A a T + BaT + CaT

= Dna

n Pm+1,n+1 n Pm+1,n n Pm+1,n-1 (n

=

2, 3, ••• , N-1)

_(3.27) where A 8

=

a v n n Pm+S,n Up _{+ a (K2_1)v p} + 9Ax(~ llNu' +1 n m+S,n n m+ S,n 3Pr

..m,

-a K 2v n Pm+S,n

[Up - {l-S)Ax(~ llNu)m nJTp

m+S,n 3 P r ' m,n b n Pm+ S,n Pm

v

~T + t.x ~ TllNu lm+s n

3Pr '

Therefore, when x ;> xçri' the di fference equations are composed of Eqs. (3.12), (3.13), (3.15), (3.16), (3.18), (3.26) and (3.27) with the assumption that the particle density is determined by

Pp

=

~p. The detailed derivations of all the finite-difference equations above are given in Appendices A and B.

It is noted that the boundary-layer equations (2.12)-(2.19) are a coupled nonlinear partial-differential system. To avoid the coupling and nonlinearity, in the process of discretizing every conservation equation, only one corresponding variable appears as an unknown in the resulting di fference equat i on and the difference expressi ons for the products of the unknown variables, functions or derivatives are chosen such that the unknown variables appear linearly in the products. This procedure leads to a linear system of algebraic equations which are not coupled. As pointed out by Blottner [22] the coupling between the equations results in a tridiagonal matrix which is somewhat more complicated to solve than in the uncoupled case. In addition, for linear algebraic equations, there are several effective methods of solution available. For example, the Thomas algorithm

is a very powerful and convenient technique to solve the linear equations with the tridiagonal matrix of the coefficients.

4. METHODS OF SOLUTION OF THE FINITE-DIFFERENCE EQUATIONS

The methods

of solution depend upon

the characters of the

finite-difference equations.

For the six-point scheme, the resulting

difference equations constitute the system of simultaneous algebraic

equations with a tridiagonal matrix of the coefficients.

Advantage can be

taken of this tridiagonal form of the coefficient matrix to solve the

algebraic equations by use of the Thomas algorithm [30].

With the Thomas

algorithm, the solution is obtained by

W~+l,n

=

Gi _ Ei

n

n

W~+1,n+1

(4.1)

Fi

_n

where

Gi

Ei

_n

₌

_A~

_,

_Fi

_n

₌

Si

_n

ei n-1

_{n -.-,}

Ei

Gi

_n

₌

Di

_n

ei n-1

n-·-F~_l

F~_l

The recurrence relation, Eq. (4.1), can be used to solve difference

equations (3.12) and (3.13), or (3.26) and (3.27). Correspondingly,

W

1

=

U,

W 2

=

T,

W 7

=

Up and W

8

=

T p.

Sy using the following procedure, the flow profiles can be determined:

(1)

~ith

the boundary conditions at the wall, calculate quantities

E~, F~,

G~

from the wall towards the outer edge.

(2) With the bOIJndary conditions at the outer edge, calculate the flow

properties

w~

from the outer edge towards the wall.

Af ter the gas temperature Tm+1 n is known, the gas density Am+1 n can be

calculated directly from the statè equation (3.15).

Then, starting at the

wall and using the gas continuity equation (3.16), the normal velocity

vm+1,n can be obtained.

When x

<

xcri' the differential equations for the particle phase are

discretized using the four-point difference schemes.

The methods for

solving the difference equations (3.17) to (3.2Ó) are not difficult. Af ter

starting either at the wall (only for the y-momentum equation) or at the

outer edge of the boundary layer, the calculations proceed consecutively

from one grid point to another in recursion, until the whole boundary layer

has been traversed:

i

W

m+

1

,n

ei

=

n _(i

₌

₄₎

_(4.3)

where W

3

₌

_{up' W}

4

₌

_{V ,}

_W

5

₌

Tp and

W6 = P

_{p• Af ter the critical point}

(x ) xcri)' the x-momenfum and energy equations of the particle phase can be

discretized with the six-point scheme and then solved by the Thomas

algorithm, as for the gas phase.

Before sol vi ng the difference equati ons numeri cally usi ng the methods

described above, some considerations are required:

(1) How to evaluate the coefficient matrix elements;

(2) How to give the boundary conditions in a form suitable for the numeri cal

computation;

(3) How to obtain the initial profiles;

(4) How to determine the value of xcri' the nondimensional coordinate for

the critical point.

Fi rst, the fi nite-difference

equa~

i

on~

can

b~

sol ved provi ded that the

values of all the coefficients

A~, B~, e~

and

O~

are known.

However, from

the expressions of these coefficients,

it

is seen that they depend on

unknown values of the variables at the grid line (m

+

1), since the

difference scheme is an implicit one.

This difficulty can be surmounted by

using an iteration procedure. Of course, the iteration technique increases

the computation time very much. The other way to overcome the difficulty ;s

to use a l;near;zation approximation:

the quantit;es appearing in the

coeffic;ents are evaluated at the prev;ous grid line (m) if these quantit;es

are still unknown at the gr;d line (m

+

1). Otherwise, they take on their

updated values. In this analysis, this linearizat;on approach was employed,

since it is easier to program and requires less computer storage.

Of

course, it is less accurate compared with the iteration approach. However,

satisfactory accuracy can be achieved by reducing the step size.

Second, in order to solve the system of simultaneous algebraic equations

at every new grid line (m

+

1), it is necessary to have the boundary

conditions in a suitable forma

In this analysis this is straightforward,

since there are no derivatives involved in the boundary conditions of the

Dirichlet type. In the finite-difference scheme, the boundary conditions at

the wall, Eq. (2.21), are written as

um+1,1

=

0,

vm+1,1

=

0,

(4.4)

v

=

0

Similarly, the boundary conditions at the outer edge, Eq. (2.22), are given in the form, um+ 1 ,N

=

1, Tm+1,N

=

1

(4.5)

u

=

1 Pm+ 1, N '

T

Pm+ 1 ,N

=

1,

Here, concerni n9 the outer-edge boundary condit i ons, another diffi culty arises: How to select the number N, the maximum value for the number of grid points at the column (m+1). In other words, in the computation process, it is required to know how far to calculate the flow variables across the boundary layer. In order to guarantee th at the value of N used represents the freestream condition, one can specify a large number ~ax for the grid points at the last grid line (mmax) far downstream, where the computati on termi nates. For all the previ ous gri d 1 i nes (0 < rn < "\na x ), the same number N

rnax is used to defi ne the outer edge of the boundary layer. This method is direct but inefficient, since it needs more computation time. There is another approach in which a special value of N ;s chosen for a given grid line (m). In the latter method, Nis determined as follows

[31J:

(1) It is assumed that ~+1 _{at the new grid line (m + I) is equal to Nm, the} number of grid points at the previous line (m).

(2) The finite-difference equations are solved with the assumed ~+1 and the values of the flow quantities at the last two consecutive grid points, Wm+1,N-1 and Wm+l,N' can be found.

(3) The difference between W_{m+l N-1 and Wm+1 N is compared with a certain} small quantity e:. If the cond~tion of smooth conjugation

I

Wm+ 1 ,N - Wm+ 1 ,N -1

I (

e: (4.6)

is satisfied, the selection of Nm+1 is correct and the computation can proceed to the next step.

(4) If the condition (4.6) is not fulfilled, it is required to assume Nm+1

=

N

m + 1 and then to obtain the new values of Wm+1 N-l and Wm+1 N. If the condition (4.6) is not fulfilled again, it is neèessary to inèrease

~+1

with unity, and so forth, until the smooth-conjugation condition is satisfied. With this method, the number of grid points across the boundary layer varies as the thickness of the boundary layer increases. Next, in order to initiate the computation, the initial flow profiles must be given. In most pure-gas boundary-layer studies, the initial profiles are obtained from similarity solutions. For dusty-gas boundary-layer flows, however, no analogous similarity solutions exist. In previous work on finite-difference solutions for incompressible dusty-gas boundary-layer equations, the initial profiles were specified in two ways:

(1) The Blasius similarity profiles were chosen for the gas phase and uniform profiles for the particle phase [14].

(2) The initial profiles were obtained by using an integral method [13, 15].

It is we" known that all the integral methods for boundary-layer analysis do not attempt to satisfy the basic equations at every point; instead, they guess or assume a suitable expression for the velocity and temperature profiles and satisfy the boundary-1ayer equations only on an average extended over the thickness of the boundary 1ayer. In general , the initia1 profil es obta i ned from i ntegral methods are quite approximate. From the studies of the behaviour of dusty-gas boundary-layer flows near the leading edge

[17, 18, 27],

it is also known that the simi1arity profiles for the gas phase and the uniform profiles for the particle phase are the zeroth-order approximation in the large-slip region. The zeroth-order asymptotic profiles are physically reasonable and were tested in this analysis. More accurate initial profiles up to the first order can be obtained from the asymptotic large-slip solution to the compressible laminar boundary-layer equations for the gas-particle flow over a semi-infinite flat plate by using a series-expansion method [27]. Thus, it is suggested here to employ the first-order asymptotic solution as the initial profiles. However, in this approach, it is required to obtain the asymptotic solution first and then to solve the difference equations starting at a given initial position which may be very near the leading edge but cannot be exactly at the leading edge, since the asymptotic solution involves a singularity at the leading edge. Wu [32] once proposed the following type of initial profiles at the leading edge in pure-gas cases:

(1) The tangential velocity u* and temperature T* have their freestream values at all the grid points across the boundary layer except at the wa 11.

(2) At the wall, the tangential velocity u~ is zero and the temperature T~ corresponds to the wall temperature.

(3) The normal velocity v* is assumed zero at all the grid points.

Clearly, this method is very advantageous for starting numeri cal computation without any preliminary calculations of initial profiles. F1û'gge-Lotz and Blottner [23] studied the possibility of using the Wu-type initial profiles in pure-gas boundary-layer cases and concluded that the Wu-type of initial profiles can give reasonable results if the proper mesh sizes are chosen. In the dusty-gas case, similar initial profiles can be set up by using Wu-type profiles for the gas phase and uniform profiles for the particle phase. Soo [24] applied such initial profiles to his analysis of incompressible, laminar boundary-layer flow of a dilute particulate suspensi on. The profil es are termed here as the extended Wu-type. In thi s report, the above three different types, i.e., the first-order

asymptotic-type, the zeroth-order asymptotic-type and the extended Wu-type, were used respectively as the initial profiles in order to compare them.

As mentioned before, the asymptotic types of initial profiles can be obtained from the large-slip solution [27J. However, because of the different notations, certain relations must be established with the asymptotic large-slip solution as follows:

x

=

{x)asy' y = {/2x Tl)asy v

=

(1

(T}J - f) )asy' 12x T

=

(T)asy' Tp

=

(T p) , asy Pp

=

~(pp) asy p

=

(p) asy {4.7}

Similarly, a set of relations can be written for connecting the finite-difference solution with the asymptotic small-slip solution. In this analysis, nondimensional slip quantities for the particle phase are defined as

(4.8)

The corresponding relations are

x

=

{x)asy' y -_ 1

-==

(/2x - Tl) asy Il+~

(4.9)

v

_s

=

1

( s )

,

11+~ I2x asy

These expressi ons are useful when compari ng the fi nite-difference sol uti on with the asymptotic small-slip solution. In Eqs. {4.7} and (4.9), the subscript lasyl denotes the asymptotic solution.

Finally, as pointed out earlier, af ter the critical point (x ) xcri)' it is assumed that the particle density is equal to the gas density times the mass loading ratio of the particles. It is equivalent to assume that there is no accumulation of particles on the surface of the plate and the flow is then mainly diffusion-controlled for the particle phase as well as for the

gas phase.

In addition, af ter the critical point, quite simple

compatibility conditions (3.24) and (3.25) are val id. The value of xcri can

be determined from the compatibility equation (3.21) and the condition

(3.23).

Equation (3.21) is an ordinary differential equation and the

solution upw(x) can be obtained numerically or analytically. For instance,

in the case of the Stokes relation (0

=

1.0), Eq. (3.21) can be integrated

analytically as

(4.10)

where upw(O)

=

1 is taken as the initial condition at x

=

O.

From the

condition uPw(Xcri)

=

0, Eq. (3.23), the critical value xcri can be

determined, say, for the Stokes case

=

1 (4.11)

Pw

when T

w

=

0.5 and

w

=

0.5, xcri

=

/2 or

x~ri

=

/2 À!,.

If

p~

=

2.5 g/cm

3 ,

d*

=

10 Ilm, T!,

=

300 K (or Il!,

=

1.80x10-~

NS/m

2)

and u!.

=

500 miS, the

typical values of relaxation parameters are obtained as:

À!.

=

0.386m and

x~ri =

0.546 m.

5. RELATIONS FOR SHEAR STRESS, HEAT TRANSFER, AND DISPLACEMENT THICKNESS

Once the gas

fl

ow profil es across the boundary

1

ayer are determi ned,

some boundary-layer characteristics of practical interest can be calculated.

The important quantities describing the behaviour of boundary-layer flows

are shear stress at the wall, heat-transfer rate at the wall and

displacement thickness. They are given in dimensional form as

(1) Shear stress at the wall:

't:'A.: = u'A.: ~~*)

w rw

'óy* w

(2) Heat-transfer rate at the wall:

q*

=

-k* ('óT*)

w

w 'óy* w

(3) Oisplacement thickness:

CD ö* =

J

(1 - p*u* )dy*

p!, u!, o

(5.1)

(5.2)

(5.3)

The corresponding nondimensional characteristic quantities are defined as

1;* 1;

=

w

w P* u* CD CD 2

(5.4)

(5.5)

ö

=

ö* !Re '11.* CD CD

(5.6)

Substituting the nondimensional transformations

(2.11)

into the expressions

(5.1)-(5.3),

the

nondimensional

boundary-layer

characteristics,

Eqs.

(5.4)-(5.6), can be written as

~

=

~ (<XJ )

oy w

(5.7) •

_tJw

_(oT)

qw

-

-Pr Ec

oy w

(5.8)

CD Ö

=

_J

(1 -

pu)dy

(5.9) 0

For numeri cal computations, it is necessary to express the above

relations in finite-difference form.

By means of polynomial fitting, the

gas velocity u and temperature T near the wall may be expressed with

sufficient accuracy as

(5.10)

(5.11)

Taking the derivatives of the above variables with respect to y and setting

y

= 0,

the formulae for shear stress and heat-transfer rate at the wall are

obtained as

•

_tJw

Pr Ec bT

The values of bu and bT can be determined by evaluating Eqs.

(5.10)

and

(5.11)

at the four grid points nearest the wall and solving the resulting

equations. Then the shear stress and heat-transfer rate can be given by the

following expressions (see Appendix C):

'"w

=

~

(K2 + K + 1) [

u3

+

u4

]

K2~Y1

u2 - K(K + 1)

K(K 2 + K + 1) 2

(5.12)

• _~

_{(K2 + K + 1) [(T}

2 _

Tl) _ T3 - Tl +

T4 - Tl

qw

-

-

]

Pr Ec

_K2~Y1

_{K(K + 1)}

_{K(K2+K+ 1)}

₂ (5.13)

where the subscripts 1,

2,

3,

4

denote the four grid points nearest the wa11

and uI

=

Uw

= 0,

Tl

=

Tw.

To ca1cu1ate the nondimensiona1 displacement thickness

ö,

a three-point

difference formu1a of integration was used. The formu1a can be app1ied to a

nonequidistant step size [33]

N-1

ti

L

!I

n -1 [3K + 2

F

+ 3K + 1 F

2

6

K + 1 n-1

K

n

(5.14)

ö

=

where F

=

1 -

pu, and the subscri pts n-1, n and n+l represent three

consecutive grid points at a given section.

6. COMPUTER PROGRAM

The computer program FOBLEP for sol vi ng boundary-1 ayer equat i ons of a

dusty gas over a semi-infinite flat p1ate was written in Fortran 1anguage on

a Perkin-E1mer computer system at UTIAS (see Appendix

0).

On1y the main features of the ca1cu1ation procedure will be reviewed

here.

A rectangu1 ar-gri d system i ndi cated in Fi g. 2 was adopted, with the

m-1ines in the y-direction (i.e., normal to the p1ate) and the n-1ines in

the x-direction (i.e., parallel to the p1ate). The y-axis is designated as

the initia1 1ine for the m-1ines and the x-axis for the n-1ines.

In the

finite-difference procedure, the flow profiles a10ng some m-1ine (say, the

initia1 1ine m

=

0)

are known and the flow parameters a10ng the (m+1)-line

have to be determined.

Given be10w are the main steps in the computation procedure:

(1) Compute Um+1,n.

Using the known solution at the m-1ine (linearized conditions) and the

boundary conditions on the (m+1)-line, the new tangentia1 velocity of

the gas at all the grid points of the (m+1) line, um+1,n' are

ca1cu1ated.

(2) Test for the outer edge of the boundary 1ayer whi1e computing u. Af ter

the boundary 1ayer has been traversed, the two consecutive va1ues of u,

um+1 N-1 and um+1 N' are compared to see if the difference between them

is tess than some small positive number

E.

The va1ue of

€:

is

determined by the desired accuracy of the computation.