Question

Computers & Fluids Vol. 17, No. 4, pp. 537-553, 1989 01M5-7930/89 $3.00 + 0.00 Printed in Great Britain. All rights reserved Copyright © 1989 Maxwell Pergamon Macmillan pie

S T E A D Y L A M I N A R F L O W I N A 9 0 D E G R E E

P L A N A R B R A N C H

R. E. HAYES, K. NANDAKUMAR and H. NASR-EL-DIN Department o f Chemical Engineering, 536 Chemical-Engineering Building, University o f Alberta,

Edmonton, C a n a d a T6G 2G6

(Received 21 June 1988; m revised form 6 December 1988)

A b s t r a c t - - T h e flow characteristics of a Newtonian fluid in a two-dimensional, planar, right angled Tee b r a n c h are studied over a range of inlet Reynolds number of 10-800 by solving the Navier-Stokes equations using a finite element discretization. The effects of the branch length and the grid size on the interior flow field are examined to assess the accuracy of the solutions. In one case the computed velocity field is compared with the Laser Doppler anemometry measurements available in the literature and excellent agreement has been obtained. The computed velocity field is believed to be accurate within about 5%. Results are presented for two types of experimentally realizable boundary conditions--viz, equal exit pressure at the outlet of each branch and specified flow split between the branches. For the case of equal exit pressures the fractional flow in the main duct increases with increasing Reynolds n u m b e r and the flow characteristics in the side b r a n c h become akin to that in a cavity. F o r the case of specified flow split, the number, size and strength of the recirculation zones increase as more fluid is forced to go into the side branch. The length of the side branch appears to have very little influence on the interior flow field, particularly at higher Reynolds number. This observation is rationalized as being due to the parabolized approximation becoming more valid at higher Reynolds numbers. The critical Reynolds n u m b e r at which the first recirculation zone appears in the side branch increases with increasing fractional flow in the side b r a n c h and with decreasing side branch width.

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

Branching flows of fluids and suspensions occur in a number o f biomechanical and engineering applications. A number o f variations in flow configurations are possible and a few common ones that have been studied in the literature are shown in Table 1. For a given geometry such as a 90 ° Tee, a 45 ° branch or a symmetric Y-branch, a number of flow configurations such as combining or dividing flows are of interest. The literature cited in Table 1 is by no means exhaustive; it merely serves to illustrate the scope o f the problem and will be used to develop a coherent picture o f the general flow characteristics in such geometries.

From a biomechanical point o f view, the formation of atherosclerotic plaque and thrombi in human cardiovascular system is believed [1] to be strongly influenced by the fluid mechanical considerations such as boundary layer separation, formation of recirculation or stagnation zones. With this perspective, various aspects of this problem have been studied by a number o f investigators. Lynn, Fox and Ross [2] studied the 2-D steady flow characteristics in a symmetric Y-branch using a finite difference approximation. More recently, Ehrlich and Friedman [3] have computed the 2-D flows through regions approximating human aortic bifurcation using boundary- fitted coordinates. Bramley and Sloan [4] have also used boundary-fitted coordinates with particular emphasis on the effect of sharp or smooth corners on the flow separation. The pulsatile flow characteristics in a symmetric, 2-D Y-branch have been examined in a numerical study by O'Brien et al. [5] while flow visualization and velocity measurements o f pulsatile flows in 3-D ducts have been presented in Siouffi et al. [6] and Fukushima e t al. [7]. Siouffi et al. [6] used a symmetric Y-branch with a branch angle o f 60 ° and a duct o f rectangular cross-section with an aspect ratio o f 0.8. Karino e t al. [8] carried out an experimental study using a model 3 mm diameter right angled T-junction with square and rounded corners. They used a particle tracking method with neutrally buoyant polystyrene beads and high speed photography to determine the velocity fields. Liepsch et al. [9] used a 2-D 90 ° branch to measure the detailed velocity field using LDA, while Kried e t al. [10] obtained similar data in a 3-D channel with circular cross-section and a dividing flow configuration. Popp and Sallet [11] present both flow visualization results and velocity measure- ments for single- and two-phase flows in a 2-D 90 ° branch.

537

538 R. E. HAVES e t a l .

Table 1. Branching flows studied previously

Flow Geometry characteristics Reference Nature of work

Steady 2-D,'~ Lynn et al. [2] Numerical, finite difference I. --,<] Newtonian J Bramley and Dennis [27] Numerical, finite difference

Ehrlich and Friedman [3] Numerical. boundary-fitted coord. Bramley and Sloan [4] Numerical, boundary-fitted coord.

Pulsatile 2-D,'~ O'Brien e t al. [5] Numerical, finite-difference 2. - - ~ Newtonian J

Pulsatile 3-D,'~ Fukushima e t al. [7] Experimental, flow visualization 3. - - ~ Newtonian J Sioutfi e t al. [6] Experimental, Velocity measurements 4. ~ Turbulent 3-D~, Pollard [14] Numerical, k ~ model

Newtonian j Heat and Mass transfer 5. ~ Steady 3-D'~ Karino et al. [8] Experimental, particle tracking

Newtonian J 6. ~X~ Steady 2-D'~ Liepsch e t al. [9] Experimental and numerical

Newtonian J 7. ~ Steady 3-D'~ Kried e t al. [I0] Experimentak L D A

Newtonian J Steady 2-D, 1- "~ Popp and Sallet [11] Experimental, LDA 8. and 2-phase flow J

9. ~ - ~ Steady 3-D'~ Sparrow and Kemink [ 2 8 ] Experimental, heat transfer Newtonian J

10. , / " ~ Steady 2-D Ukeguchi e t al. [29] Transonic flow numerical

1 I . . , / , Steady 3-D "~ Ku and Liepsch [30] Experimental, LDV Non-Newtonian J Flow visualization

There are two essential features common to this type of flow problem and both are caused by a sudden change in the flow direction. They are (i) flow separation and recirculation zone developing in the streamwise or main flow direction and (ii) secondary velocities developing in a plane normal to the main flow direction. The first one is a dominant feature in all branching flows and can be suppressed only under exceptional circumstances [6]. This has been computed in most of the 2-D Y- and T-branches (e.g. Bramley and Sloan [4] for 2-D, Y-branch; Liepsch e t al. [9] for a 2-D, T-branch) and has been observed in both 2-D and 3-D branching geometries (Poop and Sallet [11] for a 2-D branch; Karino e t al. [8] for a 3-D branch). Such flow separation has been observed and computed even in a smooth 90 ° bend [12] with a short radius of curvature. For a 90 ° branch, the amount of fluid withdrawn in the side branch can also induce a weak recirculation zone in the main duct itself [8]. The second feature of developing secondary velocities is of course dominant in coiled ducts as the centrifugal force is sustained in that geometry, but in branches the fluid experiences this force only near the branch region where the fluid is forced to turn into the branch. Under proper flow conditions the secondary flow can develop, but will decay in the downstream direction. In experiments simulating 2-D conditions (i.e. large aspect ratio rectangular ducts) [9, 11] the secondary flow can be expected to form only beyond a threshold Reynolds number (more appropriately Dean number) [13]. In shorter 3-D ducts the secondary flow is always present even at low Reynolds number and its magnitude can be up to 30% of the main flow [6]. All of the numerical simulations cited in Table 1 (with the exception of Pollard [14]) are limited to 2-D flow geometry. Pollard [14] has studied the turbulent 3-D flow with heat and mass transfer in a 90 ° branch with the k-E model.

Phase separation in branches is another aspect that has received wide attention in the literature for its obvious importance in power plants and wet steam distribution in enhanced oil recovery. Available results on phase separation of gas-liquid flow in branches are reviewed by Lahey [15]. Phase separation of solids from suspensions is of concern in developing solids sampling probes [16] and in ensuring even distribution of coal through manifolds in coal liquefaction process [17]. Phase separation of solids in blood has been the focus of an experimental study by Bugliarello and Hsiao [18]. The effect of the single phase flow field on the particle trajectories will be the subject of a follow up study. In the present study we focus attention on the steady, laminar flow of a Newtonian fluid in a right angled branch.

In the present work we have undertaken an extensive parametric study of the 2-D branching flow through a right angled branch. The steady, 2-D form of the Navier-Stokes equations for an incompressible, Newtonian fluid are solved using a finite element discretization. The flow

Steady laminar flow in a 90 degree planar branch

a--~2 = o

d

539

~=o \ ! . =o°. sp.o.,e0 J \ / v ~ = v 2 = °

F W - 1 v~ = v z = 0 ds

o I

x 2

v2 L ,-

o ~ " | a * = o " ~ l O x 2

I v l = O I v Z = p a r a b o l i c p r o f i l e

Fig. 1. Geometry, coordinate system and the finite element discretization.

characteristics are functions o f the following parameters: (i) R e y n o l d s number, R e = WUc/v, where W is the width o f the main branch, Uc is the inlet centreline velocity and v is the kinematic viscosity, (ii) width o f the side branch, ds, and (iii) specified fraction o f flow in the main branch o r specified exit pressures. F r o m the point o f numerical accuracy, the effects o f grid refinement a n d d o w n s t r e a m b o u n d a r y conditions (in particular the length) o n the u p s t r e a m flow profile are also discussed.

G O V E R N I N G E Q U A T I O N S

The geometry o f the Tee, the c o o r d i n a t e system a n d the finite element grid used in this study are s h o w n in Fig. 1. The equations governing steady viscous incompressible flow are:

1 ( V . V ) V = - - V P + v V 2 V ( 1 )

p

V. V = 0 (2)

C.A.F. 17/4---42

540 R . E . HAYES e t al.

where V denotes the velocity vector, P is the pressure, and p and v are the fluid density and kinematic viscosity, respectively. The weak variational formulation o f eqns (1) and (2) in a domain

b o u n d e d by F is [19]:

1 R(V, P ) = v(VV, V~/) + (V. VV, ~/) - - (P, V.r/) = 0 I3)

P

(~p, V. V) = 0 (4)

where R(V, P ) is the residual and q and ~o are Galerkin test functions in the appropriate spaces o f V and P. ( . , . ) denotes a scalar product, defined as:

f" (Z, ~o) = Jo Zq~

df~

The b o u n d a r y integral term arising from this formulation is [19]:

It q .(Pn - vVV.n)dF,

(5)

(6)

where n is the unit outward normal. F o r a specified velocity on F, this term vanishes. In the present problem this will occur at the entrance (OH), where a parabolic profile is imposed a n d on the walls (OA, BCD, E F H ) where the no slip condition is used. At the main outlet (AB) and the branch outlet (DE) one o f two physically realizable b o u n d a r y conditions can be specified--viz. (i) specified pressure at b o t h the main and side branch exits (e.g. both exits exposed to atmosphere) or (ii) specified flow split (defined as the ratio o f flow rate at the exit o f main branch to the inlet flow rate). In the first case the known pressure appears as d a t a through the b o u n d a r y integral eqn (6) leaving the velocity gradients as the natural b o u n d a r y conditions.

9, VI ~ V2 - - 0 on AB ( 7 a )

~x2 ~x2

0 V ~ = ~ 3 ~ = 0 on D E (7b) 9x~ ~x~

These constraints correspond to a fully developed state at the exit and, in principle, they can only be realized with sufficiently long lengths o f ducts. Increasing Reynolds number would tend to increase the development length required to reach the above state. On the other hand, the parabolization approximation o f the Navier-Stokes equation becomes more realistic with increas- ing Reynolds number, implying t h a t the effect of the above b o u n d a r y conditions on the upstream flow field should become less significant with increasing Reynolds number. In fact for a fully parabolized formulation, Gunzburger et al. [19] show that only P at the outflow needs to be specified and eqn (7) c a n n o t even be imposed on the solution. A l t h o u g h our present formulation is fully elliptic, thus requiring the b o u n d a r y conditions (7), the above mentioned effect has been observed in our simulations and the fact that their effect on the interior flow field is minimal has been verified for select cases by comparing results like velocity, stream function and vorticity profiles from short and long length ducts at a few Reynolds numbers. F o r the second case, when a specified flow split is to be realised, a parabolic velocity profile corresponding to the known split is specified in AB, hence the b o u n d a r y integral vanishes. On the b o u n d a r y DE, however, only the pressure is specified again, leaving eqn (7b) as the natural b o u n d a r y conditions for the second case also. Since, for an incompressible fluid, the absolute value o f pressure is immaterial, its specification on D E merely serves as a reference value for the entire flow field.

The discretization o f eqns (3)--(4), (6) is performed using the e n r i c h e d - - q u a d r a t i c v e l o c i t y - - linear pressure triangular element o f Crouzeix and Raviart [20]. This is a standard quadratic triangular element to which a bubble function has been added at the barycentre. At the c o m p u t a t i o n a l level, the central node is eliminated from the global matrix t h r o u g h the use o f static condensation, leaving 12 d.f. per element. The integration is performed using a six-point Gaussian quadrature.

Steady laminar flow in a 90 degree planar branch 541

After discretization the matrix problem can be written

Rm(V, P ) = A (V)V - B r P - F = 0 (8)

BV = 0 ( 9 )

where V and P are the velocity and pressure nodal values respectively, A (V) is the nonlinear Navier-Stokes matrix and B is the divergence matrix. F is a solicitation vector which incorporates the b o u n d a r y conditions. The above f o r m u l a t i o n is referred to as the mixed f o r m a t i o n since it includes the pressure explicitly as an u n k n o w n . The u n k n o w n pressure term can be eliminated t h r o u g h the use o f the penalty f o r m u l a t i o n [21], which gives the matrix system:

Rm(V, P ) = A,(V)V - F = 0

BV = 0

P = _ I B v

(lO)

( l l )

(12) E

where E is the penalty parameter and:

1 A,(V) = A(V) + - BTB. (13)

E

The non-linear system represented by eqns (10)-(13) is solved using the modified New- t o n - R a p h s o n m e t h o d [22] in which the Jacobian is computed at the first iteration only. This Jacobian is then used for all subsequent iterations, which greatly reduces the c o m p u t a t i o n a l cost. In order to obtain convergence at higher Reynolds number it is necessary to have a good initial estimate o f the velocity field. Thus, a simple c o n t i n u a t i o n in Reynolds n u m b e r is used by providing a previously converged velocity field at a lower Reynolds as an initial estimate for obtaining solutions at higher Reynolds numbers.

Once a converged velocity field is obtained, the stream function distribution is c o m p u t e d using the following defining equations for stream function, ~b, and vorticity, co.

v, = ~ , v , = - ~x--7 co = ~x~ ~ x , (14)

The stream function equation is:

V2~ = o . (15)

The corresponding variational form is:

Rm(~k) = (V$, Vr/) + (co, t/) = 0. (16)

The b o u n d a r y conditions for the stream function equation are:

~k = 0 on E F G , ~k = Qoh on O A , ~k = Qd* on B C D (17a)

d---~ = 0 on O H , AB and D E , (17b)

where Qoh is the inlet flow rate t h r o u g h b o u n d a r y O H a n d Qde is the outlet flow rate t h r o u g h DE. F o r the case o f equal exit pressures, the stream function value on the b o u n d a r y B C D is n o t k n o w n a priori, but can be determined by integrating the converged velocity field across DE.

The interior pressure field can be c o m p u t e d from the Uzawa algorithm [23] using the converged velocity field and the linear pressure elements. It can also be obtained as a solution to the following Poisson e q u a t i o n t

_ 2[ov, ov, ov, OvL-l, v ' P = kOx, Ox, e x , Fx, j (18)

tThe authors are grateful to o n e of the reviewers for pointing out the importance of pressure field in biomvchanies and for suggesting an alternate procedure for computing it.

542 R . E . HAYES et aL

where the right hand side is c o m p u t e d from a converged velocity field. The b o u n d a r y values tbr the pressure field on the solid walls are natural and those on the entrance and exit are obtained from the converged solution. Since the pressure discretization within each element is linear, the pressure gradients c o m p u t e d from the U z a w a algorithm are discontinuous across element boundaries. In that sense the pressure profiles c o m p u t e d from the Poisson equation (18) is to be preferred.

R E S U L T S A N D D I S C U S S O N

It is essential to select the proper c o m p u t a t i o n a l parameters to ensure the accuracy o f the c o m p u t e d results. In the present case these parameters are: (i) the grid discretization (N x N) a n d the length, L o f the duct d o w n s t r e a m o f the branch. In order to assess the effect o f these parameters on the accuracy o f the solution three different lengths (L = 3, 6, 26) and grid sizes (5 × 5, 7 x 7, l0 x 10) were tested over a range o f Reynolds number. A grid parameter o f (7 × 7), shown in Fig. 1, implies 98 triangular elements in a unit square o f O I G H . The lengths are scaled with respect to the entrance width, W. U p o n refining the grid from 5 x 5 to l0 x l0 for the case o f equal exit pressures with L = 3, the change in predicted values o f the split were less than 2% for inlet Reynolds n u m b e r o f up to 600. A duct length o f L = 6 was used t h r o u g h o u t and a few selected simulations were carried out with L -- 3 and L = 26 to assess the effect o f branch length on the flow field, particularly near the exit o f the side branch. The simulations for L = 6 using 7 × 7 grid were carried out on an FPS 164 and the simulations for L = 26 with 7 × 7 grid were carried out on a C Y B E R 205. It t o o k a b o u t six N e w t o n iterations to reduce the relative errors below i0 --3.

Figure 2 shows the streamlines for a flow split o f 0.5 and Reynolds numbers o f 300 and 700. Here the Reynolds number, Re = WUc/v, is defined based on the centre line velocity at the inlet, Uc and the entrance b r a n c h width, W. In Fig. 2(a) the streamlines obtained using L = 3 and 6 are superimposed over the c o m m o n flow d o m a i n o f xE (0, 4) and yE(0, 6), a l t h o u g h for the longer duct the actual flow d o m a i n is x~(0, 7) and yE(0, 9). It should also be pointed out that in the longer duct

6

(a)

Re = 300

°o ;

6i

o6

(b)

Re = 700

-0.0:

J 2 3 4 5 6 7

Fig. 2. Effect of branch length on the flow field. (a) Streamlines for Re = 300, split = 0.5 computed with L = 3 (---) and L = 6 (--). (b) Streamlines for Re = 700, split = 0.5 computed with L = 6 (---) and L = 26 (--). Note that in both the cases the flow corresponding to the larger length reached a fully

developed state at the exit.

Steady l a m i n a r flow in a 90 degree p l a n a r b r a n c h 543

9

I

6

5 A

3 r

2

0

(a)

Re = 4 1 0

Split = O . 5 6

. . . . . . . . . i . . . . . . . . . i . . . . . . . . . l . . . . . . . . . i . . . . . . . . . i . . . . . . . . . i

0 1 2 3 4 5 6 7' 0

i (b)

Fig. 3. C o m p a r i s o n o f c o m p u t e d velocity profiles ( - - ) w i t h the m e a s u r e m e n t s o f Liepsch e t al. [9] for Re = 418, split = 0.56. (a) Velocity profiles. (b) Streamlines. C o m p u t a t i o n a l p a r a m e t e r s are L = 6 a n d grid = (7 x 7). All t h e profiles are n o r m a l i z e d w i t h respect to the inlet average velocity. T h i s c o r r e s p o n d s

to Re = 496 in Ref. [19] b a s e d o n the average inlet velocity a n d h y d r a u l i c diameter.

the flow did reach a fully developed state at both the exits. In the Fig. 2(b) the results obtained using L = 6 and 26 are superimposed in a similar manner over the common flow domain and once again the flow corresponding to the longer duct with L = 26 resulted in a fully developed state at the exits. In both cases there is a recirculation zone crossing the exit boundary of the shorter duct, the existence of which is confirmed by the simulation using the longer duct. Observe that in Fig. 2(b) there is a second weak recirculation zone appearing near the top plate of the side branch, and the result from L = 6 is able to predict its existence. These results confirm our expectation that the influence of the natural boundary conditions (7) on the interior flow field becomes less important with increasing Reynolds number. Also observe the appearance o f a recirculation zone in the main branch in both cases. Flow visualization results of Popp and Sallet [22] for a 2-D flow and that o f Karino et al. [8] for a 3-D flow are in qualitative agreement with the above observations on the recirculation zones.

A measure of quantitative agreement is seen in Fig. 3 where the computed velocity profiles are compared with the Laser Doppler anemometry measurements of Liepsch et al. [9]. Excellent agreement is obtained for Re = 418 using the computational parameters o f L = 6 and a grid of (7 x 7). This corresponds to Reynolds number of 496 in Ref. [9] based on the average inlet velocity and hydraulic diameter for a rectangular duct of 10 x 80 mm. The agreement between the measured velocity profiles from Ref. [9] and our computed values is better than a similar comparison shown in Ref. [9], which is perhaps due to the coarser grid used in the finite difference calculations in Ref. [9]. The grid sensitivity test and the comparison with the experimentally measured profiles suggest that the use of the computational parameters of L = 6 and a grid of (7 x 7) is adequate in the following parametric study. Although we obtained converged results for Re as high as 1400 with a grid of (7 x 7), the results are presented for Reynolds number of only up to about 800 as the reliability of the results at higher Re decreases for both physical and computational reasons. Physically, one can expect the centrifugal forces to induce instabilities in the 2-D flow beyond a threshold Re, making the flow structure 3-D, even for a planar geometry. Numerically a much finer grid is required to resolve the sharp gradients expected at higher Reynolds numbers.

The streamlines are shown in Fig. 4 for the case o f equal exit pressures and d, = 1.0 at various Reynolds numbers. Similar results are shown in Fig. 5 for Re = 300 and various widths of the side branch, d,. Flow separation from the lower wall o f the side branch occurs at higher Reynolds number, but the critical Reynolds number at which the wake beings to form is a function of the side branch width. Since the fraction o f fluid entering the side branch is sharply decreased with

544 R . E . HAYES e t al.

:l-!,II ,., Re;lo i :tI/! l ( b ) Re = 200

3

I

t i

- - 4

6 '

5

4

3

2

1

0 o

Re = 600

Fig. 4. Streamline c o n t o u r s for the case o f equal exit pressures. Side b r a n c h width, ds = 1. In each case a fully developed state is reached at the exit as s h o w n by the c o n s t a n t level o f vorticity at the exit in

Fig. 7(a).

increasing Reynolds number, the flow characteristics in the side branch become akin to that of flow past a cavity with the side branch serving as the cavity. The depth o f penetration o f the recirculation zone is not significantly affected by the Reynolds number, but the width increases to span almost the entire cross section of the side branch. The strength and size of the recirculation zone is significantly affected by the side branch width, d,, as seen in Fig. 5.

The fractional flow in the main branch is shown as a function of Reynolds number for various d, in Fig. 6. For a side branch width of d, = 1.0, the split increases rapidly from about 0.5 at the creeping flow limit to about 0.95 at Re = 600. At each fixed Reynolds number, the split increases with decreasing side branch width. Also shown in this figure are three data points for the flow split computed with a grid of (5 x 5) to provide a measure of grid sensitivity. All these results were obtained with a length of L = 3 since a fully developed state is reached within this length as can be verified from Figs 4 and 5.

The vorticity profile (which for a 2-D flow is also the shear stress profile) on the lower wall of the side branch is shown in Fig. 7 for various Reynolds numbers and side branch widths. These profiles appear very similar to those computed by Bramley and Sloan [4] for a symmetric Y-branch.

Steady laminar flow in a 90 degree p l a n a r b r a n c h 545

0 0

(a)

ds=0.5

/ ' ~ - - - - - - - - ~ 0.045 ~ ll ~ - - 0.015

m

2 3

! (b)

ds=0.75

----~---- 0.06 ~

3 4

(c) ds=1.00

0 0 0 6 ~

3 4

Fig. 5. Streamline c o n t o u r s for the case o f equal pressures a n d Re = 300 for various d,. Even for Re as high as 600, a fully developed state is reached at the exit as s h o w n by the c o n s t a n t level o f vorticity at

the exit in Fig. 7.

1.1 0

L 1.0 ,a

=0 o,g

~ o . 8

0

0.7

0

:~ 0 . 8 0

0.5

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

900 Tee. E q u a l e x i t p r e s s u r e s . Branch l e n g t h = 3

/ ~ 1.00 (to=to) j - - 0 . 7 5

/ ---- 0.50 • ' ' 0.20

o 1 . 0 0 ( 8 = e )

0 100 200 300 400 500 600

R e y n o l d s n u m b e r

Fig. 6. Effect o f increasing Reynolds n u m b e r o n the flow split for the case o f equal exit pressures for various side b r a n c h widths.

5 4 6 R . E . HAYES e l al.

0.10

0 T: 0 >

0 . 0 8

0 . 0 8

0 . 0 4

0 , 0 2

0 . 0 0

- 0 . 0 2

- 0 . 0 4 0

Equal exit pressures, d s ~ 1.0

-- Re=lO

--- Re=200

. . R e = 6 0 0

/ / J ~ ~ ~ ~ / < l \ r e a t t a c h m e n t c o r n e r • ' p o i n t

1 2 3 4

D i s t a n c e a l o n g l o w e r wall (x2--2.0) of b r a n c h

0,10

0.08

0.06 > ,

0 . 0 4

~0 0.02 >

0.00

- 0 , 0 2

- 0 . 0 4

. . . . , . . . . , • . , , . . . .

(b) Equal exit pressures, d s = 0.75

- - Re = 50

--- Re = 200

• Re = 600

er \k.J .i"

. . . . i . . . .

0 1 2 3 4

D i s t a n c e a l o n g l o w e r wall (x2=2.125) of b r a n c h

0 . 1 0 (C) Equal exit pres res. d s = 0.50

0 . 0 8 - - R e = 5 0

0.06 ~ - - - Re = 200 • R e = 6 0 0

0 . 0 4 .o

"~0 0.02 . . . .

0.00 "7

- 0 . 0 2

- 0 , 0 4 0 1 2 3 4

D i s t a n c e a l o n g l o w e r wall (x2=2.250) of b r a n c h F i g . 7. V o r t i c i t y d i s t r i b u t i o n o n t h e l o w e r w a l l o f s i d e b r a n c h . A f u l l y d e v e l o p e d s t a t e is r e a c h e d i n e a c h

c a s e a s i n d i c a t e d b y t h e c o n s t a n t v a l u e o f t h e v o r t i c i t y a t t h e exit.

The vorticity increases sharply n e a r the corner. Th e extent o f the recirculation zo n e is determined by the sign change in the vorticity along the wall. A t Re = 10 in Fig. 7a, there is n o separation. T h e extent o f the recirculation zone decreases with decreasing d,. Th e vorticity reaches a c o n s t a n t value n e a r the exit o f the d u c t in all the cases, confirming t h at the flow reaches a fully developed state in each case even for Reynolds n u m b e r as high as 600. T h e stream fu n ct i o n profiles shown in Fig. 5 f o r Re = 300 and various b r a n c h widths also indicate t h at a fully developed state is reached at the exit f o r each case.

T h e r e a t t a c h m e n t point was determined f o r a n u m b e r o f intermediate Reynolds n u m b ers by locating the place o f sign change o f vorticity a n d these d a t a are shown in Fig. 8 for three different side b r a n c h widths. In each case the d e p t h o f p e n e t r a t i o n o f the recirculation zo n e into the side b r a n c h reaches a c o n s t a n t value at high Rey n o l d s n u m b e r s an d this extent o f p e n e t r a t i o n is

Steady l a m i n a r flow in a 90 degree p l a n a r b r a n c h 547

E q u a l e x i t p r e s s u r e s .

d s = 1 . 0 0

2 d ! - 0 . 7 5 u

a~ d m - 0 . 5 0

1 . . . . . . . . . ' . . . . . . . . . ' . . . . . . . . .

0 2 0 0 4 0 0 600

R e y n o l d s n u m b e r

Fig. 8. E x t e n t o f t h e recirculation z o n e in t h e side b r a n c h as a f u n c t i o n o f R e y n o l d s n u m b e r for v a r i o u s w i d t h s o f t h e side b r a n c h .

proportional to the width of the side branch. Since the wake must grow to its asymptotic size from an initial size of zero, it turns out to be a much more sensitive function of Reynolds number in the intermediate range of Re (50 ~ Re ~< 250) for side branches of larger size.

In the next series of simulations, the case of specified flow split boundary condition and its effect on the flow characteristics are considered. The specified flow split is always imposed at the exit of the main branch and only the pressure is specified at the exit of the side branch. Since an arbitrarily specified amount of fluid can be diverted into the side branch with this boundary specification, the flow in the side branch cannot be expected to reach a fully developed state in every case. In light of the arguments presented earlier, however on the parabolized nature of the problem at higher Reynolds numbers, reasonably accurate results are expected even if a recirculation zone crosses the exit boundary. The streamline contours are shown in Fig. 9 for various flow splits at a fixed Re = 300. At a flow split of 0.2 (Fig. 9a) 80% of the fluid flows through the side branch and hence the recirculation zone is carried far in the downstream direction within the side branch and its width is correspondingly decreased. A fairly large, but weak recirculation develops in the main branch also. As the flow split is increased, diverting more and more fluid into the main branch, the recirculation zone in the main branch becomes progressively weaker and finally disappears. The recirculation zone in the side branch becomes shorter in length and spans a wider cross section of the duct as the flow split is increased, approaching the behaviour observed earlier for equal exit pressure specification. In fact at one instance of flow split specification (flow split = 0.89 at Re = 300 from Fig. 6) the computed pressures at both the exits become equal.

The dimensionless pressure profiles along the left wall (x = 0, y), the centre line (x = 0.5, y) and the right wall (x = 1.0, y) of the main branch are shown in Figs 10 and 11 for Re = 300 and splits of 0.8 and 0.2 respectively. The dimensionless pressure is defined a s ( P i - P)/(½P(U) 2) where P i is the inlet pressure, P is the pressure at any y and ( U ) is the average inlet velocity. The solid lines are the pressures computed from the Poisson eqn (18) and the dashed lines are computed from the Uzawa algorithm. The dotted lines at the inlet and the outlet indicate the pressure gradient computed from the fully developed Poiseuille flow condition. The pressure rise near the side branch exit {ye (2, 3)} has been predicted in the finite difference computation of Liepsch et al. [9] and has been observed in the experimental measurements of Cho, Back and Crawford [24]. These figures also show that the magnitude of the pressure rise near the branch increases as the flow split is decreased (or the fraction flow in the side branch is increased). The pressure drop downstream of the branch is lower for lower flow splits which is to be expected as the flow rate in the main branch is smaller at lower flow splits. The experimentally observed pressure profiles in Ref. [24] for various flow splits also show the same trends as observed here. A quantitative comparison with Ref. [24] is not appropriate since the experiments were preformed in branches of circular cross section.

The effect of changing Reynolds number at a fixed flow split of 0.5 is shown in Fig. 12. The recirculation zones in both the branches are significantly affected with the increasing Reynolds number. In the first two cases (Figs 12a, b) a fully developed state is reached in both branches; but

548 R . E . HAVES et al.

9 ~ f (a)

6 I Split = 0 2 0

~ ~ o . ~ - ~ ~

0 2 3 4 5 6

I

(b)

Spiit = 0.50 •

0.45

2 3 4 5 6

(c)

Split = 0,80

0 I - L

0 1 2 3 4 5 6

Fig. 9. Effect o f increasing flow split on the flow field for the case o f specified flow split condition at Re = 300. (a) Split = 0.20, (b) split = 0.50, (c) split = 0.80.

this is clearly not the case for Re = 800. Since the fractional flow in each branch is equal, increasing Reynolds number stretches both the recirculation zones far in the downstream direction. The appearance of a third recirculation zone near the upper wall of the side branch and the incipient formation of a fourth vortex on the right wall of the main branch shown in Fig. 12(c), should not be surprising as they are the natural consequence of a strong inertial flow. In fact, in curved ducts, where the centrifugal force is sustained in the flow direction, spatially oscillating solutions in the flow direction have been computed recently [25]. In the present case since the centrifugal force is not sustained in the flow direction, such oscillations should decay in the flow direction. Karino et al. [26] have observed two vortices forming on opposite sides of the main branch in their experiments with a 90 ° Tee made of circular pipes (hence 3-D). Both these vortices form more readily at lower flow splits, corresponding to a larger fraction being diverted to the side branch thus increasing the magnitude of the centrifugal force. At a flow split of 0.1, the first one was observed at a Reynolds number of about 30 while the second one appearing downstream of the first one forms above a Reynolds number of 250. Because of the 2-D nature of the present computations and the 3-D nature of the experiments in Ref. [26] a direct quantitative comparison is not possible. There is however a qualitative agreement in the evolution of the flow structure with

Steady laminar flow in a 90 degree planar branch 549

0 . 2 . . . . i . . . . i . . . . i . . . . i • ' 00 t -O.S " ' . . . " , " "

° _ o . , , , , , , L . , , , , . . . : . . . . . . . 1 M . . . . i . . . . i . . . . J . . . . i • , 0,2

g) -0,2 ~ - - O ' s a ~ a a g g r t t h m ~l - - . . . . . . P o t s . e u ' t l l e J ' l .aw " , [

1 . 0 , , * , I * * , , I , i , ". I . . . . I • . I --0.4

. ~ . . . . i . . . . i . . . . i . . . . i • .

• ~ 0 . 0

- 1 . 0

2 4 6 O D i s t a n c e a l o n g m a i n b r a n c h

Fig. 10. Pressure variation along the main branch for Re = 300, flow split = 0.8

0

o~

0

0.4 . . . . i . . . . i . . . . I . . . . i ' '

0 . 2 = ~

0 . 0

- 0 . 2 0 . 4

. . . . ~ = . ~ 1 . . . . i . . . . i . . . . 1 ' ' 0.00'2

-0.2

-0.4

1 . 0

- 1 . 0

- 3 . 0

. . . . i . . . . i . . . .

, x=l.O / ~

i i . • i . . . . i . . . .

2 4

d l

. . , . t . , ]

8 S

D i s t a n c e a l o n g m a i n b r a n c h

Fig. 11. Pressure variation along the main branch for Re = 300, flow split = 0.2•

increasing Reynolds number and changing flow splits. In the 2-D, experimental work o f Liepsch e t al. [9] the velocity profiles in the side branch at such high Reynolds numbers are not given. Our computations with a longer duct (L = 26) did confirm the formation o f a second vortex in the side branch and the length in that case was sufficient for the flow to reach a fully developed state at the exit.

Figure 13 shows the vorticity distribution on the lower wall o f the side branch for various values o f flow splits and Reynolds numbers. In Fig. 13(a) the effect o f Reynolds number is shown at a fixed flow split o f 0.5. It is clear that the wake length increases significantly with increasing Reynolds number. Also s h o w n in Fig. 13(a) are the vorticity distributions obtained from a shorter duct (L = 3) for Re = 50 and 200. For Re = 50, the flow reaches a fully developed state within L = 3 and hence there is a perfect agreement between the two results. For the case o f Re = 200, a side branch length o f L = 3 is not adequate to reach a fully developed state. Hence the vorticity computed from the shorter duct ( 0 ) has not leveled off as in the case o f the profile computed from

550 R . E . HAYES e t al.

(a)

R e = 5 0

i d / ~ , = ; i l i =

: I!1111i,

0 1 2 3 4 5 6 7

i

4

3

2

[ (b) 8 ~j =

, 7

r R e : 2 0 0

e ;

0.45

i! (c)

R e = 8 ~

0 . 4 5 _ ~ ~,,.~03 ~ . ~ ~ - - - ~ O " q ' 8 " - - - ~ r " 0 . 0 3 " u . ~ . _ ~

F i g . 12. Effect o f increasing Reynolds number on flow field for the case o f specified flow split condition. (a) R e = 50, (b) R e = 200, (c) R e = 800. For the case o f R e = 800, a second recirculation zone near the upper plate o f the side branch exit is predicted with L = 6. A s i n F i g . 2 this has been confirmed with

another simulation with L = 26.

the longer duct ( - - ) . The agreement between the two profiles, including the prediction o f the reattachment point, is h o w e v e r quite remarkable over a significant portion o f the c o m m o n flow d omai n. The disagreement is confined to a few grid points from the exit o f the shorter duct. It is also clear that the flow has reached a fully developed state in the side branch even for Re = 400 with a duct length o f L = 6. In this case the length o f L = 6 is n o t adequate o n l y for Re = 600. Similar results are s h o w n in Fig. 13(b) at a fixed R e y n o l d s number o f 300 and various flow splits. In all these cases the flow has reached a fully developed state at the exit. A t a fixed R e y n o l d s number the wake length becomes a less sensitive f u n c t i o n o f the flow split o n c e the flow split is less than a b o u t 0.5. This becomes clearer in Fig. 14 where the reattachment p o i n t is s h o w n as a function o f R e y n o l d s number for three different flow splits. For flow splits o f 0.2 and 0.5 the reattachment p o i n t o f the wake varies in essentially the same manner as the R e y n o l d s number is increased. A l s o s h o w n in this figure are three data points for the location o f the reattachment p o i n t obtained from a short duct (L = 3) with a flow split o f 0.5. These provide yet another confirmation that the side branch length has very little influence o n the interior flow field.

Steady l a m i n a r flow in a 90 degree p l a n a r b r a n c h 551

0 >

( a )

0.10

0.08

0.08

0.04

0.02

0.00

- 0 . 0 2

- 0 . 0 4

- - R e = 6 0 ", R e - 2 0 0

- - - R e = 4 0 0 !

• .. Re=e00 i r e a t t a e h m e n t p o i n t s

- - ~ ' ~ . ~ . . . . / . "

, . , . , . . . '

2 3 4 5 8 7

D i s t a n c e a l o n g l o w e r wall (x2=2.0) of branch)

( b ) 0 . 0 6

Re=300.O, d s = 1.0

- - S p l i t - 0 . 2 0.04 - - - Split=0.5

. . Split=0.8

o.oo f/°°rn'r . . . . . . . t . . . . . . . " ~ N , / ~ " - , S r e a t t a h m

- 0 . 0 2 I \ ~ ~ p o i n t

- 0 . 0 4 I v , ' ' ' - - - - - ' ' ' . . . . . . . 1 3 5

Distance a l o n g l o w e r wall (x2=2.0) of branch)

Fig. 13. Vorticity d i s t r i b u t i o n o n t h e lower wall o f side b r a n c h . (a) Effect o f i n c r e a s i n g R e y n o l d s n u m b e r a t a fixed flow split o f 0.5, (b) effect o f i n c r e a s i n g flow split at a fixed Re o f 300.

The critical Reynolds number at which the recirculation zone begins to appear in the side branch was determined for a range of flow splits and side branch widths. This was accomplished by gradually increasing the Reynolds number in fixed increments at each flow split and branch width and for each run the sign change in the vorticity on the lower wall of the branch was monitored. When a sign change is noticed, indicating the existence of a recirculation zone, the search was refined by starting with the lower Reynolds number and repeating the search in smaller increments of Re. The critical Reynolds number was thus resolved to with in + 1. Figure 15 shows the variation

7 . . . . . . . . . . . . . . . . . . . , . . . . . . . . . , . . . . . . . . . . . . . . . . . . .

,~ 3 / / / I d s - 1.00 ,s / / / --- Split- o.~o 0o-,, duoO

/ / - - Split - 0.$0 (long d u c t ) • Split = 0.50 ( s h o r t d u e t )

- - - Split = 0.80 (long d u e t ) 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 2 0 0 400 600 800 1000 Reynolds n u m b e r

Fig. 14. E x t e n t o f t h e recirculation z o n e in t h e side b r a n c h a s a f u n c t i o n o f R e y n o l d s n u m b e r for t h e c a s e o f specified flow split.

552 R.E. HAYES el al.

I00

8o

80 s = 0 . 5 1

I~ 40 • " ~ , •

o s=0.7

~ 2o r..) "'-~. ds = 1.01

• Zarino et al. [ l l ] / 0 i i r i /

0.0 0 . 2 0 . 4 0 , 8 0 . 8 1 . 0

Flow split

Fig. 15. Variation of the critical Reynolds number with flow split for various side branch widths.

o f the critical R e y n o l d s n u m b e r with flow split a n d b r a n c h w i d t h . T h e critical R e y n o l d s n u m b e r d e c r e a s e s w i t h i n c r e a s i n g flow splits a n d it i n c r e a s e s w i t h d e c r e a s i n g b r a n c h w i d t h s . B o t h these t r e n d s are in q u a l i t a t i v e a g r e e m e n t with t h e m e a s u r e m e n t s o f K a r i n o e t al. [8, 26].

C O N C L U D I N G R E M A R K S

T h e l a m i n a r flow c h a r a c t e r i s t i c s in a 90 + p l a n a r b r a n c h h a v e b e e n i n v e s t i g a t e d u s i n g a finite e l e m e n t d i s c r e t i z a t i o n o f t h e g o v e r n i n g e q u a t i o n s o f m o t i o n . T h e c o m p u t e d velocity profiles a r e in e x c e l l e n t a g r e e m e n t w i t h o n e set o f e x p e r i m e n t a l l y m e a s u r e d v a l u e s a v a i l a b l e in t h e l i t e r a t u r e

[9]. T h e exit b o u n d a r y c o n d i t i o n is seen to h a v e o n l y m i n o r i n f l u e n c e o n the i n t e r i o r flow field w h i c h is r a t i o n a l i z e d as d u e to the p a r a b o l i c n a t u r e o f the flow a t h i g h R e y n o l d s n u m b e r s . T h e n u m b e r , l o c a t i o n , size a n d the s t r e n g t h o f the r e c i r c u l a t i o n z o n e s are s i g n i f i c a n t l y i n f l u e n c e d b y s u c h p a r a m e t e r s as flow split, b r a n c h w i d t h a n d R e y n o l d s n u m b e r . T h e critical R e y n o l d s n u m b e r a t w h i c h the first r e c i r c u l a t i o n z o n e is f o r m e d i n c r e a s e s w i t h d e c r e a s i n g flow splits a n d d e c r e a s i n g b r a n c h w i d t h s .

Acknowledgements--Financial support from the National Science and Engineering Research Council of Canada is gratefully acknowledged. A grant of computer time from the University of Alberta is also acknowledged.

R E F E R E N C E S

1. J. A. Fox and A. E. Hugh, Localization of Atheroma: a theory based on Boundary-Layer Separation. Br. Heart J. 28, 388-399 (1966).

2. N. S. Lynn, V. G. Fox and L. W. Ross, Computation of fluid-dynamical contributions to atherosclerosis at arterial bifurcations. Biorheology 9, 61-66 (1972).

3. L. W. Ehrlich and M. H. Friedman, Computational aspects of aortic bifurcation flows. Comput. Fluids 13, 177-183 (1985).

4. J. S. Bramley and D. M. Sloan, Numerical solution for two-dimensional flow in a branching channel using boundary-fitted coordinates. Comput. Fluids 15, 297-331 (1987).

5. V. O'Brien, L. W. Ehrlich and M. H. Friedman, Unsteady flow in a branch. J. Fluid Mech. 75, 315-336 (1976). 6. M. Siouffi, R. Pelissier, D. Farahifar and R. Rieu, The effect of unsteadiness on the flow through stenoses and

bifurcations. J. Biomechan. 17, 299-315 (1984). 7. T. Fukushima, T. Homma, T. Azuma and K. Harakaw, Characteristics of secondary flow in steady and pulsatile flows

through a symmetrical bifurcation. Biorheology 24, 3-12 (1987). 8. T. Karino, H. H. M. Kwong and H. L. Goldsmith, Particle flow in models of branching vessels: I. Vortices in 90 ~

T-junctions. Biorheology 16, 231-248 (1979). 9. D. Liepsch, S. Moravec, A. K. Rastogi and N. S. Vlachos, Measurements and calculation of laminar flow in a ninety

degree bifurcation. J. Biomechan. 15, 473-485 (1982). 10. D. K. Kreid, C.-J. Chung and C. T. Crowe, Measurements of the flow of water in a T junction by the LDV technique.

Trans. A S M E J. appl. Mech. 97, 498-499 (1975). I 1. M. Popp and D. W. Sallet, Experimental investigation of one- and two-phase flow through a tee junction. Int. Conf.

on the Physical Modelling o f Multi-Phase Flow, Paper B3, Coventry, England (1983). 12. J. A. C. Humphrey, A. M. K. Taylor and J. H. Whitelaw, Laminar flow in a square duct of strong curvature. J. Fluid

Mech. 83, 509-527 (1977). 13. W. H. Reid, On the stability of viscous flow in a curved channel. Proc. R. Soc. Lond. A244, 186-198 (1958).

Steady laminar flow in a 90 degree planar branch 553

14. A. Pollard, Computer modelling of flow in Tee-Junctions. Physico Chem. Hydrodyn. 2, 203-227 (1981). 15. R . T . Lahey Jr, Current understanding of phase-separation mechanisms in branching conduits. Nucl. Engng Design 95,

145-161 (1986). 16. H. Nasr-EI-Din, Slurry sampling from pipelines. In: Sewage and Wastewater Treatment, Vol. 3 (Edited by N. P.

Cheremisinoff). Gulf Publishing Company. In press. 17. A. Segev and K. C. Kern, Solid-liquid separation in a slurry manifold. A S M E Winter meeting, Synfuels and Coal Energy

Symposium, Dallas, U.S.A. (1985). 18. G. Bugliarello and G. C. Hsiao, Phase separation in suspensions flowing through bifurcations: a simplified

hemodynamic model. Science 141,469-471 (1964). 19. M. D. Gunzburger, R. A. Nicolaides and C. H. Liu, Algorithmic and theoretical results on computation o f

incompressible viscous flows by finite element methods. Comput. Fluids 13, 361-373 (1985). 20. M. Crouzeix and P. A. Raviart, Conforming and non-conforming finite element methods for solving the stationary

Stokes equations. R A I R O S~rie Analyse Num~rique 7, 33-76 (1973). 21. J. N. Reddy, An Introduction to the Finite Element Method. McGraw-Hill, New York (1984). 22. M. S. Engdlman, G. Strang and K. J. B~ithe, The application of quasi-Newton methods in fluid mechanics. Int. J. numer.

Meth. Engng 17, 707-718 (1981). 23. P. A. Tanguy, M. Fortin and L. Choplin, Finite element simulation o f dip coating. II. Non-Newtonian Fluids. Int.

J. num. Meth. 4, 459-475 (1984). 24. Y. I. Cho, L. H. Back and D. W. Crawford, Experimental investigation of branch flow ratio, angle and Reynolds

number effects on the pressure and flow fields in arterial branch models. Trans A S M E J. Biomech. Engng 107, 257-267 (1985).

25. S. Ravi Sankar, K. Nandakumar and J. H. Masliyah, Oscillatory flows in coiled square ducts. Phys. Fluids 31, 1348-1359 (1988).

26. T. Karino and H. L. Goldsmith, Particle flow in models of branching vessels: II. Effects of branching angle and diameter ratio on flow patterns. Biorheology 22, 87-104 (1985).

27. J. S. Bramley and S. C. R. Dennis, The numerical solution of two-dimensional flow in a branching channel. Comput. Fluids 12, 339-355 (1984).

28. E. M. Sparrow and R. G. Kemink, Heat transfer downstream of a fluid withdrawal branch in a tube. Trans. A S M E J. Heat Transfer 101, 23-28 (1979).

29. N. Ukeguchi, S. Sakata and T. Adachi, On the numerical analysis of compressible flow problems by the modified F L I C method. Comput. Fluids g, 251-263 (1980).

30. D. N. Ku and D. Liepsch, The effect of non-Newtonian viscoelasticity and wall elasticity on flow at a 90 ° bifurcation. Biorheology 23, 359-370 (1986).