
{Page 73}
School of Petroleum and Geological Engineering, University of Oklahoma, 865 Asp Avenue, FH 236, Norman, OK 73019
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.
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; YiZhou and YiHeng (11), who used a harmonic function continuation technique on the torsion problems for bars with irregular crosssections; and Heggs et al. (7), who analyzed a fin assembly heat transfer problem by separation of variables. All these examples were restricted to timeindependent or steadystate cases.
However, for problems involving transientstate or steadystate 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 timedependent or timeindependent 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 twodimensional, initial value problems. Subsequently, Mingle (9,10) solved a onedimensional initial and boundary value problem. Civan and Sliepcevich (3,4,5) extended the technique for multidimensional, steadystate and transientstate 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 timeindependent and a set of ordinary differential equations for a timedependent 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.
For the purposes of the present study consider a transientstate conduction problem in
{Page 74}
an irregular, twodimensional 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:
c(¶T/¶t) = k[(1/X)(¶/¶X)(X¶T/¶X) + (¶^{2}T/¶Y^{2})]  (1) 
subject to the initial condition
T = T_{o}, t = 0 everywhere  (2) 
and the boundary conditions in the clockwise direction
T = T_{o}, Y = 0, L_{2} £ X £ (L_{1} + L_{2})  (3) 
T = T_{o}, X = L_{1} + L_{2}, 0 £ Y £ (H_{1} + H_{2})  (4) 
T = T_{o}, Y = H_{1} + H_{2}, 0 £ X £ (L_{1} + L_{2})  (5) 
 k¶T/¶X = 0, X = 0, H_{1} £ Y £ (H_{1} + H_{2})  (6) 
 k¶T/¶X = h(T_{L}  T), Y = H_{1}, 0 £ X £ L_{2}  (7) 
 k¶T/¶X = h(T_{L}  T), X = L_{2}, 0 £ Y £ H_{1}  (8) 
in Equations 18:
T = Temperature T_{o} = outer surface temperature T_{L} = boiling cryogen temperature t = time X = radial distance Y = vertical distance = density c = specific heat capacity k = thermal conductivity h = film heat transfer coefficient L_{1}, L_{2}, H_{1}, H_{2} = 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 Ydirections as, for the p^{th} rectangular domain,
x_{p} = (X  x_{po})/L_{p}  (9) 
y_{p} = (Y  y_{po})/H_{p}  (10) 
in which X and Y are measured form the point 0(0,0) shown in Fig. 1, L_{p} and H_{p} are the width and height, and (x_{po} , y_{po} ) is the origin of the coordinates of the region p.
Define dimensionless temperature as
u = òT T_{o}k(T')dT'/ òT_{L}T_{o}k(T')dT'; 0 £ u £ 1  (11) 
Define dimensionless time as
t = t (k/c)_{max} L_{max}2  (12) 
where
L_{max} = maximum of (L_{p}; p = 1,2)  (13) 
and
(k/c)_{max} = max of (k/c) in T_{L} £ T £ T_{o}  (14) 
Thus, by using Equations 9 through 12, Equations 1 through 8 are normalized. The unsteadystate conduction equation becomes
{Page 75}
¶u_{p} /¶t = j(T)a p 2[(x_{p} + x_{po}/ L_{p})^{1} ¶u_{p} /¶x_{p} + ¶^{2}u_{p} /¶x p 2 + b p 2¶^{2}u_{p}/¶y p 2] ; p = 1,2,3 
(15) 
where
j(T) = (k/c)_{T} /(k/c)_{max}  (16) 
a_{p} = L_{max}/L_{p}  (17) 
and the socalled aspect ratio
b_{p} = (L / H)_{p}.  (18) 
The initial condition is reduced to
u(x_{p},y_{p}) = 0 , t = 0 ; p = 1,2,3  (19) 
Boundary conditions for the unit square domain p = 1 in the clockwise direction:
u_{1} = 0, 0 £ x_{1} £ 1, y_{1} = 0  (20) 
u_{1} = 0, x_{1} = 1, 0 £ y_{1} £ 1  (21) 
The continuity and the compatibility conditions along the interface of domains 1 and 3 are, respectively,
u_{1} = u_{3}  (22a) 
and
(¶u/¶y)_{1} = g_{1}(¶u/¶y)_{3}, 0 £ x_{1} £ 1, y_{1} = 1  (22b) 
in which
g_{1} = H_{1} / H_{2}.  (23) 
The convective boundary is
(¶u/¶x_{1}) = l_{1}(T_{1}  T_{L}), x_{1} = 0, 0 £ y_{1} £ 1  (24) 
where
l_{1} = L_{1}h / òT_{L}T_{o}k(T') dT'  (25) 
Boundary conditions for the unit square domain p = 3 in the clockwise direction:
Equations 22a,b hold for : 0 £ x_{3} £ 1, y_{3} = 0.
u_{3} = 0, x_{3} = 1, 0 £ y_{3} £ 1  (26) 
u_{3} = 0, 0 £ x_{3} £ 1, y_{3} = 1.  (27) 
The continuity and the compatibility conditions between the domains 3 and 2 are
u_{3} = u_{2}  (28a) 
and
(¶u/¶x)_{2} = g_{2} (¶u/¶x)_{3}, x_{3} = 0, 0 £ y_{3} £ 1  (28b) 
in which
g_{2} = L_{2}/L_{1}.  (29) 
Boundary conditions for the unit square domain p = 2 in the clockwise direction:
(¶u/¶y)_{2} = l_{2}(T_{2}  T_{L}), 0 £ x_{2} £ 1, y_{2} = 0  (30) 
where
l_{2} = H_{2}h / òT_{L}T_{o}k(T') dT'.  (31) 
Equations 28a,b hold for x_{2} = 1 and 0 £ y_{2} £ 1.
u_{2} = 0, 0 £ x_{2} £ 1, y_{2} = 1  (32) 
(¶u/¶x)_{2} = 0, x_{2} = 0, 0 £ y_{2} £ 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 = 1SNw_{ij}u_{j}  (34) 
in which x_{i} are the sample points obtained by decomposing the independent variable range into N discrete points, u_{i} are the function values at these points, and w_{ij} are the weighting coefficients attached to these points. To calculate the weighting coefficients consider a test function
u_{k} = x^{k  1} ; k = 1,2,. . . ,N  (35) 
so that
du_{k} / dx = (k1) x^{k2}  (36) 
Substituting Eqs. 35 and 36 into Eq. 34 gives
(k1) xk2i = j = 1S Nw_{ij} xk2i; i,k = 1,2,. . . ,N  (37) 
Equation 37 is solved for the weighting coefficients, w_{ij}, 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(T_{p,ij})a p 2{[x_{pi} + (x_{po} / L_{p})]^{1} • • k = 1SN^{x}w ik xu_{p,kj} + k = 1SN^{x}w ik xxu_{p,kj} + b p 2 k = 1SN^{y}w jk yyu_{p,ik}} p = 1,2,3; i = 2,3,. . . ,(N^{x}  1); j = 2,3,. . . ,(N^{y}  1) 
(38) 
{Page 76}
where
j(T_{p,ij}) = (k/c)T_{p,ij} / (k/c)_{max}  (39) 
a_{p} = L_{max} / L_{p}  (40) 
b_{p} = (L / H)_{p}.  (41) 
Initial Condition:
u_{p,ij} = 0, t = 0; p = 1,2,3; i = 1,2,. . . ,N^{x}, j = 1,2,. . . ,N^{y}  (42) 
Boundary values:
u_{1,i1} = 0; i = 1,2,. . . ,N^{x}  (43) 
u_{1,}_{Nxj} = 0; j = 2,3,. . . ,(N^{y}  1)  (44) 
u_{1,iNy} = u_{3,i1}  (45a) 
k = 1SN^{y}w_{Nyk}yu_{1,ik} = g_{1}k = 1SN^{y}w_{ 1k} yu_{3,ik}  (45b) 
where
g_{1} = H_{1} / H_{2}  (46) 
k = 1SN^{x}w_{ 1k} yu_{1,kj} = l_{1}(u_{1,1j}1); j = 1,2,. . . ,N^{y} if k and h are constant, 
(47a) 
where
l_{1} = L_{1}h / k = a constant,  (48a) 
but
k = 1SN^{x}w_{ 1k} yu_{1,kj} = l_{1,1j}(T_{1,1j}T_{L}); j = 1,2,. . . ,N^{y} if k and h are functions of T, 
(47b) 
where
l_{1,1j} = [L_{1} / òT_{L}T_{o}k(T')dT' ] h (T_{1,1j})  (48b) 
u_{3,Nxj} = 0; j = 1,2,. . . ,N^{y}  (49) 
u_{3,iNy} = 0; i = 1,2,. . . ,N^{x}  (50) 
u_{2,Nxj} = u_{3,1j}  (51a) 
k = 1SN^{x}w_{ Nxk}yu_{2,kj} = g_{2}k = 1SN^{x}w_{ 1k} yu_{3,kj} ; j = 2,3,. . . ,(N^{y}1)  (51b) 
in which
g_{2} = L_{2} / L_{1}  (52) 
k = 1SN^{x}w_{ 1k} yu_{2,ik} = l_{2}(u_{2,i1}  1); i = 1,2,. . . ,N^{x} if k and h are constant  (53a) 
where
l_{2} = H_{2}h / k = a constant  (54a) 
k = 1SN^{x}w_{ 1k} yu_{2,ik} = l_{2,i1}(T_{2,i1}  T_{L}); i = 1,2,. . .N^{x} if k and h are functions of T  (53b) 
where
l_{2,i1} = [H_{2} / òT_{L}T_{o}k(T')dT' ] h (T_{2,i1})  (54b) 
u_{2,iNy} = 0; i = 1,2,. . . ,N^{x}  (55) 
k = 1SN^{x}w_{ 1k} xu_{2,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, u_{1,i1} = 0 and u_{3,iNy} = 0, into Eq. 45b and solving Eqs. 45a,b simultaneously give
u_{1,iNy} = u_{3,i1} = [g_{1}k = 2SN^{y}1w_{ 1k} yu_{3,ik} k = 2SN^{y}1w_{ Nyk}y u_{1,ik}] / [w_{ NyNy} yg_{1}w_{ 11} y] ; i = 2,3,. . . ,(N^{x}  1 
(57) 
By substituting u_{1,Nxj} = 0 and arranging, Eq. 47a becomes
u_{1,1j} = (l_{1} + k = 2SN^{x}1w_{ 1k} xu_{1,kj}) / (w_{ 11} x l_{1}); j = 2,3,. . . ,(N^{y}  1) 
(58a) 
where
l_{1} = L_{1}h / k = a constant  (59) 
If k and h are functions of temperature, then substituting u_{1,Nxj} = 0 and comvining Eqs. 47b and 48b give
w_{ 11} xu_{1,1j} [L_{1} / òT_{L}T_{o}k(T')dT' ] h (T_{1,1j})(T_{1,1j}T_{L}) = k = 2SN^{y}1w_{ 1k} yu_{1,kj} ; j = 2,3,. . .,(N^{y}  1) 
(58b) 
This equation has to be solved for u_{1,1j} simultaneously with Eq. 11, which relates T to u, by a trialanderror method.
Substituting u_{3,Nxj} and u_{3,1j} = u_{2,Nxj} (Eq. 51a) into Eq. 51b and rearranging give
xw_{ Nx1}u_{2,1j} + ( xw_{NxNx}  g_{2} xw_{11})u_{2,Nxj} =g_{2}k = 2SN^{x}1w_{ 1k} xu_{3,kj} k = 2SN^{x}1w_{ Nxk}x u_{2,kj} 
(60) 
On the other hand Eq. 56 can be arranged as
xw_{11}u_{2,1j} + xw_{1,Nx}u_{2,Nxj}= k = 2SN^{x}1w_{ 1k} xu_{2,kj}  (61) 
Equations 60 and 61 can be solved simultaneously to obtain expressions for u_{2,1j} and u2,N^{x}j ( = u_{3,1j}); j = 2,3,. . . , (N^{y}  1).
Substituting u_{2,iNy} = 0 and rearranging Eq. 53a gives for constant l_{2} value,
u_{2,i1} = ( l_{2} + k = 2SN^{y}1w_{ 1k} yu_{2,ik}) / (w_{ 11} y  l_{2}) ; i = 2,3,. . . ,N^{x} 
(62) 
{Page 77}
where
l_{2} = H_{2}h / 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} yu_{2i} [H_{2} / òT_{L}T_{o}k(T')dT' ] h (T_{2,i1})(T_{2,i1}T_{L}) = k = 2SN^{y}1w_{ 1k} yu_{2,ik} ; i = 2,3,. . .,N^{x} 
(64) 
which can be solved for u_{2i1} by means of Eq. 11 via a trialanderror method.
Numerical solution of Eq. 38 has been accomplished using a loworder RungeKuttaFehlberg 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 L_{1} = H_{2} = 0.15 m, L_{2} = H_{1} = 0.30 m, h = 568 W/m^{2} • K, k = 8.6 W/m • K, T_{L} = 161 °C and T_{o} = 0 °C using N = 6 discrete points in the x and ydirections 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.
The method of quadrature can be applied to problems involving irregular geometry.
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}
1. R.E. Bellman, B.G. Kashef, and J. Casti, J. Comput. Phys. 10, 4052 (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, 711724 (1983).
4. F. Civan and C.M. Sliepcevich, J. Math. Anal. Appl., 93, 206221 (1983).
5. F. Civan and C.M. Sliepcevich, J. Math. Anal. Appl., 101, 423443 (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, 210212 (1982).
8. C.K. Hsieh and K.C. Su, ASME J. Heat Transfer 103, 4246 (1981).
9. J.O. Mingle, Internat. J. Numer. Methods Engr. 7, 103  116 (1973).
10. J.O. Mingle, J. Math. Anal. Appl. 60, 559569 (1977).
11. C. YiZhou and C. YiHeng, Int. J. Engr. Sci. 19, 791804 (1981).
