Bem_arbitrary_duct_2007_babel.doc

BOUNDARY ELEMENT SOLUTION OF LAMINAR VISCOUS FLOW IN ARBITRARY DUCT Dr. Karima E. Amori Mech. Eng. Dep., University of Baghdad, Iraq

ABSTRACT A boundary element scheme for 2-D steady incompressible laminar viscous flow inside arbitrary duct is presented. The solution provides computation of velocity distribution, average velocity, and local wall shear stress. The computational results obtained are verified by comparing them with the corresponding exact solution for circular and elliptical ducts, which shows an excellent agreement. Moreover the generality of this scheme is illustrated by applying it to ducts of hydraulic diameter (Dh=1.0) and different cross sections namely square, rectangular, triangular and circular tube at different orientation angle of its axis from vertical axis. The obtained results show that the local shear stress is distorted for ducts other than circular duct while the mean shear stress is constant for the same pressure gradient and hydraulic diameter. ‫الخلصاة‬ ‫طبععاقي ولععزج فععي مجععرى ذو مقطععع عععام‬, ‫ غيععر انضععغاطي‬,‫ مسععتقر‬,‫يتضمن البحث حل"ا لجريععان ثنععائي البعععد‬ .‫ ومعدل السرعة وتوزيع أجهععاد القعص الموضعععي‬,‫ يتم تحديد توزيع السرعة‬.‫بأستخدام طريقة العناصر المحيطية‬ ‫تم مقارنة النتائج المحسوبة المستحص لة معع النتائعج النظريعة لجريعان داخ ل مجعرى دائعري و بيضعوي المقطعع‬ ‫ لتوضيح الستخدام العام لهذه الطريقة المطروحععة تععم تطبيقهععا علععى مجععارى‬.‫والتي أظهرت تقارب ممتاز بينهما‬ ‫ انبوب دائري ذو زعانف داخليععة و انبععوب دائععري‬, ‫ مثلث‬, ‫ مستطيل‬, ‫ وذو مقطع مربع‬1= ‫ذي قطر هيدروليكي‬ ‫ أوضحت النتائععج أن أجهععاد القععص الموضعععي يتغيععر علععى طععول محيععط‬.‫متغير زوايا الميل مع المحور الشاقولي‬ .‫ اما متوسط أجهاد القص فل يتغير لنفس القطر الهيدروليكي و الضغط المسلط‬,‫المجرى عدا المجرى الدائري‬ Key wards: Viscous, Flow, Boundary Element method, Duct INTRODUCTION The problem of viscous fluid flow inside arbitrary ducts presents itself in many engineering applications such as compact heat exchangers, gas turbine cooling systems, cooling channels in combustion chambers, nuclear reactors and others. Analytical solution of this problem is possible if the duct has circular or elliptical cross section Bachlor, (1967), equilateral triangular cross section Tao, (1961). An experimental investigation by Carlson and Irvine (1961) has been done for isosceles triangular, a noncircular duct was studied by Barrow and Hassan (1984). Hussain,I.y. (1997) used a finite element method to solve the problem in equilateral triangular duct. The present work represents an approach to use the boundary element method to solve developed laminar viscous flow in ducts of symmetric and asymmetric cross section such as circular, elliptical, square, rectangular, triangular (equilateral, isosceles and asymmetric). The flow inside internal finned tubes, and a pipe at different orientation angle (0 to 90o) has been also investigated. MATHEMATICAL MODEL The shape of an arbitrary duct and the coordinate system used during this work is shown in Fig.(1).

1

C

y

duct boundary

z

Flow

o'

x

n

Fig.(1): Coordinate System for Fluid Flow in Arbitrary Duct

Governing equation of fully developed steady unidirectional laminar viscous incompressible flow is: Khader(1981)  2u z 

s 

(1)

where u z = fluid velocity (m/s) s

 dp  g z cos  dz

( the negative pressure gradient incorporating the effect

  fluid viscosity (N.s/ m2)

of gravity along the z- axis) ( N/m3)

  angle between duct axis and y-axis

. The solution of eq. (1), is decomposed into homogenous solution (f) and a particular solution( fp) (i.e) (2) uz  f  f p these solutions must satisfy the following conditions  2f  0

(3)  2f

let

S  S fp  v( x , y )

p



(4)



(5) where (v) is a plane function taken as 1 v  [( x  x o  ) 2  ( y  y o  ) 2 ] 4

(6)

where (xo’ ,yo’) are the coordinates of arbitrary reference point in x,y plane, such that it satisfies the following 2v  1

(7) also v  . H (8)

2

where H  h x i  h y j 

x3 y3 i j 12 12

(9)

Hence S

fp 



.H

(10) For no slip boundary condition at the solid boundary of the duct ( u z  0 ) hence f  f p

(11)

Using the boundary integral formulation for the homogenous solution (f) at a point (xo) inside the tube an integral equation for normal derivative (∂f / ∂n) can be written as:

 G ( x, x o )

c

f ( x ) s s dl( x )   v( x o )  v( x ) ( n ( x ) . G ( x , x o ) n   c



dl( x )

(12) A derivation of this equation is given in the appendix Flow Rate Computation Flow rate is defined as Q

 u z dA A



 (f

 f p ) dA 

A

 f dA   f A

p

dA

A

(13) After some mathematical manipulation for eq.(13) leads to Q   n . ( vu z  c

s



H ) dl( x )

(14) A brief derivation to this equation is given in the appendix Equation (14) is a boundary integral equation which can be solved after the normal gradient of (f) at the boundary (c ) is computed from eq.(12).The average velocity is calculated as Vav 

Q A

(15)

Computation of Shear Stress The shear stress at the duct wall is defined as u  wall   z   n . u z n substituting eq.(2) in eq.(16) to get  wall   n.(f  f

p

)

or  wall   (n . f )   ( n . f p ) (17) substituting eq.(6) in eq.(17) to get the local wall shear stress as:

3

(16)

 wall  

f  s  s  [nx ( v)  n y ( v)] n x  y 

or  wall  

f s  [ n x (x  x o )  n y ( y  y o  ) n 2

(18)

The mean wall shear stress can be evaluated as  mean 

1  wall dl( x ) L c

(19) where L is the perimeter of the duct Computation of Velocity Distribution Inside the Domain The total velocity value uz at a point (xo) inside the tube can be evaluated once the  f 

when value of  n  at each element is computed. Referring to eq. (2)  i p uz (xo) = f(xo) +f (xo)

(20)

where N f s   f (x o )       G ( x , x o ) dl( x )   v( x ) n ( x ) . G ( x , x o ) c i  1 n  i c i

(21) and f p (x o ) 



s s 1  v( x o )  *   x o  x o   2   y o  y o   2     4 



dl( x )

(22)

COMPUTATIONAL PROCEDURE The main flow parameters namely, velocity distribution, flow rate and local wall shear stress can be evaluated according to the following procedure: 1-The duct boundary (C) is divided into (N) boundary elements e i (of constant type) defined by (N+1) nodes.  f 

eq.(12) by considering  n  as a constant unknown value over  i element ei , hence eq(12) can be rewritten as: 2- Discretizing

s s N  f  G ( x , x ) dl ( x )   v ( x )      vi n( x) . G( x, xo ) dl ( x) o o    i 1 ei i 1  n  i ei (23) N

this equation forms a system of (N+1) linear equations which is solved by Gauss elimination method to evaluate local  f n  i . 3- Compute the velocity distribution inside the domain from eq. (20) after substituting the values of ( f and f p ) from eq.s (21) and (22). 4-Compute flow rate and average velocity by using the boundary integral eq.s(14) and eq.(15) respectively. 4

5-Compute the elemental and mean shear stress value using eq.s (18 and 19). COMPUTER PROGRAM A computer program was developed in Fortran-90 to perform the numerical solution formulated previously. The program consists of five parts. The first part for input data (namely hydraulic diameter=1.0, s/μ=1.0+9.81cos θ), while the second to the generation of element matrices, and the others related to the assymbler, solver and output data presentation respectively. The execution time of the developed program on a personal computer (Pentium IV ) of (256 k RAM) took about less than (1 second). DISCUSSION AND RESULTS The forgoing computational approach is applied to problems of known exact solutions as a rigorous test. Fig.2 shows a comparison between the obtained numerical results of radial velocity distribution in circular duct of radius (R=0.5unit or hydraulic diameter=1.0unit ) with the corresponding exact solution as given in Bachlor(1967) uz 

sR 4

r 2  1  ( R ) 

(24).

for different number of boundary elements. It is clear that the accuracy of solution become very well when the boundary elements be (N=16). Fig.3 shows that the maximum computational velocity value (which is equal to twice the average value) exist at the duct center which agrees with the analytical solution . Fig.4 shows the velocity distribution inside an elliptical duct of dimensions (a=1.0 unit and b=0.36 unit). The results compared with the analytical solution given by Bachlor(1967) uz 

s 2  (a

2

b

2

x 2 y 2  1  ( a )  ( b )  )

(25)

an excellent agreement between them is obtained. Also it is clear that the maximum value of velocity ratio (V/Vaverage) is traced at duct center ( the same value as that for circular duct). Fig.5 shows that the variation of predicted local wall shear stress at each boundary element number for circular and elliptical ducts. It is clear that the local value of wall shear stress is constant along the perimeter of circular duct, since every boundary point has the same distance away from the center. A comparison between the computational values and the corresponding analytical values (derived the appendix) for elliptical duct shows an excellent agreement, also it can be noticed that the values of local (τwall) be minimum at its ends of the major axis and maximum at the tips of minor axis, this is due to the difference in duct curvatures at these positions. Fig.6 shows the velocity distribution inside rectangular ducts for different aspect ratios. It predicts that the velocity values decreases as the aspect ratio increases for the same hydraulic diameter (D h =1.0 unit). The local shear stress along the longer side of this duct increases as the longer side length increases while there is no variation noticed along the shorter sides as shown in Fig.7.This is due to the weaker velocity variation obtained as the side length increases for the same flow rate. The local shear stress for ducts of triangular section was affected by the apex angle as shown in Fig.8. The figure shows that the local (τwall) increases along the side walls of the duct as the apex angle decreases, while there is no effect on the base side. Fig.9 shows the velocity variation along y-axis of triangular duct for different apex angle. It can be deduced that the maximum velocity traced at the duct center and it increases as the apex angle decreases. The above trained has been also shown for unsymmetrical 5

triangular duct as shown in Fig.10. A cross section of internal finned circular tube (6 fins) is shown in Fig.11-a. The local shear stress variation along a pie section boundary of this duct (divided into 18 segments as shown in Fig.11-b is shown in Fig.11-c. It can be noticed that the local (τwall) value vanishes at the fluid plain of symmetry and has a maximum value at the fin tips. This is due to the tube geometry distortion. The velocity variation along selected lines within the pie section is presented in Fig.12, in all cases the velocity was found to be maximum at the duct center and vanishes at the duct contour. The radial velocity distribution for a circular tube at different inclination angles with y-axis varies from (0 o-90o) is shown in Fig.13. It is clear that the velocity increases as the inclination angle increases for the same supply pressure(s).This is due to the effect of gravity force. CONCLUSION The boundary element approach derived in this work for fully developed laminar flow in ducts has been successfully applied. It is clear that the numerical solution was capable in predicting the complete velocity distribution and wall shear stress (local and mean) in arbitrary duct section, with minimum effort compared with FDM and FEM. Another primary advantage of the present method is the space reduction of the problem from 2-D to one dimension. References: Bachlor,G.K. “ An Introduction to Fluid Dynamics” Cambridge University Press. Cambridge. 1967. Barrow,H. ; Hassan,A.K.A. and Avegerinos, C. “Peripheral Temperature Variation in the Wall of a Noncircular Duct-An Experimental Investigation” Int. J. Heat Mass Transfer, vol. 27,No.7, p.p.1031-1037,1984 Brebbia, C.A. and Walkrt, S. “Boundary Element Techniques in Engineering”, McGraw-Hill,1980 Carlson,L.W.; Irvine,JR.T.F. “Fully Developed Pressure Drop in Triangular Shaped Ducts” J. Heat Transfer Nov. 1961 p.p. 441-444 Hussain,I.y. “A finite element algorithm for Laminar Forced Convection in Equilateral Triangular Duct” Proceedings of The Fourth Scientific Engineering Conference, College of Engineering, University of Baghdad,(18-20)Nov.1997 Khader,M.S. ”A Surface Integral Numerical Solution for Laminar Developed Duct Flow” J. Applied Mechanics, Transaction of ASME Dec.1981,vol.48,p.p.695-700 “Fully Developed Flow and Heat Transfer in Ducts Having Streamwise-Periodic Variations of Cross Sectional Area” J. Heat Transfer, May,Vol.99,p.p. 180-186, 1977 Partap,V.S. and Spalding D.B. “ Fluid Flow and Heat Transfer in ThreeDimensional Duct Flows” Int. J. Heat Mass Transfer, vol.19, p.p.1183-1188 Tao,L.N. ”On Some Laminar Forced-Convection Problems”Trans. ASME, J. of Heat Transfer Nov.1961, p.p.166

6

A n a ly t ic a l

0 .0 6

N o . o f e le m e n ts = 4 N o . o f e le m e n ts = 8 N o . o f e le m e n ts = 1 6

V e lo c ity (m /s )

0 .0 5

0 .0 4

0 .0 3

0 .0 2

0 .0 1

0 .0 0 0

2

4

6

8

10

12

14

16

B o u n d a ry N o d e N o .

Fig.(2): Effect of Number of Boundary Elements on Numerical Results of Circular Duct 0 .0 8 2 .0

0 .0 6

0 .0 4

1 .0

0 .0 2

0 .5 c o m p . v e lo c it y A n a ly t. v e lo c it y V / V a v e ra g e

0 .0 0

0 .0 0 .0 0

0 .2 0

0 .4 0

r /R

0 .6 0

0 .8 0

1 .0 0

Fig.(3): Radial Velocity Distribution in Circular Duct

7

V /V a v e r a g e

v e lo c ity (m /s )

1 .5

2 .0 0 .0 6

1 .5

0 .0 3

V /V a v e r a g e

1 .0

V

v e lo c ity (m /s )

0 .0 5

0 .5

0 .0 2 C o m p . v e lo c it y A n a l. v e lo c it y V / V a v e ra g e

0 .0 0

0 .0 0 .0 0

0 .2 0

0 .4 0

0 .6 0

x /a o r y /b

0 .8 0

1 .0 0

Fig.( 4 ) : Velocity Distribution in Elliptical Duct 0 .3 5 C o m p u t . ( E llip .) A n a ly . (E llip . ) C o m p u t . ( C ir c u la r )

S h e a r S tr e s s (N /m ^ 2 )

0 .3 0

0 .2 5

0 .2 0

0 .1 5

0 .1 0 0

10

20

30

40

B o u n d a ry N o d e N o .

50

60

70

Fig.(5): Local Shear Stress Distribution Along Duct Boundary (Circular and Elliptical)

8

0 .3 0

2A

2B

0 .2 5

v e lo c ity (m /s )

0 .2 0

0 .1 5

0 .1 0

0 .0 5 A =B A = 1 .5 B A =2

B

0 .0 0 -1 .0

-0 .5

0 .0

0 .5

1 .0

x /A

Fig.(6 ) : Velocity Distribution in Rectangular Duct

(a) 1 .2 0 A=B

1 .1 0

A = 1 .5 B A=2B

1 .0 0 0 .9 0

64 32

S h e a r S tr e s s (N /m ^ 2 )

y

x

0 .8 0 0 .7 0 0 .6 0 0 .5 0 0 .4 0 0 .3 0 0 .2 0

96

0 .1 0

1

0 .0 0 -0 .1 0

128

0

20

40

60

(b) Fig.(7): a) Boundary Node No., b) Local Shear Stress Distribution Along Rectangular Duct Boundary for Different Aspect Ratio

9

80

Node No.

100

120

140

1 196

0 .4 5 0 .4 0 0 .3 5

S h e a r S tr e s s (N /m ^ 2 )

0 .3 0 0 .2 5 0 .2 0 0 .1 5 0 .1 0 0 .0 5

y

a = 60 a = 30

0 .0 0

128

a = 20

64

a = 15

-0 .0 5

x

0

20

40

60

80

100

120

B o u n d a ry N o d e N o .

(a) (b) Fig.(8): a) Boundary Node No. ; b) Local Shear Stress Distribution Along Triangular Duct Boundary For The Same Hydraulic Diameter and Different Apex Angle 2 .5 a =60 a =30 a =20

2 .0

a =15

V

V / V a v e ra g e

y

b

V

1 .5

x

2b

a

1 .0

0 .5

0 .0 -2 .0

-1 .5

-1 .0

-0 .5

y /b

10

0 .0

0 .5

1 .0

140

160

180

200

Fig.( 9 ): Velocity Distribution Along y-Axis in Triangular Duct for Different Apex Angle 0 .6

0 .5

V

y

V / V a v e ra g e

V

0 .4 x

0 .3

0 .2

0 .1

0 .0 -0 .5

- 0 .4

-0 .3

-0 .2

- 0 .1

0 .0

y

0 .1

0 .2

0 .3

Fig.(10): Velocity Distribution Along y-Axis in Unsymmetrical Triangular Duct

12 13 11

Li

ne

A

14 10

16 15

L in

17

e B 9

18

60

1

(a) (b)

11

30

45

2

L in e C

5 3

4

6

7

8

1 .0

S h e a r S tr e s s (N /m ^ 2 )

0 .8

0 .6

0 .4

0 .2

0 .0 0

2

4

6

8

10

Node No.

12

14

16

18

(c) Fig.(11): a) Internal Finned Duct ; b) Pie Section of Duct ; c)Local Shear Stress Distribution Along Boundary of a Pie Section of Finned Duct 0 .0 3 0

L in e A

0 .0 2 5

L in e B L in e C

v e lo c ity (m /s )

0 .0 2 0

0 .0 1 5

0 .0 1 0

0 .0 0 5

0 .0 0 0 0 .0 0

0 .2 0

0 .4 0

0 .6 0

0 .8 0

1 .0 0

r /R o

Fig.(12): Velocity Distribution Along Selected Lines(A,B and C) in Finned Duct

12

0 .8

Q= 0 Q= 30 Q= 60

0 .7

Q= 90

v e lo c ity (m /s )

0 .6

0 .5

0 .4

0 .3

0 .2

0 .1

0 .0 0 .0

0 .2

0 .4

0 .6

0 .8

1 .0

r / R

Fig.(13): Radial Velocity Distribution in Circular Duct Located at Different Orientation Angle

APPENDIX Boundary Integral Eq. (12) f ( x o )    G (x, x o ) c

f ( x ) dl( x )   f ( x )( n ( x ) . G ( x , x o ) dl( x ) n c

(26) where G is the 2-D Green’s function defined as: G





1 1 1 ln X  X o  ln  x  x o  2   y  y o  2 2 2 2

(27) where G is a differentiable function and satisfies the singularity forced Laplaces equation:  2 G ( X, X o )   ( X  X o )  0

(28) where δ is the Dirac delta function Substituting eq. (11) in (26) when the point (xo) approaches the tube boundary (c) and using eq. (5) to get f ( x o )    G ( x, x o ) c

f ( x ) s dl( x )   v( x ) (n ( x ) . G ( x , x o ) n c

(29) 13



dl( x )

Using Green’s second identity which states that G  2 f  f  2 G   . ( G f  f G )

(30) eq. (30) can be rewritten as: G  2 f  f  ( X  X o )   . ( G f  f G )

(31) substituting eq. (3) to get f  ( X  X o )   . (G f  f G )

(32) Now integrating both sides over the domain area A



f  (X  X o ) dA 

A



  . (G f  f G ) dA

A

(33) then applying the divergence theorem to get f ( x o )    n . (Gf )dl  c

(34) using n . f 



n . (f G )dl

c

f ,Brebbia, and Walkrt(1980), to get : n f (x o )   G c

f dl   f n . G dl n c

(35) If the field point (xo) approaches the duct boundary then f ( x o )  f p ( x o )

substituting eq.(5) to get

s f v( x o , y o )    G dl   n c

 f n . G dl c

Flow Rate Eq.(16) Dealing with the first term on RHS of eq.(13) it can be rewritten as

 f dA   (f  A

2

v  v  2 f ) dA

A

(36) then using the divergence theorem

 (f 

2

v  v  2 f ) dA 

A

c n . ( v f

(37) at the boundary ( c ) using eq. (11) v f  f v  v  f  (  f p ) v

 v f  (   v ( f 

s v ) v 

s v) 

14

 f v)dl( x )

using eq. (5) to get v f  f v  n . v  ( f  f p )

  f dA  A

c n . v u z dl( x )

(38) Now using the divergence theorem on the second term on RHS of eq.(13) such as:

 f

p

dA 

A

s

   . H A

dA   n .  c

s



H dl( x )

(39) then substituting equations (38) and (39) in eq.(13) to get Q   n . ( vu z  c

s



H ) dl( x )

Wall Shear Stress of Elliptical Duct u z u u  n . u z  z n x  z n y n x y

(40) Now referring to equation (25)  wall  

u z s  2x  2y  ( nx  ny) n 2(a  2  b  2 ) a 2 b2

(41) where (x,y)coordinates of a point belongs to the duct contour

15