PROCEEDINGS OF THE
OKLAHOMA ACADEMY OF SCIENCE
Library Digitization Homepage
OAS Homepage
Copyright
Search
Volume 65—1985

{Page 73}

Application of Differential Quadrature to Solution of Pool Boiling in Cavities

Faruk Civan

School of Petroleum and Geological Engineering, University of Oklahoma, 865 Asp Avenue, FH 236, Norman, OK 73019

C.M. Sliepcevich

Flame Dynamics Laboratory, University of Oklahoma, 1215 Westheimer Drive, Norman, OK 73069

The pool boiling problem in cavities is analyzed from the point of heat transfer by dividing the solution domain into three rectangular regions which are connected through the interface boundary conditions in the form of the compatibility and the continuity relationships. The numerical solution is accomplished by using the method of differential quadrature.

INTRODUCTION

Common engineering problems frequently deal with shapes of physical domains having irregular geometries as contrasted to regular ones which consist of single, uniform shapes, such as rectangular, cylindrical or spherical. In general, regular domains in which the physical properties are constant are amenable to analytical solutions. On the other hand, analytical solutions to irregular geometries are not available except in those cases where the irregular domain can be divided into regular shapes such that the system properties are constant within these regions, which are connected through the boundaries along which the system variables are established by the compatibility and the continuity relationships. Several recent examples of these applications are reports by : Hsieh and Su (8), who used a superposition technique along with separation of variables for a cavity geometry; Yi-Zhou and Yi-Heng (11), who used a harmonic function continuation technique on the torsion problems for bars with irregular cross-sections; and Heggs et al. (7), who analyzed a fin assembly heat transfer problem by separation of variables. All these examples were restricted to time-independent or steady-state cases.

However, for problems involving transient-state or steady-state cases with variable system properties, numerical methods such as finite differences, finite elements, or quadrature are unavoidable.

The purpose of this paper is to demonstrate a general approach for analyzing irregular geometries in either time-dependent or time-independent systems. This approach, which requires only a few discrete points to produce numerical solutions, is based on the method of quadrature that was originally proposed by Bellman et al. (1,2) as a technique for rapid solutions of one- and two-dimensional, initial value problems. Subsequently, Mingle (9,10) solved a one-dimensional initial and boundary value problem. Civan and Sliepcevich (3,4,5) extended the technique for multidimensional, steady-state and transient-state boundary value problems, the models for which may involve partial derivatives and integrals.

Essentially the method of differential quadrature replaces a partial space derivative by a linear weighted sum of the function values at the discrete points. As a result, a partial differential equation reduces to a set of algebraic equations if it is time-independent and a set of ordinary differential equations for a time-dependent case. Numerical solution methods for both cases are well developed.

In the present paper the application of the method of quadrature will be demonstrated only on a problem involving an irregular geometry.

MATHEMATICAL FORMULATION

For the purposes of the present study consider a transient-state conduction problem in

{Page 74}

an irregular, two-dimensional geometry for a medium whose properties are dependent only on its temperature. For example, consider the pool boiling of a cryogen in a cylindrical cavity drilled through a solid block, shown in Fig. 1, whose external surfaces are kept at a constant temperature. The heat transfer is represented by the following model:

symbolc(T/t) = k[(1/X)(/X)(XT/X) + (2T/Y2)] (1)

subject to the initial condition

T = To, t = 0   everywhere (2)

and the boundary conditions in the clockwise direction

T = To, Y = 0, L2 £ X £ (L1 + L2) (3)
T = To, X = L1 + L2, 0 £ Y £ (H1 + H2) (4)
T = To, Y = H1 + H2, 0 £ X £ (L1 + L2) (5)
- kT/X = 0, X = 0, H1 £ Y £ (H1 + H2) (6)
- kT/X = h(TL - T), Y = H1, 0 £ X £ L2 (7)
- kT/X = h(TL - T), X = L2, 0 £ Y £ H1 (8)

in Equations 1-8:

T = Temperature
To = outer surface temperature
TL = boiling cryogen temperature
t = time
X = radial distance
Y = vertical distance
symbol = density
c = specific heat capacity
k = thermal conductivity
h = film heat transfer coefficient
L1, L2, H1, H2 = characteristic lengths and heights for the block shown in Fig. 1

For computational convenience the geometrical domain of interest is separated into three rectangular domains (designated by the letter p) through the dashed lines as shown in Fig. 1. Next, each rectangular domain is reduced to a unit square as shown in Fig. 2 by defining the dimensionless distances in X- and Y-directions as, for the pth rectangular domain,

xp = (X - xpo)/Lp (9)
yp = (Y - ypo)/Hp (10)

in which X and Y are measured form the point 0(0,0) shown in Fig. 1, Lp and Hp are the width and height, and (xpo , ypo ) is the origin of the coordinates of the region p.

Define dimensionless temperature as

u = òT Tok(T')dT'/ òTLTok(T')dT'; 0 £ u £ 1 (11)

Define dimensionless time as

t = t (k/symbolc)max Lmax2 (12)

where

Lmax = maximum of (Lp; p = 1,2) (13)

and

(k/symbolc)max = max of (k/symbolc) in TL £ T £ To (14)

Thus, by using Equations 9 through 12, Equations 1 through 8 are normalized. The unsteady-state conduction equation becomes

{Page 75}

up /t = j(T)a p 2[(xp + xpo/ Lp)-1 up /xp + 2up /x p 2
      + b p 22up/y p 2] ; p = 1,2,3
(15)

where

j(T) = (k/symbolc)T /(k/symbolc)max (16)
ap = Lmax/Lp (17)

and the so-called aspect ratio

bp = (L / H)p. (18)

The initial condition is reduced to

u(xp,yp) = 0 , t = 0 ; p = 1,2,3 (19)

Boundary conditions for the unit square domain p = 1 in the clockwise direction:

u1 = 0,     0 £ x1 £ 1,     y1 = 0 (20)
u1 = 0,     x1 = 1,     0 £ y1 £ 1 (21)

The continuity and the compatibility conditions along the interface of domains 1 and 3 are, respectively,

u1 = u3 (22a)

and

(u/y)1 = g1(u/y)3,   0 £ x1 £ 1,   y1 = 1 (22b)

in which

g1 = H1 / H2. (23)

The convective boundary is

(u/x1) = l1(T1 - TL), x1 = 0, 0 £ y1 £ 1 (24)

where

l1 = L1h / òTLTok(T') dT' (25)

Boundary conditions for the unit square domain p = 3 in the clockwise direction:

Equations 22a,b hold for : 0 £ x3 £ 1, y3 = 0.

u3 = 0,     x3 = 1,     0 £ y3 £ 1 (26)
u3 = 0,     0 £ x3 £ 1,     y3 = 1. (27)

The continuity and the compatibility conditions between the domains 3 and 2 are

u3 = u2 (28a)

and

(u/x)2 = g2 (u/x)3,     x3 = 0,     0 £ y3 £ 1 (28b)

in which

g2 = L2/L1. (29)

Boundary conditions for the unit square domain p = 2 in the clockwise direction:

(u/y)2 = l2(T2 - TL),     0 £ x2 £ 1,     y2 = 0 (30)

where

l2 = H2h / òTLTok(T') dT'. (31)

Equations 28a,b hold for x2 = 1 and 0 £ y2 £ 1.

u2 = 0,     0 £ x2 £ 1,     y2 = 1 (32)
(u/x)2 = 0,     x2 = 0,     0 £ y2 £ 1 (33)

Solution of the Problem Using Differential Quadrature

The quadrature method replaces a linear operation, such as an integration or a differentiation, on a function by a linear weighted sum of the function values at discrete sample points (2). Therefore, it can be used to obtain numerical solutions. For example, consider

du/dx |x = xi   » j = 1SNwijuj (34)

in which xi are the sample points obtained by decomposing the independent variable range into N discrete points, ui are the function values at these points, and wij are the weighting coefficients attached to these points. To calculate the weighting coefficients consider a test function

uk = xk - 1 ; k = 1,2,. . . ,N (35)

so that

duk / dx = (k-1) xk-2 (36)

Substituting Eqs. 35 and 36 into Eq. 34 gives

(k-1) xk-2i =  j = 1S Nwij xk-2i;     i,k = 1,2,. . . ,N (37)

Equation 37 is solved for the weighting coefficients, wij, as described by Civan and Sliepcevich (3,4,5). These coefficients are then used in Eq. 34 to approximate the derivative at a discrete point in terms of the function values at all of the discrete points. Weighting coefficients for other types of linear operations can be obtained similarly.

In the following the weighting coefficients associated with the first and second order space derivatives with respect to x or y variables will be designated by superscripts in the form of a single and a double x or y, respectively.

When the partial space derivatives are replaced by formulae given by Civan and Sliepcevich elsewhere (5), Equations 15 through 33 reduce to the following equations, respectively.

Transient heat conduction equation:

(du/dt)p,ij » j(Tp,ij)a p 2{[xpi + (xpo / Lp)]-1

•   k = 1SNxw ik xup,kj +   k = 1SNxw ik xxup,kj + b p 2   k = 1SNyw jk yyup,ik}

p = 1,2,3; i = 2,3,. . . ,(Nx - 1); j = 2,3,. . . ,(Ny - 1)
(38)

{Page 76}

where

j(Tp,ij) = (k/symbolc)Tp,ij / (k/symbolc)max (39)
ap = Lmax / Lp (40)
bp = (L / H)p. (41)

Initial Condition:

up,ij = 0,     t = 0;     p = 1,2,3;    i = 1,2,. . . ,Nx,      j = 1,2,. . . ,Ny (42)

Boundary values:

u1,i1 = 0;     i = 1,2,. . . ,Nx (43)
u1,Nxj = 0;     j = 2,3,. . . ,(Ny - 1) (44)
u1,iNy = u3,i1 (45a)
  k = 1SNywNykyu1,ik = g1k = 1SNyw 1k yu3,ik (45b)

where

g1 = H1 / H2 (46)


      k = 1SNxw 1k yu1,kj = l1(u1,1j-1); j = 1,2,. . . ,Ny

if k and h are constant,
(47a)

where

l1 = L1h / k = a constant, (48a)

but

   k = 1SNxw 1k yu1,kj = l1,1j(T1,1j-TL);

   j = 1,2,. . . ,Ny   if k and h are functions of T,
(47b)

where

l1,1j = [L1 / òTLTok(T')dT' ] h (T1,1j) (48b)
u3,Nxj = 0; j = 1,2,. . . ,Ny (49)
u3,iNy = 0; i = 1,2,. . . ,Nx (50)
u2,Nxj = u3,1j (51a)


  k = 1SNxw Nxkyu2,kj = g2k = 1SNxw 1k yu3,kj ;     j = 2,3,. . . ,(Ny-1) (51b)

in which

g2 = L2 / L1 (52)


  k = 1SNxw 1k yu2,ik = l2(u2,i1 - 1); i = 1,2,. . . ,Nx     if k and h are constant (53a)

where

l2 = H2h / k = a constant (54a)


  k = 1SNxw 1k yu2,ik = l2,i1(T2,i1 - TL); i = 1,2,. . .Nx     if k and h are functions of T (53b)

where

l2,i1 = [H2 / òTLTok(T')dT' ] h (T2,i1) (54b)
u2,iNy = 0; i = 1,2,. . . ,Nx (55)


  k = 1SNxw 1k xu2,kj = 0;     j = 1,2,. . . , Ny (56)

Equations 38 through 56 given above can be employed together to obtain solution of the actual model at the prescribed discrete points as a function of time. However, before proceeding any further, we express the unknown boundary values, for which derivative boundary conditions are prescribed, and the continuity and the compatibility conditions in terms of the interior domain function values.

First we consider the continuity condition, Eqs. 45a,b. Substituting the known boundary values, u1,i1 = 0 and u3,iNy = 0, into Eq. 45b and solving Eqs. 45a,b simultaneously give

u1,iNy = u3,i1 =

[g1k = 2SNy-1w 1k yu3,ik -k = 2SNy-1w Nyky u1,ik] / [w NyNy y-g1w 11 y] ;

i = 2,3,. . . ,(Nx - 1
(57)

By substituting u1,Nxj = 0 and arranging, Eq. 47a becomes

u1,1j = (l1 + k = 2SNx-1w 1k xu1,kj) / (w 11 x- l1);

j = 2,3,. . . ,(Ny - 1)
(58a)

where

l1 = L1h / k = a constant (59)

If k and h are functions of temperature, then substituting u1,Nxj = 0 and comvining Eqs. 47b and 48b give

w 11 xu1,1j- [L1 / òTLTok(T')dT' ] h (T1,1j)(T1,1j-TL)

= -k = 2SNy-1w 1k yu1,kj ;

j = 2,3,. . .,(Ny - 1)
(58b)

This equation has to be solved for u1,1j simultaneously with Eq. 11, which relates T to u, by a trial-and-error method.

Substituting u3,Nxj and u3,1j = u2,Nxj (Eq. 51a) into Eq. 51b and rearranging give

  xw Nx1u2,1j + ( xwNxNx - g2 xw11)u2,Nxj

=g2k = 2SNx-1w 1k xu3,kj -k = 2SNx-1w Nxkx u2,kj
(60)

On the other hand Eq. 56 can be arranged as

  xw11u2,1j + xw1,Nxu2,Nxj= -k = 2SNx-1w 1k xu2,kj (61)

Equations 60 and 61 can be solved simultaneously to obtain expressions for u2,1j and u2,Nxj ( = u3,1j); j = 2,3,. . . , (Ny - 1).

Substituting u2,iNy = 0 and rearranging Eq. 53a gives for constant l2 value,

u2,i1 = ( l2 + k = 2SNy-1w 1k yu2,ik) / (w 11 y - l2) ;

i = 2,3,. . . ,Nx
(62)

{Page 77}

where

l2 = H2h / k = a constant (63)

When k and h are temperature dependent, a procedure similar to that used in deriving Eq. 58b can be employed to obtain

w 11 yu2i- [H2 / òTLTok(T')dT' ] h (T2,i1)(T2,i1-TL)

= k = 2SNy-1w 1k yu2,ik ; i = 2,3,. . .,Nx
(64)

which can be solved for u2i1 by means of Eq. 11 via a trial-and-error method.

NUMERICAL RESULTS

Numerical solution of Eq. 38 has been accomplished using a low-order Runge-Kutta-Fehlberg four (five) method as described by Fehlberg (6) until the steady state is reached. For the boundaries at which the temperatures are prescribed, these values were used directly. For the boundaries at which the heat fluxes (or temperature derivatives) are prescribed, the temperatures are calculated in terms of the unknown interior grid point values using Eqs. 57 through 64, which are also used to calculate the values for the temperatures along these boundaries following the numerical solution of Eq. 38 for the interior grid points. Only the steady state numerical solution is presented in Table 1 for the case of L1 = H2 = 0.15 m, L2 = H1 = 0.30 m, h = 568 W/m2 • K, k = 8.6 W/m • K, TL = -161 °C and To = 0 °C using N = 6 discrete points in the x and y-directions of each of the three domains.

This steady state numerical solution has been obtained directly by solving Eq. 38 after dropping the time derivative term and solving the resulting set of algebraic equations simultaneously.

CONCLUSION

The method of quadrature can be applied to problems involving irregular geometry.

ACKNOWLEDGMENTS

The authors acknowledge the support of University Technologists, Inc., the Graduate College, the School of Petroleum and Geological Engineering, and the Merrick Computing Center of the University of Oklahoma.

{Page 78}

REFERENCES

1.   R.E. Bellman, B.G. Kashef, and J. Casti, J. Comput. Phys. 10, 40-52 (1972).

2.   R.E. Bellman, Methods of Nonlinear Analysis, Vol. II, Chap. 16, Academic Press, New York, 1973.

3.   F. Civan and C.M. Sliepcevich, Internat. J. Numer. Methods Engr. 19, 711-724 (1983).

4.   F. Civan and C.M. Sliepcevich, J. Math. Anal. Appl., 93, 206-221 (1983).

5.   F. Civan and C.M. Sliepcevich, J. Math. Anal. Appl., 101, 423-443 (1984).

6.   E. Fehlberg, NASA TR R 315, July 1969.

7.   P.J. Heggs, D.B. Ingham, and M. Manzour, ASME J. Heat Transfer 104, 210-212 (1982).

8.   C.K. Hsieh and K.C. Su, ASME J. Heat Transfer 103, 42-46 (1981).

9.   J.O. Mingle, Internat. J. Numer. Methods Engr. 7, 103 - 116 (1973).

10.   J.O. Mingle, J. Math. Anal. Appl. 60, 559-569 (1977).

11.   C. Yi-Zhou and C. Yi-Heng, Int. J. Engr. Sci. 19, 791-804 (1981).

SYMBOLISM AND NOMENCLATURE

c = specific heat capacity
d = ordinary differentiation operator
h = film heat transfer coefficient
H = characteristic height
k = thermal conductivity
L = characteristic length
t = time
T = temperature
u = dimensionless temperature
x = dimensionless radial distance
X = radial distance
y = dimensionless vertical distance
Y = vertical distance
w = weighting coefficients
= partial differentiation operator
ò = integration operator
Greek Symbols
a = Lmax/ Lp
b = aspect ratio
g = a constant
l = constant
symbol = density
t = a dimensionless time
j = diffusivity ratio
S = summation convention
Subscripts
i,j = indices indicating the location of a discrete point
L = boiling point
max = maximum
N = number of discrete points
o = outer surface
p = region index
Superscripts
x = associated with the first-order partial derivatives in the x-direction
xx = associated with the second-order partial derivatives in the x-direction
y = associated with the first-order partial derivatives in the y-direction
yy = associated with the second-order partial derivatives in the y-direction