Heat Transfer with *solidification and freezing* is of importance in many instances, for example, in engineering ice making, in freezing of foods, in icing of electric power plants, in frost and ice formation of aircraft, during the casting of the metal, and in many geophysical problems such as frost heaving (see also Frosting).

Recently, solidification has been spotlighted from the viewpoint of latent heat thermal energy storage (LHTES) like ice storage. This type of *energy storage* system contributes to reduced carbon dioxide emission, since the coefficient of performance is considerably improved by employing the LHTES concept. Another application is manufacturing of semiconductor waifer of silicon.

Mathematically, solidification is associated with the so-called *moving boundary problem* (MBP) which is characterized by the possession of moving interfaces dividing the relevant field into two regions. Such a problem becomes perfectly nonlinear because the position of the moving interfaces are neither fixed in space nor known *a priori.* Due to this nonlinearity, the analytical solution can be found only in limited situations, for example, in *Neumann's classical solution* for a one-dimensional problem.

In this entry, only recent advances in experimental, analytical, and numerical aspects in solidification and freezing problem will be introduced.

The reader is referred to the recent books and reviews; e.g. *Freezing and Melting Heat Transfer in Engineering* (eds. K. C. Cheng and N. Seki), Hemisphere, 1991; Solidification of Pure Liquids and Liquid mixtures Inside Ducts and over External Bodies, by S. Fukusako and M. Yamada, *Applied Mech.,* Review, 47-12 (1994).

An application of the *coordinate transformation method* to the multidimensional MBP was proposed by Saitoh. This method can be regarded as an extension of the Landau transformation method to the multidimensional case, which was presented qualitatively in his paper of 1950 for the one-dimensional problem. The most remarkable feature of the multidimensional MBP is that the moving interface exists in an arbitrarily shaped domain. Therefore, the general numerical method should be one which can consider the arbitrariness of the shapes of both the moving front and the domain boundary.

The boundary fixing method (BFM) considers the arbitrary geometry of both the moving interface and the domain boundary via a change of an independent variable. BFM was also studied by Duda et al., almost concurrently with the author's study.

Computed results for an ellipse with a short diameter of 0.08 m, an aspect ratio of 3/2, and a cooling rate of 3.333 × 10^{−3} K/s is shown in Figure 1. The cooling rate C_{1} is defined as

Here, T_{w}(t) is the cooling surface temperature in t time. Next, as an example of the case in which natural convection in the liquid is combined with heat conduction in the solid, the two-dimensional freezing around a horizontal circular cylinder in quiescent water was solved by the BFM. For the analysis, it was assumed that:

Flow is laminar.

Properties are constant.

*Boussinesq approximation*is valid.

For brevity, only computed results will be shown. Figure 2 shows a comparison between numerical and experimental freezing front contours at different times. The water temperature is 6.9°C, the diameter of the cylinder is 0.04 m, and the cylinder surface temperature is –12.4°C. Experimental data obtained by Saitoh and Hirose are shown by the circle. The agreement between the two is excellent except near the vicinity of the downstream edge, where separation may occur in the real flow. This was not considered in the analysis.

As shown in the above, the boundary fixing method is a very simple tool for solving multi-dimensional MBP's. For this reason, this method has been widely used to date.

The enthalpy method (previously called the Conventional method) was first developed for one-dimensional MBP's by Eyres, et al. and Baxter. However, a thorough formulation was presented by Katayama and Hattori (1974), and Shamsunder and Sparrow (1975) at almost the same time. In this method, the latent heat is treated as a part of the thermophysical properties: namely, as heat content.

According to Shamsunder and Sparrow the energy conservation equation which holds true in the domain of the phase change is expressed by the following (see Figure 3):

The validity of the above equation is easily confirmed by the fact that the usual energy equation is reducible directly from the above equation.

Katayama and Hattori have introduced an apparent heat content method employing the concept of the Dilac delta function for the freezing and melting problem with a discrete temperature. Temperature function is defined in Figure 4.

For example, the finite difference formulation for the one-dimensional problem is represented as

Here, the relationship between c, ρ*, and L is given by

It is assumed that the latent heat is liberated in the small temperature range 2ε. The temperature function is the determined so that the total heat content coincides with the latent heat L.

A numerical example is shown in Figure 5 in which computed isotherms are shown when the porous region around two cooled cylinders starts freezing.

Although two well-used methods, i.e., the Boundary fixing method and the Enthalpy method, have many advantages over other existing methods, these methods still have severe disadvantages. Namely, BFM largely depends on the smoothness of the boundary shape function, and with the Enthalpy method it is impossible to include a natural convection effect in the liquid phase.

An alternative method was recently developed by Saito and Kato (1986). This method is based on the multilateral element (MEM) method, in which the representative point is placed on the barycenter of a triangular or quadrilateral element. Note that it is not the circumcenter but the barycenter where the representative point is positioned. The existing triangular element method (TEM) has used the circumcenter. The barycenter method makes it possible to consider the arbitrary geometry of both the fixed and the moving boundaries, thereby eliminating difficulties encountered with, for example, BFM. Another advantage of this method is that the usual explicit finite difference can be incorporated in the formulation. This considerably reduces the troublesome numerical procedures that are necessary prior to computation. New freezing front positions are calculated by using temperature distribution in the vicinity of the freezing front. Then, the adjacent element is enlarged by the same amount as the newly frozen area thus obtained. This procedure is continued in a step-by-step manner until the size of the frozen area reaches nearly the same element size as that of the old element.

Since the freezing front produced in this method closely resembles the growth rings of trees, it is called the "Growth Ring Method (GRM)."

As an example, the numerical results for the two-dimensional freezing of water in a regular prism initially at its fusion temperature, has been taken up as a standard 2-D freezing problem.

The computed results by the GRM are shown in Figure 6 in which freezing front contours and generated elements (quadrilateral element in this case) are simultaneously designated. The side length of the regular square, Stefan number Ste, and the time step used in the computation, were 0.079 m, 0.49, and 0.003, respectively. The initial freezing front position in the normal direction was set to be 0.92. Therefore, every circumferential contour beyond this initial line denotes the freezing front and the quadrilateral element the mesh generated.

The computer running time (CPU) for this case was some 100s on NEAC 2200/ACOS 1000 (30 MFLOPS) at Tohoku University. Note that the use of the SOR scheme instead of the usual explicit scheme is particularly desirable since the mesh size in this example gets smaller as time elapses.

As seen from the above illustrative examples, the Growth Ring Method is a quite simple, but very powerful method applicable to arbitrary multidimensional domains. The GRM can also handle problems involving natural convection in the liquid phase.

It is advantageous that the mesh may be regenerated halfway through the computation in order to reduce the number of elements and resultant computer time.

#### REFERENCES

Saitoh, T. (1978) Numerical method for multi-dimensional freezing problems in arbitrary domains, *J. Heat Transfer*, Trans. ASME. 100, 294.

Shamsunder, N. and Sparrow, E. M. (1975) Analysis of multidimensional conduction phase change via the enthalpy model, *J. Heat Transfer,* Trans. ASME, 97, 333.

Katayama, K. and Hattori, M. (1974) A study of heat conduction with Freezing (1st Report, Numerical Method of Stefan Problem), *Trans. JSME*, 40-33, 1404.

Saitoh, T. and Kato, K. (1986) Numerical analysis of multidimensional freezing problems via growth ring method, *Proc. 23rd Heal Transfer Symp. Japan*, 695.