C.M. Sliepcevich

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; 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.

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:

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/c)maxLmax2

(12)

where

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

(13)

and

(k/c)max = max of (k/c) inTL£T£To

(14)

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

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 = 1SNwij 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.

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.

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.

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.