A random element method for a thermal boundary layer

-II

2%!3

ATHEMATICS AND COMPUTERS N SIMULATION

83 %

ELSEVIER

Mathematics and Computers in Simulation

A random

element

method

37 (1994) 527-538

for a thermal

boundary

layer

Gary A. Sod Tulane University, Department of Mathematics, New Orleans, LA 70118, USA

Abstract A new grid-free method for approximating incompressible thermal boundary layers is introduced. The computational elements are segments of vortex sheets which are heated. The method solves the general Prandtl boundary layer equations. More efficient methods of treating the diffusion and source terms are incorporated. The method is applied to the parallel flow past a heated horizontal flat plate and a flow created entirely by buoyancy forces (free-convection flow) on a heated vertical flat plate.

1. Introduction

In [l] Chorin presented a grid-free numerical method that approximated the incompressible boundary layer equations in the absence of temperature variations. The algorithm approximated the vortex sheets within the boundary layer by segments of vortex sheets. The equations solved are the Prandtl boundary layer equations. In this paper we extend the vortex sheet algorithm to treat the more general incompressible boundary layer equations which allow for temperature variations and bouyancy affects (using the Boussinesq approximation). The equations of motion and energy are coupled by the buoyancy term in the momentum equation and by the advection and friction terms in the energy equation. The diffusion process is modeled by a grid free method that uses a judious combination of random walks. The sourse terms are modeled by grid-free divided differences that introduce either no or a negligible amount of numerical diffusion. 2. Principle of the method

The boundary layer equations can be written in the form (e.g., Schlichting [2]) a,u + (u. Y)U = vay%4 - a,p/p

+g,P(T-

T,),

a,u + ayu = 0,

0378.4754/94/$07.00 0 1994 Elsevier Science SSDI 0378.4754(94)00099-6

B.V. AI1 rights reserved

528

G.A. Sod /Mathematics

and Computers in Simulation 37 (1994) 527-538

where u = (u, u) is the velocity, u is tangential to the boundary and u is normal to the boundary, x is the spatial coordinate tangential to the boundary and y is the spatial coordinate normal to the boundary, p is the density, p is the pressure, T is the temperature, v is the kinematic viscosity, k is the thermal diffusivity, and cp is the specific heat at constant pressure. The term g,P(T - T,) represents a body force due to bouyancy forces caused by temperature differences, that is, thermal expansion, where g, is the gravitational force in the x direction and p is the coefficient of thermal expansion. We assume a wall is located at y = G and that the fluid occupies the half-space defined by y 2 0. The boundary conditions are _u=O

at

y=O,

(la)

u(x, 03) = K(% t>

(lb)

and T=T,,,(x,

t)

at

y=O,

(2a)

T(x, m) = T,(x, t).

(2b) Introduce four dimensionless quantities, Re = U,L/v, the Reynolds number where L is a characteristic length, Pr = v/k, the Prandtl number, EC = U,*/c,AT, the Eckert number where AT = T, - T,, and Gr = g/3L3AT/v2, the Grashof number where g is the gravitational constant. This gives rise to the dimensionless form of the equations a,u + (_u- _O)u = Reel azu -

a,p

+

Gr Re-* T, (3)

a,u + a,u = 0, a,T + (u. y)T = (Pr Re))‘a:T

+ EC Re-‘($u)*.

Introduce the vorticity (within the boundary layer) 5 = -i3,,u and the heat flux Taking the derivative with respect to y of Eqs. (3) and (4) we obtain

(4) q =

a,T.

a,( + (44. p)& = Re -l a,‘,$+ Gr Re-* d,,T,

(5)

a,u + a,u = 0,

(6)

at4 + (g

- y)q

=

(Pr Re))‘azg

+&u

- @,T - 2E, Rep1 @,“u,

(7)

where ay($,u)* = - 2@j$. Integration of the vorticity field 5 yields U(X,

Y) = u,(x,

t) + /“6(x,

s> ds

Y

Using this along with the continuity Eq. (6) we obtain u(x, Y) =

-axLyu(x, Z) dz

(8)

G.A. Sod /Mathematics

=

-a,[

Urn@,

= -a,

t)~

and Computers in Simulation 37 (1994) 527-538

- /,i[m+,

Um(x, t)y -/os(x,

529

s) ds dz]

s) min(s,

(9)

y) ds].

1

Thus if the vorticity field 5 is known, the velocity field u can be determined. Similarly, we may integrate heat flux 4 to obtain T(x,

Y) = X(x,

t> - /=c+,

s) ds.

(10)

Y

Let At denote the time step. At a time t = nAt consider a collection of N segments of vortex sheets with intensities ,.$/’ and heat fluxes qlF, i = 1,. . . , N. The Si are straight line segments parallel to the x-axis having equal length h and center (x1!, yl). The intensity and heat flux are uniformly distributed over the length of the segment Si. The difference in u across the segment Si is [r, that is, 5; = -(U(x;,

yin+) - U(XY, y;-))

and the difference qi” =

in T across the segment

41,

that is,

T(xi”, y;+) - T(xi”, y;-).

Evaluating (8) at the center the velocity of this segment u(x;,

Si is

yi”) = Um(xi”, ndt)

of the i-th segment

+ jm[(x;,

(xr,

y:) gives the tangential

s) ds.

component

of

(8’)

Y:

Similarly,

evaluating

(10) at (x:,

T(xi”, y;) = T,(x:,

ndt)

yr) gives the temperature

- /b(xf,

of the i-th segment

s) ds.

(10’)

Y:

Let dn=

’

(1-

Ix:-xj”I/h,

if lxr-x;I

\o,

then the discrete


otherwise analogue

u; = Um(x;, ndt)

of the integrals c

- is,! -

j#i, Y:

&%f;

(8’) and (lo’>, respectively,

(11)

>Y:

and T” = Tm(x;, nAt)

-

c j#i, YY>Yl

q;di”.

(12)

530

G.A. Sod /Mathematics

and Computers

in Simulation

37 (1994) 527-538

The number of interactions among the segments is small since for a segment Si to influence a given segment Si it must intersect the semi-infinite strip Ri = {(x, y) Ixi” - h/2 G x G x1 + h/2}, where x1! -h/2 and xl? + h/2 are the ends of the segment Sj. This term dy is a weight which corresponds to the length of the portion on of the segment Sj which lies in the strip Ri.If Sj does not intersect the strip Ri then dr = 0. Below the sheet Si the velocity is increased by 51 and the temperature is increased by sr. The velocity of a segment at a height y is increased by all segments of sheets which lie above y. Similarly, the temperature of a segment at a height y is increased by all segments of sheets which lie above y. Use a centered difference approximation to a, evaluated at the center (xr, y,“> of the i-th segment to obtain the normal component of velocity of the i-th segment, ‘.:=-~[~~~(~~+~,y)dy-~~iijx:‘-q,y)dy]+O(h’) h

= -i

U, xy+T,nAt li

- ( Lml(

( X;

+

i,

S)

i

-U,

(

h xF-?>nAt

min(s, yr) ds - iat(

ii

yi

x1 - s, s) min(s, YF) ds)] + 0(h2).

The integral l$nu(xr + h/2, y) dy is the flux of u across a segment of the vertical line at xi” + h/2 between 0 and yr. Each segment which intersects this vertical line segment affects the velocity of this line and hence the flux. If a segment Sj lies above yn, that is, y,? > yr then Sj increases the velocity of the enters vertical segment from 0 to yCy and hence the flux is increased by -yn[,?. If a segment Sj lies below yr, that is, y,:’< yn then Sj increases the velocity of the vertical line segment by sj below Sj and hence increases the flux by Y,?[~!.This is summarized as follows

jti )I c S;d;y,+’

/h

,

where

and y,Fi” = min( yn, yf ) . We now will describe the algorithm. Operator splitting is used to write Eqs. (5)-(7) in terms of advection, diffusion and source terms. The following equations state that the segments are

G.A. Sod /Mathematics

advected 45 +

by the velocity

and Computers in Simulation

531

37 (1994) 527-538

field (8)~(9),

(u . _v)5 = 0,

at4 + (u. c= -a+,

_~)cI=

0,

4 = a,T, a,u + a,u = 0. These equations can be solved in a deterministic manner. Let x(t) denote the position i-th segment as a function of time and consider the ordinary differential equations

of the

Given the position of i-th segment at time t = nk, its position at time t = (n + l)k is obtained by solving the ordinary differential equation using Euler’s method (which has O(At) truncation error) x;+’

=x3 + Ami”,

y;+’ =y; +Atu; where u” and U: are determined using (11) and (13). The diffusion process shall be approximated by a random walk of the segments. However Eqs. (5) and (7a) may have different diffusion coefficients, while the heat flux and vorticity are assigned to the same sheets. In the case where Pr = 1 no further diffusion takes place. This problem can be remedied by splitting the diffusion term in (5) or (7a) into the sum of two diffusion terms atq = (Pr Re)-l$q a,5 = Re-’ Consider,

= Reel

$9 + (1 - Pr)(Pr

a,‘[ = (Pr Re))‘a:[

without

+ (Pr - l)(Pr

a loss of generality,

a,[ = Re-’

a,‘,_$,

a,q = Re-’

a;+

Re))‘azq, Re)-‘a:t,

the case where

if Pr < 1, if Pr > 1.

(14)

Pr < 1,

These can be solved deterministically by choosing a random variable qi drawn from a gaussian distribution with mean 0 and variance 2At/Re. This yields the algorithm x;+’

=x;

+Atu;,

w4

y;+’ =y; + At@ + qi. The remaining

diffusion

a,q = (1 - Pr)(Pr is treated

by a separate

ow term

Re))‘azq random

(16) walk.

532

G.A. Sod /Mathematics

and Computers

in Simulation

37 (1994) 527-538

For the i-th segment, i = 1,. . . , N, divide the heat flux qf among a number, M, of phantom segments of length h, with center (xl?“, yn+’ ) (as obtained by (15)), and heat flux q*!‘/M. To solve Eq. (16) we shall let each of those phantom segments undergo P (an integer 2 1) random walks (rather than a single random walk) in the y direction. Given a value of p and m, where p = 1,. . . , P and m = 1, 2,. . . , M, let qi(rn, p) denote a random variable drawn from a gaussion distribution with mean 0 and variance 2At(l - Pr)(Pr Re)-l/P, the position of the phantom segments (xl+‘(m), y~“(m>>, is obtained as follows y;+’ (m,p)=yin+‘(m,p-l)+rl,(m,

P),

p=l,...,P,

where yF+‘(rn, 0) = yF+l, the y component position of the i-th segment given by (15b). Define x:+‘(m) =x)+l (the x component of the i-th segment given by (Ha)) and y:+‘(m) = yn+‘(m,

P).

We now use the phantom sheets to compute a new heat flux for the i-th segment at location due to the phantom (x;+l, y;+‘), i = l)...) N. To this end, we compute a temperature segments at the segment locations (x:+‘, yf+‘) using Eq. (12)

Tin+’ = T,(x;+‘, ndt) - m=l f

F

(17)

($)G(m)

y;+‘(m)>y:+’

where -x;+‘(m)

I/h,

if I xi”+’ -x;+‘(m)

I
otherwise. We now must use this to modify the heat flux of the segments. First observe that T;+l= 7(x;+‘, y n+ ‘-). To see this consider (12) where y is so large that there are no segments above this value of y, then T = T,(x, nAtI. Suppose the N segments are ordered so that y;+’ a~,“+’ 2 ..a >yE+l. Then T(x, y;+’ ‘) = T,(x, ndt) and T(x, y) = T,(x, rdt) q;dy+’ for yz+l
T(x,

y) = T,(x,

nAt) - qFdr+’ = T(x,

Y;“-).

Y-Y;+'

Similar results hold for T(x, yr+’ -). The following algorithm is inspired by this observation. With the ordering y;+’ >yi+l> *. * >yi+l, q13’ = T(x;+‘,

y;+l+)

-

7):“+l,

(18)

where Tin+’ is given by (17) and i-l T(x;+~,

y:+l+)

= T@+‘,

nAt) -

c qj”&fj”+‘.

(19)

j=l

It should be noted that the ordering of the segments by a sorting algorithm will also minimize the number of decisions involved in carrying out various summations.

533

G.A. Sod /Mathematics and Computers in Simulation 37 (1994) 527-538

It remains to describe the treatment differential equation

of the source terms in (5) and (7). From (51, the partial

a,[ = Gr Re-*a,T is discretized in space using a centered divided difference ‘Tjx,

Y+$)-+,

y-q)

+o(Ay2)

a,[ = Gr Re-* AY I

\

The resulting ordinary differential intensity ti of the N segments

equation

is solved, using Euler’s method,

to modify the

AY

~P+~,yl’+l-~

I

5;”

= 5: +

AtGr Re-* AY i

From Eq. (71, the partial differential cretized in space

equation

++I, -+-;, y)

T(x+;,

Y)

44

=4

dtq = qd,u - &T

- 2Ec Re-’ &y*u

y) -T(x-;,

is dis-

y)

h

h

‘U(X. Y+;)

-2u(x~y)+u(~~y-:)’

- 2 EC Re-l[

AY

+O(h2)+O(Ay2)

*

i-i 2

\

I

The resulting ordinary differential equation is solved, using Euler’s method, to modify the heat flux qi of the N segments q;+’

zqi”J

xl+‘+

+Atq;,’

h

.,+I

yl?+’

2,

h

-

2,

Y;+’

/h ))

x;+‘+

;,

yin+’ -T ) (

XT+‘,

yr+’

+

q

x;+‘s, yr” -

2U(Xl+‘,

/h yr+’

Xr+‘,

yin+’

-

9

/AY*.

We must describe how Ay is chosen. We want the divided differences in the y direction to involve points that lie within the boundary layer, so we introduce a scale of Re-‘I* which approximately represents the thickness of the boundary layer. Define Ay = A yRe- ‘I2 where

534

G.A. Sod /Mathematics

and Computers in Simulation 37 (1994) 527-538

Ay= O(h). Thus the 0(Ay2) truncation error of the centered differences is O(h2 Re-‘1 = O(Re-‘1 which is much smaller than the standard deviation of the random walk dw. If ,;+I - Ay/2 < 0, so that the point falls outside the physical domain, we set AY

(using the boundary condition (2a)) and u(xp+‘, yr+’ - Ay/2) = 0 (using the tangential component of the boundary condition (la)).

3. Vorticity and heat flux creation Using (11) we see that boundary condition (lb) is satisfied and (13) we see that the normal component on the “no-slip” condition (la) is satisfied. Similarly using (12) we see that boundary condition (2b) is satisfied. However, the tangential component of the “no-slip” condition (la) and the boundary condition (2a) are not satisfied. Divide the boundary into segments of length k with centers (Xi, 0). At time t = nk, using (11) we can compute u(Xi, Of, ndt)

= U,(X, ndt)

The “no slip” condition (X,, 0) of strength &p = u(Xi, O+, ndt)

requires

- ~.$%I;.

u(Xi, O-, nAt) = 0. So a vortex sheet must be created

- u(Xi, 0-, ndt)

at

= u(Xi, O+, ndt)

per unit length. Similarly, using (10) T(i,

O+, nAt)

= Tm(Xi, nAt)

- &l”d;

while the boundary condition (2a) requires that T(ii, of the vortex sheet created at (Xi, 0) is 4; = T(ii,

O+, ndt)

- T(&,

0-, nAt)

O-, ndt)

= T(.Ei, O+, ndt)

= Tw(ii, ndt). - T&,

So the heat flux

ndt)

per unit length. The vortex sheet at (Xi, 0) is broken up into mj elements, each having intensity c//rni and heat flux 4”/rni where m, is chosen so that 141: I /mi < qmax and I ,$r I/mj < 5”” simultaneously. Consider the case in which U, = 0, that is, free-convection flow, so that ti = 0 for the sheets. Solving the ordinary differential equation stemming from the source term in (5) a,< = Gr Rep2 ayT, results in the production of vorticity. Suppose 5” = 0 then [“‘l= AtGr Rep2 a,Ti” is the strength of the i-th sheet. If gLp+’> 5”” then the sheet is subdivided into some number of

G.A. Sod /Mathematics

Fig. 1. Boundary

and Computers

Layer of a parallel

in Simulation

flow past a heated

37 (1994) 527-538

horizontal

535

flat plate.

sheets each with center (xi, yi) and strength less than 5”“” and sum of strengths equal to [(p+‘. The heat flux 41 is divided equally among these sheets so that the total heat flux has not been altered. Each of these subdivided sheets will undergo a random walk and will each have different locations at subsequent times.

4. Application to flow past a heated flat plate Consider a semi-infinite flat plate located on the positive x-axis. The fluid at temperature T, occupies the half plane y > 0. At time t = 0 the plate is instantaneously and uniformly heated to temperature T, and the fluid is impulsively set in motion with velocity U, (see Fig. 11. We shall compare the results of our method to known solutions (see Schlichting [2]). Heat due to friction is important only if EC = O(1). In the case of a gas, this is equivalent to requiring that the temperature rise due to adiabatic compression be of the same order of magnitude as the difference in temperature between the plate and the fluid. As a test we consider the case when the physical parameters are Re = 4.56 X 106, Pr = 0.71 (air), EC = 7.186, Gr = 0 (b uoyancy neglected), L = 1, U, = 1, T, = 1, and T, = 0. This represents the case of a heated flat plate where the affects of frictional heating are important, that is, EC > 0. The numerical parameters are At = 0.2, h = 0.2, A y = 0.2, 5”“” = 0.025, qrnu = 0.025, M = 4, and P = 5. In Fig. 2, we display the temperature profiles taken at the point x = 0.7, which are averaged over 10 time steps for 16 < t G 18. The total number of segments at t = 18 was 1348. The results are in excellent agreement with the similarity solutions superimposed on the figures. Let (Ie 112denote the /,-norm of the relative error in the temperature profile taken at the point x = 0.7 which is averaged over 10 steps for 16 < t G 18 and let N denote the total number of segments present at time t = 18. In Table 1 we present a comparison of the effect of Ay on

Table

aL Ilell2

N

1 0.2

0.1

0.05

0.36~10~’ 1397

0.30x 10-l 1866

0.29x 10-l 2542

G.A. Sod /Mathematics

536

and Computers in Simulation 37 (1994) 527-538

1.4 1.2

A

Boundwy LeyerTheory Computed Velues

1.0

0.8

T-T_ T,-Tc.,

0.6

0.4

0.2

Fig. 2. Temperature profiles in a forced convection horizontal flat plate in a parallel stream.

boundary

layer on a heated (E, > 0) and a cooled (E, < 0)

the error and the total number of segments. All other parameters are as above. We see that 11e II2 is not very sensitive to changes in dy. This is consistent which the order of magnitude of the difference approximations to the source terms. Note, however, the number of segments increases as A y decreases.

x

T

1

9

Thermal Boundary

Layer

TCU

Velocity Boundary

Layer

Fig. 3. Boundary layer of a natural convection flow a vertical flat plate with T, > T,.

G.A. Sod/Mathematics

and Computers

in Simulation 37 (1994) 527-538

531

1.0 -Boundary A Computed

LayerTheory Values

0.8

0.6

T-T, T,-T 0.4

0.2

rl=Y

Fig. 4. Velocity profile in a natural

5. Application

to natural-convextion

convection

boundary

boundary

layer on a heated

vertical flat plate.

layers

Consider a semi-infinite vertical flat plate placed in a fluid at rest. Assume the coordinate system is oriented in such a way that the plate lies along the x-axis (see Fig. 3). The fluid at

4

rl=Y

Fig. 5. Temperature

profile

in a natural

convection

g(T,-T,)

I/-

boundary

“2X

layer on a heated

vertical flat plate.

G.A. Sod /Mathematics

538

and Computers in Simulation 37 (1994) 527-538

temperature T, = 0 occupies the half plane q = 0. Assume the plate is heated to a temperature T, > T,, in which case the fluid adjacent to the plate is heated, becomes lighter, and rises, resulting in the formation of a velocity and thermal boundary layers (see Fig. 4). There is no free stream, so there is no reference velocity. One possible reference velocity is

which contains the buoyant terms. With this reference velocity we can define the Reynolds number Re=

KL

-

V

=

vigs(T,-L3’2 =

(p’

V

so that Gr Rep2 = 1. We consider the case when friction is negligible, that is, E, = 0. The parameters are Pr = 0.71, Re = 1.68 X 105, T, = 1, T, = 0, and L = 1. The numerical paramters are At = 0.2, h = 0.2, Ay= 0.2, 5”“” = 0.025, q”‘” = 0.025, M = 4 and P = 5. In Figs. 4 and 5 we display the temperature and the profile of the tangential component of velocity profile, respectively, taken at the point x = 0.7 and averaged over 10 time-steps for 16 < t < 18. The results are in excellent agreement with the similarity solutions superimposed on the figures.

References [l] A.J. Chorin, Vortex sheet approximations of boundary [2] H. Schlichting, Boundary-Layer Theory (McGraw-Hill,

layers, J. Comp. Phys. 27 (1978) 428. New York, 1968, 6th Ed.).