PROCEEDINGS OF THE OKLAHOMA ACADEMY OF SCIENCE 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

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home

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.

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home

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:

c(¶T/¶t) = k[(1/X)(¶/¶X)(X¶T/¶X) + (¶2T/¶Y2)]

(1)

subject to the initial condition

and the boundary conditions in the clockwise direction

in Equations 1-8:

T = Temperature To = outer surface temperature TL = boiling cryogen temperature t = time X = radial distance Y = vertical distance = 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 p^th rectangular domain,

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

Define dimensionless time as

t = t (k/c)max Lmax2

(12)

where

and

(k/c)max = max of (k/c) 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}

where

j(T) = (k/c)T /(k/c)max

(16)

and the so-called aspect ratio

The initial condition is reduced to

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

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

and

in which

The convective boundary is

where

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

Equations 22a,b hold for : 0 £ x₃ £ 1, y₃ = 0.

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

and

in which

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

where

Equations 28a,b hold for x₂ = 1 and 0 £ y₂ £ 1.

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

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

so that

Substituting Eqs. 35 and 36 into Eq. 34 gives

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:

{Page 76}

where

j(Tp,ij) = (k/c)Tp,ij / (k/c)max

(39)

Initial Condition:

Boundary values:

where

where

but

where

in which

where

where

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,iN^y = 0, into Eq. 45b and solving Eqs. 45a,b simultaneously give

By substituting u_1,N^xj = 0 and arranging, Eq. 47a becomes

where

If k and h are functions of temperature, then substituting u_1,N^xj = 0 and comvining Eqs. 47b and 48b give

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

Substituting u_3,N^xj and u_3,1j = u_2,N^xj (Eq. 51a) into Eq. 51b and rearranging give

On the other hand Eq. 56 can be arranged as

Equations 60 and 61 can be solved simultaneously to obtain expressions for u_2,1j and u2,N^xj ( = u_3,1j); j = 2,3,. . . , (N^y - 1).

Substituting u_2,iN^y = 0 and rearranging Eq. 53a gives for constant l₂ value,

{Page 77}

where

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

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

NUMERICAL RESULTS

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home

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 L₁ = H₂ = 0.15 m, L₂ = H₁ = 0.30 m, h = 568 W/m² • 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 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.

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home

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}

REFERENCES

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home

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

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home

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 = 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

Introduction	Mathematical Formulation	Numerical Results	Conclusion	References	Symbolism and Nomenclature	Table of Contents	Home