Back to Transport Phenomena in Multiphase Systems Home
Table of Contents
6 MELTING AND SOLIDIFICATION 6.1 Introduction Melting and solidification find application in the geophysical sciences; industrial forming operations such as casting and laser drilling; latent heat energy storage systems; and food and pharmaceutical processing. Any manmade metal products must undergo liquid forms at some point during manufacturing processes and solidify to form intermediate or final products. Melting and solidification processes can be classified as one of three types: one-region, two-region, and multiple-region. The classification depends on the properties of the phase change material (PCM) involved and the initial conditions. For a single-component PCM, melting or solidification occurs at a single temperature. Pure water, for example, melts at a uniform temperature of 0 °C, while pure n-Octadecane (C18H38) melts at 28 °C. For the solid-liquid phase change process of a PCM with a single melting point, the solid-liquid interface appears as a clearly-observable sharp border. Initial conditions for the solid-liquid phase change process of a single-component PCM determine whether the problem will be classified as a one- or two-region problem. For the melting (or solidification) process, if the initial temperature of the PCM, Ti, equals the melting point, Tm, the temperature in the solid (liquid) phase remains uniformly equal to the melting point throughout the process. In this case, only the temperature distribution in the liquid (solid) phase needs to be determined. Thus, the temperature of only one phase is unknown and the problem is called a one-region problem. Fig. 6.1 shows the temperature distribution of one-dimensional, one-region melting and solidification problems. The surface temperature, T0, is greater than Tm for melting and is below Tm for solidification. For the melting process, if the initial temperature of the PCM, Ti, is below the melting point of the PCM, Tm, (or above, for solidification), the temperature distribution of both the liquid and solid phases must be determined; this is called a two-region problem. Fig. 6.2 shows the temperature distribution of onedimensional two-region melting and solidification problems. Chapter 6 Melting and Solidification 421 T T0 T Ti, Tm Ti, Tm T0 s(t) (a) Melting ( T0 > Tm ) Figure 6.1 One-region melting and solidification. x s(t) (b) Solidification ( T0 < Tm ) x T T0 Tm Ti T Ti Tm T0 s(t) (a) Melting ( T0 > Tm > Ti ) Figure 6.2 Two-region melting and solidification. x s(t) (b) Solidification ( T0 < Tm < Ti ) x For a multi-component PCM, the solid-liquid phase change process occurs over a range of temperatures (Tm1, Tm2), instead of a single temperature. The PCM is liquid if its temperature is above Tm2 and solid when its temperature is below Tm1. Between the solid and liquid phases there is a mushy zone where the temperature falls within the phase change temperature range (Tm1, Tm2). Successful solution of phase change problems involving these substances requires determination of the temperature distribution in the liquid, solid, and mushy zones; therefore, these are referred to as multiregion problems. The temperature distribution of one-dimensional solidification in a multicomponent PCM is shown in Fig. 6.3, where it can be seen that the solution requires tracking of two interfaces. In solid-liquid phase change problems, the location of the solid-liquid interface is unknown before the final solution is obtained and this presents a 422 Transport Phenomena in Multiphase Systems T Ti Tm2 Tm1 x Figure 6.3 Solidification of a multicomponent PCM. Th Ti Th Ti Figure 6.4 Conduction- and convection-controlled melting. special difficulty. Since the interface also moves during melting or solidification, such problems are referred to as moving boundary problems and always have time as an independent variable. Natural convection induced by the temperature gradient in the liquid and mushy zones can have a significant effect on melting and solidification processes. The importance of natural convection in solid-liquid phase change problems was not recognized until the late 1970s. Fig. 6.4 shows the effect of natural convection on a melting process in a rectangular cavity. The left wall temperature is at Th ( Th > Tm ) and right wall temperature is at the initial temperature, Ti ( Ti ≤ Tm ). When conduction controls the melting process, the melting front moves parallel to the heated vertical wall, while natural convection in the liquid causes the front to accelerate in the upper portion of the cavity; this occurs because the warmer liquid near the heated wall rises while the colder fluid Chapter 6 Melting and Solidification 423 near the interface descends. As a result, the shape of the melting front is inclined, and the rate of melting increases due to the improved heat transfer in the liquid. Convection has very strong influence on the growth microstructure during solidification of pure metal or alloys. This chapter presents a number of prototypical problems of melting and solidification; its goal is to establish the physics of this form of phase change and to demonstrate the variety of tools available for its solution. Section 6.2 discusses boundary conditions at the solid-liquid interface, including continuity of temperature and conservation of energy at the interface. Section 6.3 presents the exact solutions of one-dimensional solid-liquid phase problems, including both one-region and two-region problems. Integral approximate solutions of various solid-liquid phase change problems are presented in Section 6.4. Integral approximate solution is applied to one-region and two-region problems under different boundary conditions and coordinate systems, as well as to solidification of binary solutions. Since the known exact analytical solutions cover a limited number of situations, numerical simulation methods are also addressed in Section 6.5 through three examples: the enthalpy method, the equivalent heat capacity method, and the temperature-transforming model. Analytical and numerical solutions for solidification of a binary solution are presented in Section 6.6. Analysis of contact melting in a cavity – a topic of importance in latent heat thermal energy storage applications – is covered in Section 6.7. Section 6.8 discusses melting and solidification in porous media, followed by a detailed discussion on application of solid-liquid phase change. This chapter closes with a discussion on microscale phase change (Section 6.9), which is essential for understanding advanced materials processing methods. 6.2 Boundary Conditions at the Solid-Liquid Interface For a solid-liquid phase change of a PCM with a single melting temperature, the solid-liquid interface clearly delineates the liquid and solid phases. The boundary conditions at this interface must be specified in order to solve the problem. The boundary conditions to be specified in this section are special cases of those discussed in Chapter 3. As shown for the one-dimensional melting problem illustrated in Fig. 6.2, the solid-liquid interface separates the liquid and solid phases. According to eq. (3.195), the temperatures of the liquid and solid phases near the interface must equal the temperature of the interface, which is at the melting point, T m . Therefore, the boundary conditions at the interface can be expressed as TA ( x, t ) = Ts ( x, t ) = Tm , x = s (t ) (6.1) where TA ( x, t ) and Ts ( x, t ) are the temperatures of the liquid and solid phases, respectively. The density of the PCM usually differs between the liquid and solid phases; therefore, density change always accompanies the phase change process. The 424 Transport Phenomena in Multiphase Systems solid PCM is usually denser than the liquid PCM, except for a few substances such as water and gallium. For example, the volume of paraffin, a very useful PCM for energy storage systems, expands about 10% when it melts. Therefore, the density of the liquid paraffin, ρA , is less than the density of the solid paraffin, ρ s . When water freezes, however, its volume increases, so the density of ice is less than that of liquid water. The density change that occurs during solid-liquid phase change will produce an extra increment of motion in the solid-liquid interface. For the melting problem in Fig. 6.2(a), where the liquid phase velocity at x = 0 is zero, if the density of the solid is larger than the density of the liquid (i.e., ρ s > ρA ) the resulting extra motion of the interface is along the positive direction of the x axis. Assume that the velocity of the solid-liquid interface due to phase change is u p = ds / dt , while the extra velocity of the solid-liquid interface due to density change is uρ . The density change must satisfy the conservation of mass at the interface, i.e., ρ s (u p − u ρ ) = ρ A u p (6.2) which follows from eq. (3.166). The extra velocity induced by the density change can be obtained by rearranging eq. (6.2) as ρ − ρA uρ = s up (6.3) ρs which is also valid for case ρ s > ρA , except the extra velocity becomes negative. Another necessary boundary condition is the energy balance at the solidliquid interface, which is a special case of eq. (3.180). If the enthalpy of the liquid and solid phases at the melting point are hA and hs , the energy balance at the solid-liquid interface can be expressed as: ′′ ′′ qA − qs = ρ A u p hA − ρ s (u p − uρ )hs x = s (t ) (6.4) ′′ ′′ where qA and qs are the heat fluxes in the x-direction in the liquid and solid phases, respectively. Substituting eq. (6.3) into eq. (6.4), one obtains ds ′′ ′′ (6.5) qA − qs = ρA hsA x = s (t ) dt where hsA = hA − hs is the latent heat of melting. If convection in the liquid phase can be neglected and heat conduction is the only heat transfer mechanism in both the liquid and solid phases, the heat flux in both phases can be determined by Fourier’s law of conduction: ∂T ( x, t ) ′′ (6.6) qA = −kA A ∂x x = s ( t ) ′′ qs = − k s ∂Ts ( x, t ) ∂x (6.7) x = s (t ) Chapter 6 Melting and Solidification 425 The energy balance at the solid-liquid interface for a melting problem can be obtained by substituting eqs. (6.6) and (6.7) into eq. (6.5), i.e., ∂T ( x, t ) ∂T ( x, t ) ds (t ) ks s − kA A = ρ A hsA x = s (t ) (6.8) ∂x ∂x dt For a solidification process, the energy balance equation at the interface can be obtained by a similar procedure: ∂T ( x, t ) ∂T ( x, t ) ds (t ) (6.9) ks s − kA A = ρ s hsA x = s (t ) ∂x ∂x dt The only difference between eqs. (6.8) or (6.9) is the density on the righthand side of the equation. If the temperature distributions in the liquid and solid phases are known, the location of the solid-liquid interface can be obtained by solving eqs. (6.8) or (6.9). It should be noted that density change causes advection in the liquid phase, which further complicates the problem. For the melting problem, if the heat transfer mechanism in the liquid phase is convection, the heat flux in the liquid phase can be obtained by Newton’s law of cooling: ′′ qA = h(Tw − Tm ) (6.10) where h and T∞ in eq. (6.10) are the convective heat transfer coefficient and the bulk temperature of the liquid phase, respectively. The energy balance at the solid-liquid interface is ∂T ( x, t ) ds (t ) ks s + h(T∞ − Tm ) = ρ A hsA x = s (t ) (6.11) ∂x dt For the solidification process, the heat flux in the liquid phase is ′′ qA = h(Tm − T∞ ) (6.12) The energy balance equation at the solid-liquid interface becomes ∂T ( x, t ) ds (t ) − h(T∞ − Tm ) = ρ s hsA (6.13) ks s x = s (t ) ∂x dt For the natural-convection controlled melting in a rectangular cavity, the solid-liquid interface is inclined along the height of the rectangular cavity as shown in Fig. 6.4(b). For solidification in a rectangular cavity, on the other hand, the interface will also be inclined, because natural convection accelerates solidification at the lower portion but decelerates solidification in the upper portion. Assuming n is a unit vector along the normal direction of the solidliquid interface, the boundary conditions at the interface can be expressed as TA ( x, y, t ) = Ts ( x, y, t ) = Tm x = s ( y, t ) (6.14) ∂T ( x, y, t ) ∂T ( x, y, t ) (6.15) ks s − kA A = ρ A hsA vn x = s ( y, t ) ∂n ∂n where vn is the solid-liquid interface velocity along the n-direction. It is apparent that eq. (6.15) is not convenient for numerical solution because it contains temperature derivatives along the n-direction. Suppose the shape of the solid-liquid interface can be expressed as x = s ( y, t ) (6.16) 426 Transport Phenomena in Multiphase Systems Equation (6.16) can then become the following form (see Problem 6.6; Ozisik, 1993): ª § ∂s · 2 º ª ∂T ∂T º ∂s x = s ( y, t ) (6.17) «1 + ¨ ¸ » « k s s − kA A » = ρA hsA ∂x ¼ ∂t « © ∂y ¹ » ¬ ∂x ¬ ¼ Similarly, for a three-dimensional melting problem with an interface described by z = s ( x, y , t ) (6.18) the energy balance at the interface is ª § ∂s · 2 § ∂s · 2 º ª ∂T ∂T º ∂s z = s ( x, y, t ) (6.19) «1 + ¨ ¸ + ¨ ¸ » « k s s − kA A » = ρA hsA ∂x ¼ ∂t « © ∂x ¹ © ∂y ¹ » ¬ ∂x ¬ ¼ For solidification problems, it is necessary to replace the liquid density ρA in eqs. (6.17) and (6.19) with the solid-phase density ρ s . The density change in solid-liquid phase change is often neglected in the literature in order to eliminate the additional interface motion discussed earlier. 6.3 Exact Solution 6.3.1 Governing Equations of the Solidification Problem The physical model of the solidification problem to be investigated in this subsection is shown in Fig. 6.2 (b), where a liquid PCM with a uniform initial temperature Ti , which exceeds the melting point Tm , is in a half-space x > 0. At time t = 0, the temperature at the boundary x = 0 is suddenly decreased to a temperature T0 , which is below the melting point of the PCM. Solidification occurs from the time t = 0. This is a two-region solidification problem as the temperatures of both the liquid and solid phases are unknown and must be determined. It is assumed that the densities of the PCM for both phases are the same. Natural convection in the liquid phase is neglected, and therefore the heat transfer mechanism in both phases is pure conduction. The temperature in the solid phase must satisfy ∂ 2T1 1 ∂T1 0 < x < s (t ), t > 0 (6.20) = ∂x 2 α1 ∂t T1 ( x, t ) = T0 x = 0, t > 0 (6.21) For the liquid phase, the governing equations are ∂ 2T2 1 ∂T2 s (t ) < x < ∞, t > 0 (6.22) = 2 α 2 ∂t ∂x T2 ( x, t ) → Ti x → ∞, t > 0 (6.23) Chapter 6 Melting and Solidification 427 T2 ( x, t ) = Ti x > 0, t = 0 (6.24) The boundary conditions at the interface are T1 ( x, t ) = T2 ( x, t ) = Tm x = s (t ), t > 0 (6.25) ∂T ∂T ds k1 1 − k2 2 = ρ hsA x = s (t ), t > 0 (6.26) dt ∂x ∂x Before obtaining the solution of the above problem, a scale analysis of the energy balance eq. (6.26) is performed. At the solid-liquid interface, the scales of derivatives of solid and liquid temperature are ∂T1 Tm − T0  (6.27) s ∂x ∂T2 Ti − Tm  (6.28) ∂x s The scale of the interface velocity is ds s (6.29)  dt t Substituting eqs. (6.27) – (6.29) into eq. (6.26), one obtains T −T T −T s k1 m 0 − k2 i m  ρ hsA (6.30) s s t Rearranging eq. (6.30), the scale of the location of the solid-liquid interface is obtained as follows: § k T −T · s2 (6.31)  Ste ¨1 − 2 i m ¸ α1t © k1 Tm − T0 ¹ where (6.32) hsA is the Stefan number. Named after J. Stefan, a pioneer in discovery of the solidliquid phase change phenomena, the Stefan number is a very important dimensionless variable in solid-liquid phase change phenomena, The Stefan number represents the ratio of sensible heat, c p1 (Tm − T0 ), to latent heat, hsA . For a latent heat thermal energy storage system, the Stefan number is usually very small because the temperature difference in such a system is small, while the latent heat hsA is very high. Therefore, the effect of the sensible heat transfer on the motion of the solid-liquid interface is very weak, and various approximate solutions to the phase change problem can be introduced without incurring significant error. It can be seen from eq. (6.31) that the effect of heat conduction in the liquid phase can be neglected if Ti − Tm  Tm − T0 or k2  k1 . In that case, eq. (6.31) can be simplified as s2  Ste (6.33) α1t Ste = c p1 (Tm − T0 ) 428 Transport Phenomena in Multiphase Systems which means that the interfacial velocity increases with increasing ΔT = Tm − T0 or decreasing hsA . 6.3.2 Dimensionless Form of the Governing Equations The governing eqs. (6.20) – (6.26) can be nondimensionalized by introducing the following dimensionless variables: T −T T −T α t½ x s S= θ= m θi = m i X = τ = 12 ° ° Tm − T0 Tm − T0 L L L° ° ¾ (6.34) ° c p1 (Tm − T0 ) α2 k2 ° Ste = Nα = Nk = ° ° α1 k1 hsA ¿ where L in eq. (6.34) is a characteristic length of the problem and can be determined by the nature of the problem or requirement of the solution procedure. The dimensionless governing equations are as follows: ∂ 2θ1 ∂θ1 = 0 < X < S (τ ), τ > 0 (6.35) ∂X 2 ∂τ θ1 ( X ,τ ) = 1 X = 0, τ > 0 (6.36) ∂ 2θ 2 1 ∂θ 2 = S (τ ) < X < ∞, τ > 0 2 ∂X Nα ∂τ θ 2 ( X ,τ ) → θi X → ∞, τ > 0 θ 2 ( X ,τ ) = θ i X > 0, τ = 0 θ1 ( X ,τ ) = θ 2 ( X ,τ ) = 0 X = S (τ ), τ > 0 ∂θ ∂θ 1 dS X = S (τ ), τ > 0 − 1 + Nk 2 = ∂X ∂X Ste dτ θ 1.0 (6.37) (6.38) (6.39) (6.40) (6.41) θ1(x,τ) s(τ) θi θ2(x,τ) x Figure 6.5 Dimensionless temperature distribution in the PCM. Chapter 6 Melting and Solidification 429 Dimensionless temperature distribution in a PCM can be qualitively illustrated by Fig. 6.5. It can be seen that the dimensionless temperature distribution is similar to that of a melting problem. Equations (6.35) – (6.41) are also valid for melting problems, provided that the subscripts “1” and “2” represent liquid and solid, respectively. The following solutions of one-region and two-region problems will be valid for both melting and solidification problems. 6.3.3 Exact Solution of the One-Region Problem If the initial temperature of a PCM equals its melting point of the PCM ( Ti = Tm ) , the melting or solidification problem is a one-region problem referred to as a Stefan problem. The temperature distribution in one phase is unknown while the temperature in the other phase remains at the melting point. Therefore, the temperature in only one phase needs to be solved. The mathematical description of a Stefan problem can be obtained by simplifying eqs. (6.35) – (6.41), i.e., ∂ 2θ ∂θ = 0 < X < S (τ ), τ > 0 (6.42) ∂X 2 ∂τ θ ( X ,τ ) = 1 X = 0, τ > 0 (6.43) θ ( X ,τ ) = 0 X = S (τ ), τ > 0 (6.44) ∂θ 1 dS X = S (τ ), τ > 0 − = (6.45) ∂X Ste dτ where the subscript “1” for the liquid phase has been dropped for ease of notation. The temperature distribution of this problem can be constructed on the basis of the heat conduction solution in a semi-infinite body (Ozisik, 1993), i.e., θ ( X ,τ ) = 1 + Berf ( X / 2τ 1/ 2 ) (6.46) where B is an unspecified constant. The error function in eq. (6.46) is defined as: 2 z − z2 erf ( z ) = (6.47) ³ e dz π 0 It can be verified that eq. (6.46) satisfies eqs. (6.42) and (6.43). The constant B in eq. (6.46) can be determined by substituting eq. (6.46) into eq. (6.44), i.e., 0 = 1 + Berf ( S / 2τ 1/ 2 ) (6.48) Since B is a constant, S / 2τ 1/ 2 must also be a constant in order for eq. (6.48) to be satisfied. This constant can be represented by λ , so λ = S / 2τ 1/ 2 (6.49) Thus, 1 (6.50) B=− erf (λ ) 430 Transport Phenomena in Multiphase Systems By substituting eq. (6.50) into eq. (6.46), the dimensionless temperature distribution in phase 1 is obtained: erf ( X / 2τ 1/ 2 ) θ ( X ,τ ) = 1 − (6.51) erf (λ ) Substituting eq. (6.51) and eq. (6.49) into eq. (6.45), the following equation is obtained for the constant λ : 2 Ste λ eλ erf (λ ) = (6.52) The constant λ is a function of the Stefan number only; this constitutes the appeal of solving the problem by using the dimensionless form of governing equations. Once the constant λ is obtained, the temperature distribution θ ( X ,τ ) and the location of the solid-liquid interface S (τ ) can be obtained by eqs. (6.51) and (6.49) respectively. The exact solution of the one-region phase change problem can also be obtained by using a similarity solution (see Problem 6.8). Example 6.1 A solid PCM with a uniform initial temperature equal to the melting point, Tm , is in a half-space, x > 0. At time t = 0, a constant ′′ heat flux q0 is suddenly applied to the surface of the semi-infinite body. Assume that the densities of the PCM for both phases are the same, and that natural convection in the liquid phase is negligible. Find the transient location of the solid-liquid interface. Solution: The melting problem under consideration is shown in Fig. 6.6. Since the temperature of the solid phase is equal to the melting point and only the temperature in the liquid phase is unknown, this is a one-region π T ′′ q0 T(x,t) Ti=Tm s(t) Figure 6.6 Melting in a semi-infinite region under constant heat flux x Chapter 6 Melting and Solidification 431 melting problem. The governing equation, and the corresponding initial and boundary conditions of the problem are as follows ∂ 2T 1 ∂T 0 < x < s (t ), t > 0 (6.53) = ∂x 2 α ∂t ∂T ( x, t ) ′′ −k = q0 x = 0, t > 0 (6.54) ∂x T ( x, t ) = Tm x = s (t ), t > 0 (6.55) ∂T ds −k = ρ hsA x = s (t ), t > 0 (6.56) ∂x dt where the subscript 1 for liquid phase has been dropped to simplify the notation. Introducing the following dimensionless variables: ′′ c p (q0 L/k ) T − Tm x s αt (6.57) , X = , S = , τ = 2 , Ste = θ= ′′ q0 L/k L L L hsA where L is a characteristic length of the problem, eqs. (6.53) – (6.56) are nondimensionalized as ∂ 2θ ∂θ 0 < X < S (τ ), τ > 0 (6.58) = ∂X 2 ∂τ ∂θ ( X ,τ ) X = 0, τ > 0 (6.59) = −1 ∂X X = S (τ ), τ > 0 θ ( X ,τ ) = 0 (6.60) ∂θ 1 dS − = X = S (τ ), τ > 0 (6.61) ∂X Ste dτ It is assumed that the temperature distribution in the liquid phase region can be constructed in a fashion similar to the exact solution of conduction in a semi-infinite body with boundary conditions of the second kind: §X· θ ( X ,τ ) = A + 2 τ ierfc ¨ (6.62) ¸ ©2 τ ¹ where ierfc is a differential complementary error function: 1 − z2 ierfc( Z ) = e − Zerfc( Z ) (6.63) π The unknown variable A in eq. (6.62) can be obtained by substituting eq. (6.62) into eq. (6.60), i.e., §S· A = −2 τ ierfc ¨ (6.64) ¸ ©2 τ ¹ Thus, the temperature distribution in the liquid phase is ª §X· § S ·º θ ( X ,τ ) = 2 τ «ierfc ¨ (6.65) ¸ − ierfc ¨ ¸» ©2 τ ¹ © 2 τ ¹¼ ¬ 432 Transport Phenomena in Multiphase Systems Substituting eq. (6.65) into eq. (6.61) yields an ordinary differential equation about the location of the solid-liquid interface, i.e., dS §S· = Ste ⋅ erfc ¨ (6.66) ¸ dτ ©2 τ ¹ which is subject to the following initial condition: S (τ ) = 0 τ =0 (6.67) A closed-form expression of S (τ ) cannot be obtained from eq. (6.66) because the separation of variables in eq. (6.66) is impossible. Equation (6.66) can be solved numerically to determine the transient location of the solid-liquid interface. The above exact solution in dimensional form was first obtained by El-Genk and Cronenberg (1979). Cho and Sunderland (1981) pointed out that there is a contradiction in this solution, since eq. (6.62) does not satisfy the partial differential equation (6.58) unless A in eq. (6.62) is a constant. However, A in eq. (6.62) is not a constant because it is a function of dimensionless time τ [see eq. (6.64)]. The interested reader can find the detailed evaluation of El-Genk and Cronenberg’s (1979) solution in Cho and Sunderland (1981). Nevertheless, this is a first attempt to obtain the exact solution of melting/solidification in a semiinfinite body with boundary conditions other than those of the first kind. In addition to the melting/solidification under boundary conditions of the first and second kinds discussed above, the boundary condition of the third kind (convective heating or cooling at surface) is a more generalized boundary condition for many applications, such as freezing of biological materials. Melting and solidification under the boundary condition of the third kind can be solved by a procedure similar to that in Example 6.1 (see Problem 6.9). A very simplified solution, which will be given in Example 6.2, can provide the first order estimation to the phase change problem. Example 6.2 A slab of biological material with a thickness of 2L is cooled by a fluid with temperature of T∞ , and the convective heat transfer coefficient is h. It is assumed that the biological material can be treated as singlecomponent PCM with well defined melting point, Tm, and its initial temperature is at Tm. The frozen process is so slow that heat transfer in a frozen layer can be regarded as a pseudo-steady-state process. Estimate the time it takes for the entire slab to be frozen. Solution: Since both sides of the slab are cooled, only half of the problem needs to be considered. The latent heat released by the biological material due to freezing is ds q′′ = ρ hsA (6.68) dt Chapter 6 Melting and Solidification 433 Symmetric line Solid Liquid Tm T0 h , T’ 0 s L x Figure 6.7 Freezing of a slab under convective cooling. Assuming steady-state conduction in the frozen layer, the latent heat released by freezing must be conducted through the frozen layer, i.e., −1 s1 ′′ = (Tm − T∞ ) § + · q (6.69) ¨ ¸ ©k h¹ Combining eqs. (6.68) and (6.69), we have Tm − T∞ § s 1· dt = ¨ + ¸ ds (6.70) ρ hsA ©k h¹ The time it takes to freeze the entire slab, tf, can be obtained by integrating eq. (6.70) tf T − T L§ s 1· ∞ m (6.71) ³0 ρ hsA dt = ³0 ¨ k + h ¸ ds © ¹ i.e., ρ hsA § L2 L · (6.72) +¸ tf = ¨ Tm − T∞ © 2k h ¹ which can be nondimensionalized as 1 §1 1 · (6.73) τf = ¨+ ¸ Ste © 2 Bi ¹ where τ f = α t f / L2 , Ste = c p (Tm − T∞ ) / hsA , Bi = hL / k are dimensionless time, Stefan number and Biot number. When Bi → ∞ , the problem 434 Transport Phenomena in Multiphase Systems becomes the boundary condition of the first kind [ T (0, t ) = T∞ ] and eq. (6.73) becomes 1 τf = (6.74) 2Ste 6.3.4 Exact Solution of the Two-Region Problem If the initial temperature of the PCM is not equal to its melting point of the PCM ( Ti ≠ Tm ) , the melting or solidification problem becomes a two-region problem, referred to as a Neumann problem in the literature. Temperature distributions in both the liquid and solid phases are unknown and must be solved. Equations (6.35) – (6.41) provide the complete mathematical description of a Neumann problem. Based on the heat conduction solution of a semi-infinite body, the temperature distribution in the PCM can be constructed as follows: θ1 ( X ,τ ) = 1 + Aerf(X / 2τ 1/ 2 ) (6.75) (6.76) where A and B in eqs. (6.75) and (6.76) are unspecified constants, and erfc is the complementary error function, defined as erfc( z ) = 1 − erf ( z ) (6.77) It should be noted that eqs. (6.75) and (6.76) satisfy eqs. (6.35) – (6.39). The constants A and B can be determined by using boundary condition (6.40), i.e., 1 + Aerf ( S / 2τ 1/ 2 ) = 0 (6.78) θ 2 ( X ,τ ) = θi + Berfc[ X / 2( Nατ 1/ 2 )] θi + Berfc[ S / 2( Nατ 1/ 2 )] = 0 (6.79) Since A and B are constants, S / 2τ 1/ 2 must also be a constant in order for eqs. (6.78) and (6.79) to be satisfied. This constant can be represented by λ , so λ = S / 2τ 1/ 2 (6.80) Thus, the constants A and B can be determined as 1 A=− erf (λ ) B=− (6.81) (6.82) θi 1 erfc(λ/Nα/ 2 ) Substituting eqs. (6.81) and (6.82) into eqs. (6.75) and (6.76), the temperature distributions in both phases are determined as follows: erf[ X / 2τ 1/ 2 ] θ1 ( X ,τ ) = 1 − (6.83) erf (λ ) θ 2 ( X ,τ ) = θi «1 − ¬ ª erf[ X / 2( Nατ )1/ 2 ] º » 1 erf (λ/Nα/ 2 ) ¼ (6.84) Chapter 6 Melting and Solidification 435 Substituting eqs. (6.83), (6.84), and (6.80) into eq. (6.41), the following equation is obtained for the constant λ : Nθ e−λ e − λ/Nα λπ + k1/ 2i = (6.85) 1/ 2 erf (λ ) Nα erfc(λ/Nα ) Ste Equation (6.85) can be solved by using an iterative method with underrelaxation because its left-hand side is a very complicated function of λ . Once λ is obtained, the temperature distributions – θ1 ( X ,τ ) and θ 2 ( X ,τ ) – and the location of the solid-liquid interface S (τ ) can be obtained from eqs. (6.83), (6.84), and (6.80), respectively. 2 6.4 Integral Approximate Solution The melting and solidification problems that can be solved by exact solution are very limited, so it is necessary to introduce some approximate solution techniques. The integral approximate method proposed by Goodman (1958) is one of the most attractive techniques, because it is very simple and its physical concept is very clear. After the integral approximate solution technique is introduced by solving heat conduction in a semi-infinite body, its application to various melting/solidification problems will be discussed. 6.4.1 Heat Conduction in a Semi-Infinite Body The early applications of integral approximate solutions to heat transfer problems included integral approximate solutions of boundary-layer momentum and energy equations. The method can also be used to solve linear or nonlinear transient conduction problems. The integral approximate solution will be employed to solve the heat conduction problem in a semi-infinite body with a boundary condition of the first kind. The physical model of the problem is illustrated in Fig. 6.8, and the governing equation of the heat conduction problem and the corresponding initial and boundary conditions are: ∂ 2T ( x, t ) 1 ∂T ( x, t ) = x > 0, t > 0 (6.86) ∂x 2 α ∂t T ( x, t ) = T0 x = 0, t > 0 (6.87) T ( x, t ) = Ti x > 0, t = 0 (6.88) For transient heat conduction in a semi-infinite body, the initial temperature is uniformly at Ti . At time t = 0, the surface temperature of the semi-infinite body is suddenly increased to a temperature T0 . The temperature near the surface of the semi-infinite body will increase as a result of the surface temperature change, while the temperature far from the surface of the semi-infinite body is not affected and remains at the initial temperature Ti . It is useful here to 436 Transport Phenomena in Multiphase Systems T T0 Ti δ(t) x Figure 6.8 Heat conduction in a semi-infinite body with constant wall temperature. introduce a concept similar to the thermal boundary layer for convective heat transfer, the thermal penetration depth. Assuming the thickness of the thermal penetration depth at time t is δ , the temperature of the semi-infinite body at x < δ will be affected but the temperature at x > δ will remain unchanged (see Fig. 6.8). It should be pointed out that the thermal penetration depth, δ , increases with increasing time. According to the definition of the thermal penetration depth, the temperature at the thermal penetration depth should satisfy ∂T ( x, t ) =0 (6.89) x = δ (t ) ∂x T ( x, t ) = Ti x = δ (t ) (6.90) Integrating eq. (6.86) in the interval ( 0, δ ), one obtains ∂T ∂x − x =δ ( t ) ∂T ∂x = x =0 α³ 1 δ (t ) 0 ∂T ( x, t ) dx ∂t (6.91) The right-hand side of eq. (6.91) can be rewritten using Leibnitz’s rule, i.e., ∂T ∂T 1ªd δ dδ º − =« (6.92) ³0 Tdx − T x =δ dt » ∂x x =δ ( t ) ∂x x = 0 α ¬ dt ¼ which represents the energy balance within the thermal penetration depth. Substituting eq. (6.89) and (6.90) into eq. (6.92) yields ∂T d −α = (Θ − Tiδ ) (6.93) ∂x x = 0 dt where ( ) Θ(t ) = ³ δ (t ) 0 T ( x, t )dx (6.94) Chapter 6 Melting and Solidification 437 Equation (6.93) is the integral energy equation of the conduction problem, and this equation pertains for the entire thermal penetration depth. It follows that a temperature distribution that satisfies eq. (6.93) does not necessarily satisfy differential eq. (6.86), which describes the energy balance at any and all points in the domain of the problem. In order to use the integral equation (6.93) to solve the conduction problem, the temperature profile in the thermal penetration depth is needed. The next step of the integral approximate solution is to assume a temperature distribution in the thermal penetration depth. The assumed temperature distribution can be any arbitrary function provided that the boundary conditions at x = 0 and x = δ are satisfied. Let us assume that the temperature distribution in the thermal penetration depth is a third-order polynomial function of x, i.e., T ( x, t ) = A0 + A1 x + A2 x 2 + A3 x 3 (6.95) where A0 , A1 , A2 , and A3 are four constants to be determined using the boundary conditions. Since there are only three boundary conditions available – eqs. (6.87), (6.89) and (6.90) – one more condition must be identified so that all four constants in eq. (6.95) can be determined. The surface temperature of the semi-infinite body, T0 , is not a function of time t , so ∂T ( x, t ) =0 x=0 (6.96) ∂t Substituting eq. (6.86) into eq. (6.96) yields ∂ 2T ( x, t ) =0 x=0 (6.97) ∂x 2 Substituting eq. (6.95) into eqs. (6.87), (6.89), (6.90) and (6.97) yields four equations for the constants in eq. (6.95). Solving for the four constants and substituting the results into eq. (6.95), the temperature distribution in the thermal penetration depth becomes 3 T ( x, t ) − Ti 3§ x · 1§ x · =1− ¨ ¸ + ¨ ¸ (6.98) 2©δ ¹ 2©δ ¹ T0 − Ti where the thermal penetration depth, δ , is still unknown. Substituting eq. (6.98) into eq. (6.93), an ordinary differential equation for δ is obtained: dδ 4α = δ (6.99) t >0 dt Since the thermal penetration depth equals zero at the beginning of the heat conduction, eq. (6.99) is subject to the following initial condition: δ =0 t =0 (6.100) The solution of eqs. (6.99) and (6.100) is δ = 8α t (6.101) which is consistent with the result of scale analysis, δ  α t [see eq. (1.200)]. 438 Transport Phenomena in Multiphase Systems The temperature distribution in the thermal penetration depth can be obtained by eq. (6.98), and the temperature in the semi-infinite body beyond the thermal penetration depth equals the initial temperature, Ti . The temperature distribution in the thermal penetration depth obtained here, as well as the thermal penetration depth thickness, depend on the assumed temperature distribution. The degree of the polynomial function for the temperature distribution in the thermal penetration depth should not be higher than four, or the result obtained will oscillate around the actual temperature profile, giving erroneous results. From the above analysis, we can summarize the procedure of the integral approximate solution as follows: 1. Obtain the integral equation of the problem by integrating the partial differential equation over the thermal penetration depth. 2. Assume an appropriate temperature distribution – usually a polynomial function – in the thermal penetration depth, and determine the unknown constants in the polynomial by using the boundary conditions at x=0 and x= δ . Additional boundary conditions, if needed, can be obtained by further analysis of the boundary conditions and the conduction equation. 3. Obtain an ordinary differential equation of thermal penetration depth by substituting the temperature distribution into the integral equation. The thermal penetration depth thickness can be obtained by solving this ordinary differential equation. 4. The temperature distribution in the thermal penetration depth can be obtained by combining the thermal penetration depth thickness from step 3 into the temperature distribution from step 2. 6.4.2 One-Region Problem The one-region phase change problem to be studied here is the same as the oneregion problem in Section 6.3.1. The dimensionless governing equations for melting and solidification problems have the same form if the dimensionless variables defined in the preceding section are used. The procedure for solving one-region phase change problems will be demonstrated by solving a one-region melting problem. The dimensionless governing equations of the melting problem are eqs. (6.42) – (6.45). To solve the problem by integral approximate solution, the thermal penetration depth must be specified. For a one-region melting problem, only the temperature distribution in the liquid phase needs to be solved, because the temperature in the solid phase remains uniformly equal to the melting point of the PCM. Therefore, the thickness of the liquid phase is identical to the thermal penetration depth. Integrating eq. (6.42) with respect to X in the interval of (0, S ) , one obtains ∂θ ( X ,τ ) ∂X − X = S (τ ) ∂θ ( X ,τ ) ∂X = X =0 d Θ(τ ) dτ (6.102) Chapter 6 Melting and Solidification 439 where Θ(τ ) = ³ S (τ ) 0 θ ( X ,τ ) dX (6.103) Substituting eq. (6.45) into eq. (6.102) yields the integral equation of the one-region problem: 1 dS ∂θ ( X ,τ ) d Θ(τ ) − − = (6.104) Ste dτ ∂X dτ X =0 Assume now that the temperature distribution in the liquid phase is the following second-degree polynomial: 2 § X −S · § X −S · (6.105) θ ( X ,τ ) = A0 + A1 ¨ ¸ + A2 ¨ ¸ ©S¹ ©S¹ where A0, A1, and A2 are three unknown constants to be determined. Equations (6.43) and (6.44) can be used to determine the constants in eq. (6.105). However, eq. (6.45) is not suitable for determining the constant in eq. (6.105) because dS /dτ in eq. (6.45) is unknown. Another appropriate boundary condition at the solid-liquid interface is needed. Differentiating eq. (6.44), one obtains ∂θ ∂θ dθ = dX + dτ = 0 X = S (τ ) (6.106) ∂X ∂τ i.e., ∂θ dS ∂θ X = S (τ ) + =0 (6.107) ∂X dτ ∂τ Substituting eqs. (6.42) and (6.45) into eq. (6.107) yields 2 ∂ 2θ § ∂θ · Ste¨ = (6.108) ¸ 2 © ∂X ¹ ∂X Equation (6.108) is an additional boundary condition at the solid-liquid interface; it can be used to determine the coefficients in eq. (6.105). Thus, the constants in eq. (6.105) can be determined by eqs. (6.43), (6.44) and (6.108). After the constants in eq. (6.105) are determined, the temperature distribution in the liquid phase becomes · § X − S ·2 1 − 1 + 2 Ste § X − S · § 1 − 1 + 2 Ste θ ( X ,τ ) = + 1¸ ¨ ¨ ¸+¨ ¸ © S ¸ (6.109) Ste Ste ©S¹¨ ¹ © ¹ Substituting eq. (6.109) into the integral equation (6.104) leads to an ordinary differential equation for S (τ ) : dS 1 − 1 + 2 Ste + 2 Ste =6 (6.110) dτ 5 + 1 + 2Ste + 2Ste The initial condition of eq. (6.110) is τ =0 (6.111) S (τ ) = 0 The right-hand side of eq. (6.110) is a constant and its solution is 1 S = 2λτ 2 (6.112) S 440 Transport Phenomena in Multiphase Systems 0.8 0.7 0.6 0.5 Integral Exact λ 0.4 0.3 0.2 0.1 0 0 0.4 0.8 1.2 1.6 2.0 2.4 2Ste Figure 6.9 Comparison between integral and exact solutions of the one-region conduction controlled melting problem. where ª 1 − 1 + 2 Ste + 2Ste º 2 λ = «3 (6.113) » ¬ 5 + 1 + 2Ste + 2 Ste ¼ Figure 6.9 shows the comparison of λ obtained by the exact solution and the integral approximate solution. The integral approximate solution agrees very well for a small Stefan number, but the difference increases along with increasing Stefan number. For a latent heat thermal energy storage system, the Stefan number is usually less than 0.2, so the integral approximate solution can provide sufficiently accurate results for that case. Example 6.3 A solid PCM with a uniform initial temperature at its melting point, ′′ Tm , is in a half-space, x > 0. At time t = 0, a variable heat flux, q0 (t ), is suddenly applied to the surface of the semi-infinite body. Assume that the densities of the PCM for both phases are the same, and that natural convection in the liquid phase is negligible. Find the transient location of the solid-liquid interface. Solution: This problem is the same as Example 6.1 except for the heat flux, which is ′′ q′′(t ) = f (t ) q0 (6.114) 1 Chapter 6 Melting and Solidification 441 ′′ where q0 is a reference heat flux and f (t ) is a given function. The dimensionless governing equation and the corresponding initial and boundary conditions of the problem are ∂ 2θ ∂θ 0 < X < S (τ ), τ > 0 = (6.115) ∂X 2 ∂τ ∂θ X = 0, τ > 0 = − f (τ ) (6.116) ∂X θ ( X ,τ ) = 0 (6.117) X = S (τ ), τ > 0 ∂θ 1 dS (6.118) X = S (τ ), τ > 0 − = ∂X Ste dτ where the nondimensional variables in the above eqs. (6.115) – (6.118) are defined by eq. (6.57). Since eq. (6.116) is a nonhomogeneous boundary condition, the exact solution in Example 6.1 cannot be directly applied here. The integral approximate solution, on the other hand, has the capability of dealing with nonlinear and non-homogeneous boundary conditions. Integrating eq. (6.115) with respect to X in the interval (0, S), and considering the boundary conditions at the surface, the integral equation of the problem is obtained: 1 dS (τ ) d Θ(τ ) − + f (τ ) = (6.119) Ste dτ dτ where Integrating eq. (6.119) with respect to τ in the interval (0, τ ) yields τ S (τ ) + Θ(τ ) = ³ f (τ )dτ (6.121) 0 Ste Assuming that the temperature distribution in the liquid phase is 2 § X −S · § X −S · (6.122) θ ( X ,τ ) = A0 + A1 ¨ + A2 ¨ ¸ ¸ ©S¹ ©S¹ where the constants, A0, A1, and A2 are three unspecified constants that can be determined from eqs. (6.116), (6.117) and (6.108). After all unknown constants are determined, the temperature distribution in the liquid phase becomes 1 «1 − (1 + 4 μ ) θ ( X ,τ ) = ¬ 2 Ste ª 1/ 2 º » ¼ Θ(τ ) = ³ S (τ ) 0 θ ( X ,τ ) dX (6.120) § X − S · 1 «1 − (1 + 4 μ ) ¬ ¨ ¸+ S¹8 Ste © ª 1/ 2 º 2 » ¼ § X −S · ¨ ¸ ©S¹ (6.123) 2 where μ = f (τ ) S (τ ) Ste Substituting eq. (6.123) into eq. (6.121), one obtains (6.124) 442 Transport Phenomena in Multiphase Systems (6.125) 6 For the case of constant heat flux, i.e., f (τ ) = 1, eq. (6.125) can be simplified to S ª SteS + 5 + (1 + 4SteS )1/ 2 º = 6 Steτ (6.126) « » ¬ ¼ 0 Ste2 f (τ ) ³ f (τ ) dτ = τ μª « ¬ μ + 5 + (1 + 4μ )1/ 2 º » ¼ 6.4.3 Two-Region Problem The physical model of a melting problem is shown in Fig. 6.10: a solid PCM with a uniform initial temperature Ti, which is below its melting point Tm , is in a ′′ half-space, x > 0. At time t = 0, a constant heat flux, q0 , is suddenly applied to the surface of the semi-infinite body. Because the initial temperature of the PCM is below its melting point, melting does not begin until after the wall temperature reaches the melting point. Therefore, the problem can be divided into two subproblems: (1) heat conduction over the duration of preheating, and (2) the actual melting process. An integral approximate method will be employed to solve both sub-problems (Zhang et al., 1993). Duration of Preheating At the beginning of heating, no melting occurs, and the problem is a pure conduction problem with a boundary condition of the second kind. Its mathematical description is as follows: ∂ 2T2 1 ∂T2 0 < x < ∞, 0 < t < tm = (6.127) 2 α 2 ∂t ∂x T T1(x,t) ′′ q0 T2(x,t) T2(x,tm) Tm s(t) δm δ(t) x Figure 6.10 Melting in a subcooled semi-infinite body under constant heat flux Chapter 6 Melting and Solidification 443 ∂T2 1 ′′ = − q0 x = 0, 0 < t < tm (6.128) k2 ∂x T2 ( x, t ) → Ti x → ∞, 0 < t < tm (6.129) T2 ( x, t ) = Ti , 0 < x < ∞, t = 0 (6.130) where tm is the duration of preheating. This problem is solved by integral approximate method. Assuming the temperature profile is a second-degree polynomial, one obtains the temperature profile: 2 ′′ q0δ § x· T2 ( x, t ) = Ti + (6.131) ¨1 − ¸ 2k2 © δ ¹ where δ is the thermal penetration depth. It can be obtained by substituting eq. (6.131) into the integral equation, which can in turn be obtained by integrating eq. (6.127) in the interval of (0, δ ). The result is (6.132) The highest temperature of the semi-infinite body occurs at its surface (x=0) and can be expressed as q′′δ Ts (t ) = Ti + 0 (6.133) 2k 2 Melting occurs when the surface temperature reaches the melting point, Tm , and the corresponding thermal penetration depth is 2k (T − T ) (6.134) δm = 2 m i ′′ q0 Then the duration of preheating can be obtained by combining eqs. (6.132) and (6.134), 2k 2 (T − T ) 2 tm = 2 m 2 i (6.135) ′′ 3α 2 q0 The temperature distribution at time tm is § x· T2 ( x, tm ) = Ti + (Tm − Ti )¨ 1 − ¸ © δm ¹ Governing Equations for the Melting Stage 2 δ = 6α 2t (6.136) After melting starts, the governing equations in the different phases must be specified separately. The temperature in the liquid phase satisfies ∂ 2T1 1 ∂T1 = (6.137) 0 < x < s (t ), t > tm ∂x 2 α1 ∂t ∂T1 ( x, t ) 1 ′′ = − q0 x = 0, t > tm (6.138) ∂x k1 444 Transport Phenomena in Multiphase Systems The governing equation and corresponding boundary and initial conditions for the solid phase are ∂ 2T2 1 ∂T2 s (t ) < x < ∞, t > tm = (6.139) 2 ∂x α 2 ∂t T2 ( x, t ) → Ti x → ∞, t > tm (6.140) § x· (6.141) T2 ( x, t ) = Ti + (Tm − Ti )¨1 − ¸ x > 0, t = tm © δm ¹ At the solid-liquid interface, the following boundary conditions are necessary to link solutions in the liquid and solid phases: T1 ( x, t ) = T2 ( x, t ) = Tm x = s (t ), t > tm (6.142) ∂T ∂T ds (6.143) k2 2 − k1 1 = ρ hsA x = s (t ), t > tm ∂x ∂x dt By defining the dimensionless variables as follows: c (T − T ) c (T − T ) c p 2 (Tm − Ti ) ½ ° ° Sc = θ1 = p1 1 m θ2 = p2 2 m ° hsA hsA hsA ° 2 x X= ′′ α1 ρ1hsA /q0 Δm = s S= ′′ α1 ρ1hsA /q0 Δ= δ δm ′′ α1 ρ1hsA /q0 τ= α1t ′′ (α1 ρ1hsA /q0 ) 2 ′′ α1 ρ1hsA /q0 α Nα = 2 α1 ° ° ° ¾ ° ° ° ° ° ° ° ¿ (6.144) where Sc is subcooling parameter that signifies the ratio of the sensible preheat need to raise the temperature of the body before melting takes place, eqs. (6.137) – (6.143) become ∂ 2θ1 ∂θ1 0 < X < S (τ ), τ > τ m = (6.145) ∂X 2 ∂τ ∂θ1 ( X ,τ ) = −1 X = 0, τ > τ m (6.146) ∂X ∂ 2θ 2 1 ∂θ 2 S (τ ) < X < ∞, τ > τ m = (6.147) 2 Nα ∂τ ∂X θ 2 ( X ,τ ) → − Sc X → ∞, τ > τ m (6.148) 2 ª§ º X· θ 2 ( X ,τ ) = Sc «¨1 − X > 0, τ = τ m ¸ − 1» «© Δ m ¹ » ¬ ¼ θ1 ( X ,τ ) = θ 2 ( X ,τ ) = 0 X = S (τ ), τ > τ m ∂θ 1 ∂θ 2 dS − 1+ = X = S (τ ), τ > τ m ∂X Nα ∂X dτ where Δ m and τ m can be obtained by substituting eq. (6.144) into and (6.135), i.e., (6.149) (6.150) (6.151) eqs. (6.134) Chapter 6 Melting and Solidification 445 Δ m = 2 Nα Sc 2 τ m = Nα Sc 2 3 Integral Approximate Solution (6.152) (6.153) Figure 6.11 shows the physical model represented by the above dimensionless governing equations. After the melting process starts, the solid-liquid interface moves along the positive x-axis, and the thermal penetration depth continuously increases. The solid phase temperature over the range of S (τ ) < X < Δ (τ ) is affected by the boundary condition at the surface. The temperature in the region beyond the thermal penetration depth, Δ, is not affected and remains at –Sc. Integrating eq. (6.147) over the interval ( S , Δ ), then applying the definition of the thermal penetration depth and eq. (6.150), yields the integral equation of the solid phase: ∂θ 1d −2 = (6.154) ( Θ2 + ScΔ ) ∂X X = S Nα dτ where Θ2 = ³ θ 2 dX S Δ (6.155) The temperature distribution in the solid phase is assumed to be a seconddegree polynomial, i.e., θ 2 = B0 + B1 ( X − S ) + B2 ( X − S ) 2 (6.156) Figure 6.11 Dimensionless temperature distribution for melting in a semi-infinite body under constant heat flux. 446 Transport Phenomena in Multiphase Systems The constants in eq. (6.156) can be obtained from the boundary condition, eq. (6.150) and the definition of the thermal penetration depth ( θ 2 X =Δ = − Sc and ∂θ 2 / ∂X phase is ª§ Δ − X · 2 º (6.157) ¸ − 1» «© Δ − S ¹ » ¬ ¼ Substituting eq. (6.157) into eq. (6.154) yields a relationship between the location of the solid-liquid interface and the thermal penetration depth: 6 Nα dS d Δ =2 + (6.158) dτ dτ Δ−S Integrating the differential equation of liquid phase eq. (6.145) over the interval of (0, S) and applying the boundary conditions eqs. (6.146) and (6.150) yields the integral equation of the liquid phase: ∂θ1 d Θ1 +1 = (6.159) dτ ∂X X = S where X =Δ = 0). The final form of the temperature distribution in the solid θ 2 ( X ,τ ) = Sc «¨ Θ1 = ³ θ1dX 0 S (6.160) Assuming that the temperature profile in the liquid phase is a second-degree polynomial function, § X −S · § X −S · (6.161) θ1 = A0 + A1 ¨ ¸ + A2 ¨ ¸ ©S¹ ©S¹ where A0, A1, and A2 need to be determined using appropriate boundary conditions. Equations (6.146) and (6.150) provided two conditions; and the third condition needs to be obtained by differentiating eq. (6.150), i.e., ∂θ ∂θ dθ1 = 1 dX + 1 dτ = 0 X = S (τ ) (6.162) ∂X ∂τ i.e., ∂θ1 dS ∂θ1 + =0 X = S (τ ) (6.163) ∂X dτ ∂τ Substituting eqs. (6.145) and (6.151) into eq. (6.163) yields 2 1 ∂θ1 ∂θ 2 ∂ 2θ1 § ∂θ1 · = (6.164) ¨ ¸− 2 © ∂X ¹ Nα ∂X ∂X ∂X which is an additional boundary condition at the solid-liquid interface, and can be used to determine the coefficients in eq. (6.161). After the constants in eq. (6.161) are determined, the temperature distribution in the liquid phase becomes 2 θ1 ( X ,τ ) = S § X − S · 1 X 2 − S2 ¨ ¸− p 2© S ¹ 2 S2 2 (6.165) Chapter 6 Melting and Solidification 447 where § N ScS 1 · § N ScS 1 · p=¨ α − ¸+ ¨ α − ¸ +S (6.166) © Δ−S 2¹ © Δ−S 2¹ Substituting the integral eqs. (6.159) and (6.154) into eq. (6.151) yields dΘ dS d (6.167) + ( Θ2 + ScΔ ) + 1 = 1 dτ dτ dτ Integrating both sides of eq. (6.167) with respect to τ within ( τ m ,τ ), one obtains S + [Θ2 (τ ) − Θ 2 (τ m )] + Sc(Δ − Δ m ) + [Θ1 (τ ) − Θ1 (τ m )] = (τ − τ m ) (6.168) Substituting eqs. (6.155), (6.157), (6.160), and (6.165) into eq. (6.168), the resulting interface equation is 12 S + ( p + 3 + 2 Sc) S + Sc(Δ − Δ m ) = 3(τ − τ m ) (6.169) 2 Equation (6.158) is a first-order ordinary differential equation, while eq. (6.169) is an algebraic equation. They can be solved simultaneously with an implicit numerical method to determine the location of the solid-liquid interface and the thermal penetration depth. If there is no subcooling, Sc=0, the solidliquid interface location can be obtained by the following expression: § · S ¨ S + 5 + 1 + 4S ¸ = 6τ (6.170) ¨ ¸ © ¹ 2 Ablation under Constant Heat Flux Heating Ablation is defined as the removal of material from the surface of an object by vaporization, chipping, or other erosive process. One application is when tire manufacturers design tires for automotive race events, such as NASCAR races. They have to balance the expected tire wear on the track versus the amount of heat that will be removed via ablation. If they underestimate the wear, the tire overheats and fails as too little heat is removed. If they overestimate the wear, the tire temperature never reaches its optimal operating condition to maximize adhesion. Ablation is also an effective means of protecting the surfaces of missiles and space shuttles from high-rate aerodynamic heating during atmospheric re-entry. It is a sacrificial cooling method, because the protective layer is partially destroyed. The advantage of the ablative cooling process is its self-regulation: the rate of ablation is automatically adjusted in response to the heating rate. The most commonly used materials for ablative cooling are PCMs with higher melting points (such as glass, carbon, or polymer fiber) in combination with an organic binder. From the heat transfer point of view, ablation is a special case of melting in which the liquid phase is completely removed immediately upon its production. Therefore, ablation can be considered as a one-region melting problem where heating occurs directly on a solid-liquid interface. If the spacecraft maintains a 448 Transport Phenomena in Multiphase Systems constant velocity during re-entry, one can assume that the entire process occurs under constant heat flux heating. Similar to the melting of a subcooled solid under constant heat flux, ablation does not start simultaneously with heating, because the initial temperature of the ablating material, Ti , is usually below its melting point, Tm. The melting will start only after a period of preheating that allows the surface temperature to reach the melting point of the ablating material, Tm. The preheating problem is a pure conduction problem and can be solved analytically. Since the surface of the ablating material is constantly moving and the liquid phase does not accumulate, ablation can be described by using a coordinate system, x, that is attached to and moves with the ablating velocity, Ua (see Fig 6.12). The ablation material moves with a velocity –Ua in the moving coordinate system. The energy equation in a moving coordinate is (Eckert and Drake, 1987): ∂ § ∂T · ∂T ∂T (6.171) + ρcp ¨k ¸ = −U a ρ c p ∂x © ∂x ¹ ∂x ∂t Since the ablating materials usually have very low thermal conductivity, it is reasonable to assume that the ablation occurs in a semi-infinite body. Therefore, the initial and boundary conditions of eq. (6.171) are T = Ti , t = 0 (6.172) T = Tm , x = 0 (6.173) ∂T T = Ti , = 0, x→∞ (6.174) ∂x The energy balance at the surface is ∂T , x=0 (6.175) q′′ − ρU a hsA = −k ∂x x = 0 Figure 6.12 Physical model for ablation. Chapter 6 Melting and Solidification 449 Shortly after ablation begins, the process enters steady-state and the ablation velocity becomes a constant. Equation (6.171) can be simplified as U dT d 2T =− a (6.176) 2 dx α dx where the thermal properties are assumed to be independent of temperature. Equation (6.176) can be treated as a first-order differential equation of dT/dx, and its solution is U − ax dT = C1e α (6.177) dx where C1 is an integrating constant. Equation (6.177) can be further integrated to yield U C1α − αa x T =− e + C2 (6.178) Ua The constants C1 and C2 can be determined from eqs. (6.173) – (6.174), and the temperature distribution in the ablating material becomes T = Ti + (Tm − Ti )e α Substituting eq. (6.179) into eq. (6.175) yields q′′ − ρU a hsA = ρ c pU a (Tm − Ti ) − Ua x (6.179) (6.180) which can be rearranged to obtain the ablation velocity: q′′ (6.181) Ua = ρ hsA (1 + Sc) where Sc is the subcooling parameter defined as c p (Tm − Ti ) (6.182) Sc = hsA The fraction of heat removed by ablation over the total heat flux is ′′ qa ρU a hsA = (6.183) q′′ q′′ Substituting the ablation velocity from eq. (6.181) into eq. (6.183), one obtains ′′ qa 1 = (6.184) q′′ 1 + Sc If the subcooling parameter is Sc = 0.2, 83% of the total heat flux will be removed by ablation. 6.4.5 Solidification/Melting in Cylindrical Coordinate Systems Application of the integral approximate solution to one-dimensional solid-liquid phase change problems – including ablation – in a Cartesian coordinate system has been discussed in the preceding sections. Since tubes are widely used in 450 Transport Phenomena in Multiphase Systems shell-and-tube thermal energy storage devices, it is necessary to study melting and solidification in cylindrical coordinate systems as well. The polynomial temperature distribution is a very good approximation of the one-dimensional problem in the Cartesian coordinate system, but it can result in very significant error if it is used to solve for the phase change heat transfer in a cylindrical coordinate system. This is because heat transfer area for a cylindrical coordinate system varies with the coordinate r instead of remaining constant. Therefore, the temperature distribution in the coordinate has to be modified by taking into account the variation of the heat transfer area. Solidification around a cylinder with radius of ri, as shown in Fig. 6.13, will be investigated in order to demonstrate application of the integral approximate solution in the cylindrical coordinate system. An infinite liquid PCM has a uniform initial temperature equal to the melting point of the PCM, Tm . At time t = 0, the inner surface of the cylinder suddenly decreases to a temperature T0 , which is below the melting point of the PCM. Since the temperature of the liquid PCM equals the melting point of the PCM and the temperature in the solid phase is unknown, this is a one-region solidification problem. Conduction controls the solidification process, because the temperature in the liquid phase is uniformly equal to the melting point of the PCM. The governing equation and the initial and boundary conditions of this problem are 1 ∂ § ∂T · 1 ∂T (6.185) ri < r < s (t ) t > 0 ¨r ¸= r ∂r © ∂r ¹ α ∂t Figure 6.13 Solidification of an infinite liquid PCM around an ID-cooled cylinder. Chapter 6 Melting and Solidification 451 T (r , t ) = T0 r = ri t > 0 (6.186) T (r , t ) = Tm r = s t > 0 (6.187) ∂T ds = ρ hsA k r =s t >0 (6.188) ∂r dt where the subscript 1 for solid phase has been dropped for ease of notation. Defining dimensionless variables, c (T − T ) T −T r s αt R= S= (6.189) θ= m τ = 2 Ste = p m 0 Tm − T0 ri ri ri hsA eqs. (6.185) – (6.188) become 1 ∂ § ∂θ · ∂θ (6.190) 1 < R < S (τ ) τ > 0 ¨R ¸= R ∂R © ∂R ¹ ∂τ θ ( R, τ ) = 1 R = 1 τ > 0 (6.191) θ ( R,τ ) = 0 R = S (τ ) τ > 0 (6.192) ∂θ 1 dS − = R=S τ >0 (6.193) ∂R Ste dτ The above eqs. (6.190) – (6.193) are also valid for melting around a hollow cylinder when used with the appropriate dimensionless variables. From elementary heat transfer, it is know that the logarithmic function is the exact solution of the steady-state heat conduction problem in a cylindrical wall. Therefore, we can assume that the temperature distribution has a second-order logarithmic function of the form (Zhang and Faghri, 1996a): § ln R · § ln R · (6.194) ¸ − (1 + ϕ )¨ ¸ © ln S ¹ © ln S ¹ where ϕ is an unknown variable. Equation (6.194) automatically satisfies eqs. (6.191) and (6.192), and ϕ can be obtained by differentiating eq. (6.192) ∂θ ∂θ dθ = dR + dτ = 0, R = S ∂R ∂τ which can be rearranged to get ∂θ dS ∂θ + =0 (6.195) ∂R dτ ∂τ Substituting eqs. (6.190) and (6.192) into eq. (6.195) yields the following expression: 2 θ =1+ϕ ¨ ϕ: § ∂θ · 1 ∂ § ∂θ · (6.196) − Ste¨ ¸+ ¨R ¸=0 R=S τ >0 © ∂R ¹ R ∂R © ∂R ¹ Substituting eq. (6.194) into eq. (6.196) gives the following expression for 2 2 +ϕ = 1 + 2Ste − 1 Ste (6.197) 452 Transport Phenomena in Multiphase Systems Substituting eqs. (6.194) and (6.197) into eq. (6.193) leads to an ordinary differential equation for the locations of the solid-liquid interface dS 1 + 2Ste − 1 = (6.198) dτ S ln S which is subjected to the following initial condition: S (τ ) = 1, τ =0 (6.199) Integrating equation (6.198) over the time interval (0,τ ) results in the following equation for the location of the solid-liquid interface: 12 1 (6.200) S ln S − ( S 2 − 1) = 1 + 2 Ste − 1 τ 2 4 which is valid for solidification around a cylinder with a constant inner temperature. However, in practical applications, cooling inside the tube is usually achieved by the flow of cooling fluid through the tube. Therefore, the boundary condition at the inner surface of the tube should be a convective boundary condition instead of an isothermal boundary condition. The following example demonstrates application of integral approximate method to the solution of solidification around a cylindrical tube convectively cooled from inside. ( ) Example 6.4 An infinite liquid PCM has a uniform initial temperature equal to the melting point of the PCM, Tm . At time t = 0, a cooling fluid with temperature Ti flows inside the tube. The heat transfer coefficient between the cooling fluid and the inner surface of the tube is hi . Find the transient location of the solid-liquid interface. Solution: The governing equations of the problem can also be represented by eq. (6.185) – (6.188) except that eq. (6.186) needs to be replaced by the following expression: ∂T k = hi (T − Ti ) r = ri t > 0 (6.201) ∂r The dimensionless governing equations of the problem are the same as eqs. (6.190) – (6.193), but eq. (6.191) needs to be replaced by the dimensionless form of eq. (6.201), i.e., ∂θ = Bi (θ − 1), R = 1 τ = 0 (6.202) ∂R where Bi in eq. (6.202) is the Biot number defined as hr Bi = i (6.203) k and θ is defined as T −T θ= m (6.204) Tm − Ti Chapter 6 Melting and Solidification 453 It is also assumed that the temperature distribution has a secondorder logarithmic function of the form § ln R · § ln R · θ = A + B¨ (6.205) ¸ − ( A + B)¨ ¸ © ln S ¹ © ln S ¹ The two unknown variables A and B in eq. (6.205) can be determined by substituting eq. (6.205) into eqs. (6.202) and (6.196), with the result that B = − Bi (1 − A) ln S (6.206) 1 + 2 ASte − 1 (6.207) Ste Substituting eq. (6.206) into eq. (6.207), an equation for A is obtained as follows: 1 + 2 ASte − 1 2 A − Bi (1 − A) ln S = (6.208) Ste Substituting eq. (6.205) into eq. (6.193), the ordinary differential equation for the location of the solid-liquid interface is obtained: dS 1 + 2 ASte − 1 = (6.209) dτ S ln S which is subject to the initial condition specified by eq. (6.199). The temperature of the inner surface of the tube is θ (1,τ ) = A (6.210) The solid-liquid interface location can be obtained by numerical solution of eq. (6.209) and (6.210). It should be noted that A becomes 1 if the Biot number becomes infinite [see eq. (6.208)]. In that case the temperature of the inner surface becomes 1 and eq. (6.209) reduces to eq. (6.198). 2A + B = 2 6.5 Numerical Simulation 6.5.1 Overview Due to the strong nonlinearity of solid-liquid phase change phenomena and the moving boundary, the problems that can be solved via the analytical method are very limited. Exact solutions and some approximate solutions for melting and solidification have been introduced in preceding sections. The limitation of these methods is that they apply to conduction-controlled one-dimensional problems. Although some investigators have attempted to solve two-dimensional phase change problems by using analytical methods (Poots, 1962; Rathjen and Jiji, 1971; Budhia and Kreith, 1973), the cases investigated represented very simple and special geometries, such as a semi-infinite corner or a semi-infinite wedge. 454 Transport Phenomena in Multiphase Systems Generally, the analytical method will not work for two- or three-dimensional problems; this is particularly true for convection-controlled phase change processes. Most real applications of phase change occur in complex two- or three-dimensional geometries. In addition, natural convection in the liquid phase often plays a significant role in the phase-change processes. For such complex solid-liquid phase-change processes, analytical solutions will not work. Therefore, it is necessary to employ numerical methods. A large number of numerical methods have been developed and they can be divided into two groups. The first group is called strong numerical solution, or classical solution. For this group, transformed coordinate systems or deformed grids are employed to deal with the location of the solid-liquid interface. In this methodology, complex geometric regions of solid or liquid phases are transformed into fixed, simple geometric regions through the coordinate transformation technique. At the same time, the governing equations and the boundary conditions become more complicated. The governing equations of the problem can be solved using the finite difference method. Such methods also deal successfully with multidimensional one-region or two-region problems with or without natural convection in the liquid phase. The disadvantage of these methods is that they are often difficult to program and thus require a significant amount of computer time. The second group is called the weak solution, or the fixed-grid solution. In these methods, the governing equations are written for the entire phase change region, including liquid and solid phases. The location of the solid-liquid interface is determined after the converged temperature distribution is obtained. The three principal important methods in this group are the enthalpy method (Shamsunder and Sparrow, 1975; Crank, 1984), the equivalent heat capacity method (Bonacina et al., 1973), and the temperature-transforming model (Cao and Faghri, 1990a). The temperature transforming model has the advantages of both the enthalpy method and the equivalent heat capacity method and can also account for the effect of natural convection in the liquid phase. The most significant advantage of the weak solution, as compared with the strong numerical solution, is its simplicity. In this section, various fixed-grid numerical solutions of the solid-liquid phase problem will be introduced. 6.5.2 Enthalpy Method In this methodology, the governing energy equation is written for the entire region of the PCM, including solid and liquid phases and the interface. The enthalpy method is introduced by analyzing a conduction-controlled, two-region melting problem in a finite slab, as shown in Fig. 6.14. It is assumed that the densities of the liquid and solid phase are identical ( ρ s = ρA ). The energy equation can be written as ∂h ∂ § ∂T · (6.211) ρ = ¨k ¸ ∂t ∂x © ∂x ¹ Chapter 6 Melting and Solidification 455 T T0 T(x,t) Tm Ti x=0 s x=L x Figure 6.14 Conduction-controlled melting in a finite slab. h h (T ) hsA T 0 Tm Figure 6.15 Enthalpy of the PCM as a function of temperature. The enthalpy, h, is a function of temperature, T: ­ T < Tm ° c ps (T − Tm ) ® h(T ) = ° ° c (T − T ) + h m sA T > Tm ° pA ¯ (6.212) The variation of enthalpy h with temperature T can be plotted as shown in Fig. 6.15. At the melting point of the PCM, the enthalpies of solid and liquid phases at the melting point are 0 and hsA , respectively; there is a jump of enthalpy equal to the latent heat of the PCM. The thermal conductivity of the PCM, k, can be expressed as ­ T < Tm ° ks k (T ) = ® (6.213) °k ° A T > Tm ¯ 456 Transport Phenomena in Multiphase Systems It can be verified that the energy equation for liquid and solid phases can be obtained simply by substituting eqs. (6.212) and (6.213) into eq. (6.211). The energy balance at the solid-liquid interface can be obtained by analyzing the energy balance for a control volume, which includes the solid-liquid interface. Solving for temperature, T, from eq. (6.212) yields ­ h≤0 °Tm + h / c ps ® T = °Tm ° ° °m ¯ ° 0 < h < hsA h ≥ hsA (6.214) T + ( h − hsA ) / c pA The initial and the boundary conditions of the melting problem are T ( x, t ) = Ti < Tm t =0 (6.215) T ( x, t ) = T0 > Tm x=0 (6.216) ∂T ( x, t ) =0 x=L (6.217) ∂x The discretization of space and time is shown in Fig. 6.16, in which location and time are represented by j and n. The temperature at x = ( j − 1)Δx and time t = nΔt can be represented by the symbol T jn . Equation (6.211) is discretized by using an explicit scheme: n n n n h n +1 − h n k j + 1 (T j +1 − T j ) − k j − 1 (T j − T j −1 ) j j 2 2 = j = 2, ..., N − 1 (6.218) ρ 2 Δt ( Δx ) which uses forward difference in time and central difference in space. The thermal conductivities at the half-grid, k j + 1 and k j − 1 , are calculated 2 2 using the harmonic mean method: 2k j k j +1 k j+ 1 = 2 k j + k j +1 (6.219) Figure 6.16 Discretization of space and time for the solution of the two-region melting problem. Chapter 6 Melting and Solidification 457 k j−1 = 2 2k j −1k j k j −1 + k j (6.220) Equation (6.218) can be rewritten as Δt § · ¨k h n +1 = T n + k j − 1 T jn−1 ¸ ¸ j 2 ¨ j + 1 j +1 2 ¹ ρ ( Δx ) © 2 +h− n j § ¨ ¨ ¨ ¨ ¨ ¨ © Δt ¨ k j + 1 + k j − 1 ¸ ¨ ¸ © 2 2 § · ¹ ρ ( Δx ) 2 T · ¸ n¸ ¸ j¸ ¸ ¸ ¹ (6.221) j = 2, ..., N − 1 After the enthalpy distribution of the (n+1)th time step is obtained, the temperature distribution of the (n+1)th time zone can be obtained from eq. (6.214), i.e., ­ Tm + h n +1 / c ps h n +1 ≤ 0 (solid) j j ° n +1 n +1 Tj = ® Tm 0 < h j < hsA (interface) j = 2,..., N − 1 (6.222) n +1 °T + (h − h ) / c (liquid) h n +1 ≥ hsA j sA pA j ¯m The initial and boundary conditions are discretized as T j0 = Ti j = 1, 2, ..., N (6.223) T1n +1 = T0 n +1 N n +1 N −1 (6.224) T =T (6.225) The location of the solid-liquid interface at the (n+1)th time zone can easily be determined according to the temperature distribution. For example, if the n enthalpy at a certain point j = M satisfies 0 < hM+1 < hsA , the solid-liquid interface will be n Δx hM+1 + Δx s n +1 = ( M − 1)Δx − (6.226) 2 hsA The solution procedure can be summarized as follows: 1. Determine the initial enthalpy at every node h0 by using eq. (6.212). j 2. Calculate the enthalpy after the first time step at nodes j = 2, ..., N − 1 by using eq. (6.221). 3. Determine the temperature after the first time step at node j = 1, ..., N by using eqs. (6.224), (6.222), and (6.225). 4. Find a control volume in which the enthalpy falls between 0 and hsA and determine the location of the solid-liquid interface by using eq. (6.226). 5. Solve the phase-change problem at the next time step with the same procedure. The solution procedure using the explicit scheme is attractive in its simplicity. However, there is a major drawback of this scheme – the limitation of 458 Transport Phenomena in Multiphase Systems stability. The limitation of stability should be familiar from prior experience in undergraduate heat transfer dealing with the numerical solution of unsteady state heat conduction. To satisfy the stability condition, the following Neumann stability criterion has to be satisfied: max(α s ,α A )Δt 1 ≤ (6.227) 2 2 ( Δx ) where α s and α A are thermal diffusivity of solid and liquid, respectively. Equation (6.227) is valid for one-dimensional problems only. In additional to the above explicit scheme, one can use an implicit scheme – which is unconditionally stable – to overcome the drawback of potential instability of the explicit scheme. However, the solution process for an implicit scheme will be more complex than the solution for an explicit scheme because two unknown variables (enthalpy and temperature) are involved. The detailed procedure for the implicit scheme can be found in Crank (1984) or Alexiades and Solomon (1993). The above enthalpy model treated enthalpy as a dependent variable in addition to the temperature, and discretized the energy equation into a set of equations which contain both h and T. For the implicit schemes, they actually treated all of the terms containing T = T (h) as a constant heat source term in the energy equation during iterations at each time step. This may cause some problems with convergence when T = T ( h) is complicated and physical properties change significantly, as is the case for frozen heat pipe start-up, or when the boundary conditions are severe. Furthermore, when the energy equation contains a convective term, the previous methods have difficulties in handling the relationship between the convective and diffusive terms because of the twodependent-variable nature of the equation. Cao and Faghri (1989) proposed a simple strategy for transforming the energy equation into a nonlinear equation with a single dependent variable h. Thus, solving a phase-change problem is equivalent to solving a nonlinear enthalpy equation, and existing algorithms are readily applicable with some modifications. The energy equation governing three-dimensional incompressible laminar flow with no viscous dissipation and with incorporation of the continuity equation in the Cartesian coordinate system is [see eq. (3.84)] ∂ ( ρ h) ∂ ( ρ uh) ∂ ( ρ vh) ∂ ( ρ wh) ∂ § ∂T · ∂ § ∂T · ∂ § ∂T · + + + = ¨k ¸+ ¨k ¸ + ¨k ¸ ∂t ∂x ∂y ∂z ∂x © ∂x ¹ ∂y © ∂y ¹ ∂z © ∂z ¹ (6.228) with the state equation dh = c(T ) (6.229) ∂T In the cases with constant specific heats for each phase, in which phase change occurs at a single temperature, the temperature and enthalpy are related Chapter 6 Melting and Solidification 459 by eq. (6.214). Substituting eq. (6.214) into eq. (6.228), the energy equation can be transformed into the following form (Cao and Faghri, 1989): ∂ ( ρ h) ∂ ( ρ uh) ∂ ( ρ vh) ∂ ( ρ wh) + + + ∂t ∂x ∂y ∂z (6.230) 2 2 2 ∂ ( Γ h ) ∂ ( Γh ) ∂ ( Γ h ) ∂ 2 S ∂ 2 S ∂ 2 S = + + + 2+ 2+ 2 ∂x 2 ∂y 2 ∂z 2 ∂x ∂y ∂z where h≤0 ­k s / c ps ° 0 < h < hsA Γ ( h ) = ®0 (6.231) °k / c h > hsA ¯ A pA ­0 h≤0 ° 0 < h < hsA (6.232) S ( h ) = ®0 °hk / c h > h sA ¯ A pA Cao and Faghri (1989) demonstrated that the above enthalpy transforming model is capable of handling complicated phase-change problems occurring both at a single temperature and over a temperature range with fixed grids. Due to the one dependent variable nature of the transformed equation, the convection and diffusion situations can be handled with appropriate algorithms. Comparisons have been made with the numerical results existing in the literature with a good agreement, showing that the model can properly predict the phase-change processes. The advantages of this model based on enthalpy are that it allows us to avoid paying explicit attention to the nature of the phase-change front, and that it can be extended to accommodate complicated multidimensional problems with convective terms without involving cumbersome mathematical schemes. 6.5.3 Equivalent Heat Capacity Method During the solid-liquid phase-change process, the PCM can absorb or release heat at a constant temperature Tm . This means that the temperature of the PCM does not change while it absorbs or releases heat, implying that the heat capacity of the PCM at temperature Tm is infinite. In the equivalent heat capacity method, it is assumed that the melting or solidification processes occur over a temperature range (Tm − ΔT , Tm + ΔT ) instead of at a single temperature Tm . For a multicomponent system, ΔT can be chosen based on the range of phase change temperature. For a single-component with well-defined melting point, ΔT should be as small as possible. Also, the latent heat is converted to an equivalent heat capacity of the PCM in the assumed temperature range. Thus, the specific heat of the PCM can be expressed as 460 Transport Phenomena in Multiphase Systems T < Tm − ΔT ­c ps ° c ps + c pA °h (6.233) c p (T ) = ® sA + Tm − ΔT < T < Tm + ΔT 2ΔT 2 ° T > Tm + ΔT ° c pA ¯ which assumes that the temperature of the PCM is changed from Tm − ΔT to Tm + ΔT when latent heat, hsA , is absorbed by the PCM during melting. During the solidification process, the PCM releases the latent heat and its temperature decreases from Tm + ΔT to Tm − ΔT . The equivalent specific heat in the mushy zone ( Tm − ΔT < T < Tm + ΔT ) includes the effect of both latent heat (the first term) and sensible heat (the second term). The relationship between specific heat and temperature in the equivalent heat capacity method is plotted in Fig. 6.17. For a three-dimensional conductioncontrolled melting/solidification problem in the Cartesian coordinate system, the energy equation for the entire region of the PCM can be expressed as ∂T ∂ § ∂T · ∂ § ∂T · ∂ § ∂T · ρcp (6.234) = ¨k ¸+ ¨k ¸ + ¨k ¸ ∂t ∂x © ∂x ¹ ∂y © ∂y ¹ ∂z © ∂z ¹ where the thermal conductivity in eq. (6.234) is a function of temperature, T. Assuming that the thermal conductivity of the PCM in the two-phase region is a linear function of temperature, one obtains T < Tm − ΔT ­ks ° k − ks ° (T − Tm + ΔT ) Tm − ΔT < T < Tm + ΔT (6.235) k (T ) = ®k s + A 2ΔT ° T > Tm + ΔT ° kA ¯ Figure 6.17 Dependence of specific heat to temperature for equivalent heat capacity model. Chapter 6 Melting and Solidification 461 The advantage of the equivalent heat capacity model is its simplicity. Equation (6.234) is simply the nonlinear heat conduction equation, and it appears that a conventional computational methodology for conduction problems is adequate for solving a solid-liquid phase change problem. However, many studies have revealed difficulties in the selection of time step, t, grid size, ( x, y, z) and the phase change temperature range, T. If these variables cannot be properly selected, the predicted location of the solid-liquid interface and the temperature may include some unrealistic oscillation. Therefore, although the equivalent heat capacity model leads to simple code development, it is not used as widely as the enthalpy model. 6.5.4 Temperature-Transforming Model The temperature-transforming model proposed by Cao and Faghri (1990a) combines the advantages of the enthalpy and equivalent heat capacity models. For a three-dimensional conduction-controlled phase change problem, the governing equation in enthalpy form is ∂h ∂ § ∂T · ∂ § ∂T · ∂ § ∂T · ρ = ¨k (6.236) ¸ + ¨k ¸ + ¨k ¸ ∂t ∂x © ∂x ¹ ∂y © ∂y ¹ ∂z © ∂z ¹ For a phase change occurring over a temperature range ( Tm − ΔT , Tm + ΔT ) with the specific heats assumed to be constant for each phase, the relationship between enthalpy and temperature can be plotted as in Fig. 6.18. This relationship can be analytically expressed as (see problem 6.28) ­ c ps (T − Tm ) + c ps ΔT T < Tm − ΔT ° c ps + c pA h· h °§ c ps + c pA h(T ) = ®¨ + sA ¸ (T − Tm ) + ΔT + sA Tm − ΔT < T < Tm + ΔT 2 2ΔT ¹ 2 2 °© ° c pA (T − Tm ) + c ps ΔT + hsA T > Tm + ΔT ¯ (6.237) By defining specific heat in the mushy zone as c ps + c pA cm = 2 eq. (6.237) becomes ­c ps (T − Tm ) + c ps ΔT ° h· h °§ h(T ) = ®¨ cm + sA ¸ (T − Tm ) + cm ΔT + sA 2ΔT ¹ 2 °© °c pA (T − Tm ) + c ps ΔT + hsA ¯ (6.238) T < Tm − ΔT Tm − ΔT < T < Tm + ΔT (6.239) T > Tm + ΔT (6.240) Equation (6.239) can be rewritten as h(T ) = c p (T )(T − Tm ) + b(T ) 462 Transport Phenomena in Multiphase Systems Figure 6.18 Dependence of enthalpy on temperature for phase change occurring over a range of temperature. where c p (T ) and b(T ) can be determined from eq. (6.239): T < Tm − ΔT ­c ps ° h ° c p (T ) = ®cm + sA Tm − ΔT < T < Tm + ΔT (6.241) 2ΔT ° T > Tm + ΔT ° c pA ¯ T < Tm − ΔT ­c ps ΔT ° h ° b(T ) = ®cm ΔT + sA Tm − ΔT < T < Tm + ΔT (6.242) 2 ° T > Tm + ΔT °c ps ΔT + hsA ¯ Substituting eq. (6.240) into eq. (6.236) yields ∂ (c pT ) ∂ § ∂T · ∂ § ∂T · ∂ § ∂T · ∂b (6.243) = ¨k ρ ¸ + ¨k ¸ + ¨k ¸−ρ ∂t ∂x © ∂x ¹ ∂y © ∂y ¹ ∂z © ∂z ¹ ∂t where the thermal conductivity, k, is a function of temperature and can be obtained by eq. (6.235). The energy equation has been transformed into a nonlinear equation with a single dependent variable T. The nonlinearity of the phase change problem results from the large amount of heat released or absorbed at the solid-liquid interface. Also, in either the liquid or solid phases at some distance away from the interface, eq. (6.243) reduces to a linear equation. The temperature-transforming model and the equivalent heat capacity model differ in significant ways. Comparing eq. (6.243) and eq. (6.234), one can conclude that the equivalent heat capacity model is a special case of eq. (6.243), with ∂b/∂t = 0 and cp independent of the spatial variables x, y, z and time. This is Chapter 6 Melting and Solidification 463 the underlying reason why many studies using the equivalent heat capacity model have encountered difficulty in selecting the grid size and time step and have often produced physically unrealistic oscillatory results. In order to satisfy ∂b/∂t = 0 and ∂ (c pT ) /∂t = c p ∂T /∂t , the time step has to be small enough to assure that cp and b are independent of both time and space variables – a difficult criterion to satisfy. Equation (6.243) can be solved by many numerical methods, including the finite volume approach by Patankar (1980). To verify the accuracy of the temperature-transforming model, a numerical solution of the Stefan problem is presented. The physical model of the Stefan problem has been given in the previous section. Figure 6.19 shows a comparison between the exact solution and the numerical solution using the temperature transforming model for the dimensionless location of the melting front as a function of the square root of dimensionless time. The definition of dimensionless melting front location and dimensionless time is eq. (6.34). The temperature-transforming model eliminates the time-step and grid-size limitations and is insensitive to the phase-change temperature range. Therefore, the model can properly handle phase change occurring over a temperature range (as in the case of multicomponent PCM phase change) or at a single temperature (as in the case of single-component PCM phase change). The temperaturetransforming model has been used to solve many different solid-liquid phase change problems (Cao and Faghri 1990a, 1991, 1992; Zeng and Faghri 1994a, 1994b; Zhang and Faghri 1996b). Ste = 0.2 K sA = 1 C sA = 1 S 1/2 Figure 6.19 Comparison of interface location for the Stefan problem obtained by exact and numerical solutions using temperature transforming model. 464 Transport Phenomena in Multiphase Systems 6.6 Solidification of a Binary Solution System 6.6.1 Overview Solidification of multicomponent PCMs has been investigated due to its importance in the fields of metallurgy (Viskanta, 1988), crystal growth (Bardsley et al., 1979), oceanography (Hobbs, 1974), and modeling (Beckermann and Wang, 1995; Bennon and Incropera, 1987a, b; Zeng and Faghri, 1994a, b), as well as many others. It is characterized by the presence of a variety of microscopically complex interfacial structures; as a result, the macroscopic solidliquid interface is not smooth, as is the case in the solidification of singlecomponent PCMs. The most popular microstructure is known as the dendrite, which can exist in either cocolumnar or equiaxed forms (Beckermann and Wang, 1995). From the heat transfer point of view, the microscopic region where solid and liquid coexist is considered a mushy zone that exists between pure solid and pure liquid phase. The temperature at the interface between the mushy zone and the liquid region is the liquidus line temperature, while the interfacial temperature between the solid and the mushy zone is the solidus temperature (see Fig. 2.11). Therefore, phase change of multicomponent PCMs can be considered to occur over a range of temperatures. The presence of a mushy zone in solidliquid phase change system not only characterizes the phase change of multicomponent mixtures, but also poses the primary challenge in analyzing and modeling of binary and multicomponent solid-liquid phase change systems. As indicated in Table 1.10, melting and solidification of multicomponent PCMs can be viewed as dispersed phase. Therefore, the volume-averaged governing equations presented in Chapter 4 can be used to describe the transport phenomena in the mushy zone (Beckermann and Wang, 1995). Some mathematical models have been developed to describe the solidification process of a binary system. In general, such models can be divided into two classes: phase-averaged (multi-fluid model) and overall mass- or volume-averaged (mixture or homogeneous) models. Since the mushy zone can be described as a dispersed phase region, the phase-averaged models (multi-fluid model in Chapter 4) develop a set of conservation equations for each phase, as well as the coupled relationship between the phases (Beckermann and Wang, 1995). The advantage of phase-averaged model is that the nonequilibrium between the two phases can be addressed. The disadvantage is that they require modeling of the interface exchange terms that are related to the microstructure in the mushy zone. The volume-averaged model (homogeneous or mixture model in Chapter 4), on the other hand, reduces the general system of two-phase flow equations to a set of continuum conservation equations that can be applied to the entire solidification domain, including the solid, liquid, and mushy zones (Bennon and Incropera, 1987a, b; Zeng and Faghri, 1994a, b). Using the phase-averaged model, computational studies have identified important features such as liquidus Chapter 6 Melting and Solidification 465 interface irregularity, remelting channeling in the mushy zone, and double diffusive convection in the liquid. An experimental investigation was presented by McDonough and Faghri (1993) concerning the solidification of an aqueous sodium carbonate solution around a vertical cylinder. The liquid-mushy and mushy-solid interface locations are determined by the pulse-echo ultrasonic measurement technique. Ultrasonically measured interface locations were compared with visual observations of the liquid-mushy interface, and temperature measures of the mushy-liquid interface. The ultrasonically measured liquid-mushy interface was shown to be in good agreement with visual observations. High attenuation and scattering of the ultrasonic pulse in the mushy region made the location of the mushy-solid interface difficult to determine ultrasonically. Many researchers conduct experiments using a transparent PCM, such as NH4Cl-H2O solution, because its solidification is quite similar to the solidification of alloys and is easy to observe. An integral approximate solution of solidification of a binary solution on a cold isothermal surface will be presented in Section 6.6.2. Numerical simulation of a convection-controlled solidification problem using a mixture or multi-fluid model will be presented in Sections 6.6.3 and 6.6.4. 6.6.2 Integral Approximate Method The mushy zone formation in binary solutions on a horizontal cold surface will be studied using an integral approximate method. The physical model of the solidification problem is shown in Fig. 6.20 (Zhang and Faghri, 1997). An Figure 6.20 Solidification of a binary solution (Zhang and Faghri, 1997). 466 Transport Phenomena in Multiphase Systems ammonium chloride water solution with initial temperature Ti and mass fraction ωi fills a half-space x > 0 . At time t=0, the wall temperature at x=0 is suddenly reduced to a temperature Tw , which is higher than the eutectic temperature of the ammonium chloride water solution. The solidification process starts from the cold wall and the mushy zone grows upward. There is no solid phase because the cold wall temperature, Tw , is above the eutectic temperature. Therefore, this is a two-region problem with temperatures in the mushy and liquid zones, as well as the location of interface between these two zones, as unknowns. During the solidification process, a denser and colder solution appears in the mushy zone because the solution near the ice surface rejects the solute (Braga and Viskanta, 1990). Thus, the liquid phase is hydrodynamically stable and no natural convection occurs in the liquid phase. Furthermore, the following assumptions are made in order to simplify the analysis: 1. The growth of the mushy zone is controlled by heat conduction because the thermal diffusivity of the salt solution is 100 times greater than the mass diffusivity (Fang et al., 1984). 2. The properties of the solid and liquid are constant within the liquid and solid, but different from each other, and the properties of the mushy zone are weighted according to the local solid fraction. 3. The densities of the solid and liquid phases are the same, i.e., the density change during the solidification process is neglected. Based on the above assumptions, the governing energy equations of the solidification problem can be given as follows: Mushy zone: ∂T · ∂ ∂§ ∂f ρ mu (cmuTmu ) = ¨ kmu mu ¸ + ρ mu hsA 0 ≤ x ≤ s (6.244) ∂t ∂x © ∂x ¹ ∂t T = Tw x=0 (6.245) Liquid zone: ∂TA ∂ 2T = kA 2A x≥s ∂t ∂x TA = Ti x→∞ TA = Ti t =0 Mushy zone and liquid zone interface: Tmu = TA = Ts x=s ∂T ∂T kmu mu = kA A x=s ∂x ∂x ρ A c pA (6.246) (6.247) (6.248) (6.249) (6.250) Chapter 6 Melting and Solidification 467 s=0 t =0 (6.251) where the last term in eq. (6.244) represents the latent heat released during solidification. Since phase change occurs within the entire mushy zone, the latent heat appears as a source term in the energy equation (6.244), instead of appearing in a boundary condition at interface, as was the case for solid-liquid phase change of single-component PCMs. f is the local solid mass fraction in the mushy zone, which is the same as the volume fraction of solid if the density is not changed during the phase change process. It can be determined from the phase diagram by ω (Tmu ) − ω f= (6.252) ω (Tmu ) where ω (Tmu ) is the mass fraction obtained by liquidus 1 in the phase diagram, Fig. 2.11. Ts in eq. (6.249) is a liquidus line temperature distinguishing the mushy zone and the liquid zone. The liquidus line temperature can be determined from the liquidus line equation of the phase diagram of ammonium chloride water solution based on the initial concentration (see Example 2.4). The properties in the mushy zone are mass weighted according to the local solid fraction as follows: ρ mu = ρA = ρ s (6.253) c pmu = (c ps − c pA ) f + c pA (6.254) kmu = (ks − kA ) f + kA By defining the following dimensionless variables, T −T T −T T −T θ mu = mu s θ A = A s θi = i s Ts − Tw Ts − Tw Ts − Tw c pA (Ts − Tw ) c ps αt s S = τ = A2 Ste = Rc = L L hsA c pA (6.255) X= x L ½ ° ° ° ° ¾ ° s° ° ° A¿ k Rk = k (6.256) the governing equations and boundary conditions become ∂ [(( Rc − 1) f + 1)θ mu ] ∂τ ∂θ mu º 1 ∂f ∂ª = 0≤ X ≤S «(( Rk − 1) f + 1) ∂X » + Ste ∂τ ∂X ¬ ¼ θ mu = −1 X =0 ∂θ A ∂ θ A = X >S ∂τ ∂X 2 θA = θi X →∞ θA = θi τ =0 θ mu = θ A = 0 X =S 2 (6.257) (6.258) (6.259) (6.260) (6.261) (6.262) 468 Transport Phenomena in Multiphase Systems ∂θ mu ∂θ A = X =S (6.263) ∂X ∂X (6.264) S =0 τ =0 The exact analytical solution of eq. (6.259) can be obtained by using the exact solution of heat conduction in a semi-infinite body, i.e., ª erfc( X / 4τ ) º (6.265) θ A = θi «1 − » ¬ erfc( S / 4τ ) ¼ which exactly satisfies eqs. (6.259) – (6.262). For solidification on a cold isothermal surface, the thickness of the mushy zone, S, can be expressed as (Braga and Viskanta, 1990; Ozisik, 1993; Tien and Geiger, 1967) S=2λ τ (6.266) where is a constant. Integrating eq. (6.257) over the interval of (0,S) yields the integral equation of the mushy zone: dS [(( Rc − 1) f + 1)θmu ]dX dτ ³0 (6.267) ∂θ mu ∂θ mu 1dS = − [( Rk − 1) f + 1] + fdX ∂X X = S ∂X X = 0 Ste dτ ³0 where ( ∂θ mu / ∂X ) X = S can be found from eq. (6.264) via eq. (6.260). Assuming that the temperature distribution in the mushy zone is a secondorder polynomial function, and determining the coefficients, one obtains the temperature distribution in the mushy zone: 2 2 ª 2λ e − λ θ i º § X − S · ª 2λ e − λ θ i º § X − S · 2 (6.268) +« − 1» ¨ θ mu = « »¨ ¸ ¸ « erfc(λ ) » © S ¹ « erfc(λ ) »© S ¹ ¬ ¼ ¬ ¼ To solve the solidification problem using the integral approximate method, it is necessary to know the distribution of the solid mass fraction f in the mushy zone. The expression of the solid fraction as a single value of temperature in the mushy zone – as in Braga and Viskanta (1990) – is impossible to use here because the integral term in eq. (6.268) will be very difficult to obtain. In their integral solution of the solidification in a semi-infinite region without liquid superheat, Tien and Geiger (1967) assume a linear distribution of the solid fraction in the mushy zone. In that instance, the solid fraction distribution is found to have no significant effect. Cao and Poulikakos (1991) assumed that the solid fraction, and its derivative with respect to temperature, are constants when solving freezing problems of a binary alloy saturating a packed bed of spheres. For the sake of simplicity, it is also assumed here that the solid fraction is a linear function in the mushy zone, i.e., Chapter 6 Melting and Solidification 469 X· § (6.269) f = f w ¨1 − ¸ S¹ © where fw in eq. (6.269) is the solid mass fraction at the cold isothermal surface, which depends on the cold wall temperature Tw and the initial concentration of the ammonium chloride water solution ωi . It can be determined from the phase diagram by the lever rule (see Example 2.4): fw = 1 − ωi ω (Tw ) (6.270) where ω (Tw ) can be determined by the liquidus line equation: ω (Tw ) = 1.678 × 10−3 − 1.602 × 10−2 Tw − 2.857 × 10−4 Tw2 − 4.491 × 10−6 Tw3 (6.271) Substituting eqs. (6.268) and (6.269) into eq. (6.267) and considering eq. (6.266), an algebraic equation of is obtained: ­ ½2 2λ e − λ θ 0 [( Rk − 1) f w + 2] ° ° 2[( Rk − 1) f w + 1] − erfc(λ ) ° ° (6.272) λ=® ¾ −λ 2 ° 2λ e θ 0 f w + 2 + 3( Rc − 1) f w + 4 + f w ° ° erfc(λ ) 12 12 2Ste ° ¯ ¿ which is simply an algebraic equation that can be solved iteratively. After the value of is obtained, the dimensionless thickness of the mushy zone can be obtained from eq. (6.266), and the temperature distribution in the liquid and mushy zones can be obtained by eqs. (6.265) and (6.268). In order to compare these results based on the integral method with Braga and Viskanta’s (1990) experimental results, the following dimensionless variable should be introduced: X η= (6.273) 2τ The temperature distributions in the liquid and mushy zones then become ª erfc(η ) º θ A = θi «1 − (6.274) » ¬ erfc(λ ) ¼ 2 1 ª 2λ e − λ θi º § η − λ · ª 2λ e − λ θ i º § η − λ · 2 (6.275) − 1» ¨ θm = « »¨ ¸+« ¸ « erfc(λ ) » © λ ¹ « erfc(λ ) »© λ ¹ ¬ ¼ ¬ ¼ The thermal properties of the ammonium chloride water solution at specified concentrations of 5%, 10%, and 15% can be found in Zeng and Faghri (1994a) and Cao and Poulikakos (1991), and will not be repeated here. In Braga and Viskanta’s (1990) experimental investigation, the solidification of NH4Cl-H2O was performed for three different cases. Table 6.1 compares the mushy zone thickness obtained using the integral method by Zhang and Faghri (1997) and by Braga and Viskanta (1990). It can be seen that the result obtained by the integral approximate solution agreed well with Braga and Viskanta’s similarity solution and the experimental results. 2 2 470 Transport Phenomena in Multiphase Systems Table 6.1 Comparison of the mushy zone thickness obtained by different methods (Zhang and Faghri, 1997) s (mm) ωi (%) t (min) Ti ( DC ) Tw ( DC ) Braga and Viskanta (1990) Eq. (6.266) Analytical Experimental 5 15.3 -13.3 500 51.0 48.0 50.7 10 20.2 -14.6 540 36.6 32.0 33.6 15 19.5 -14.7 540 16.3 13.0 14.8 6.6.3 Mixture Model The purpose of this subsection is to develop and improve the temperature transforming model so as to deal more effectively with binary solid-liquid phasechange problems, with emphasis on the inclusion of the effect of the mushy zone in the energy equation. The essential feature of the model proposed by Zeng and Faghri (1994a, b) is its separation of the coupled effects of temperature and concentration on the latent heat evolution within the mushy zone. The assumptions made by Zeng and Faghri (1994a) using the mixture model (mass-averaged) include: (1) the thermophysical properties of two phases are constant, although not necessarily equal, and (2) the solid and liquid in the mushy zone are assumed to be in thermal equilibrium. The governing equations will be given for solidification of a binary solution in a two-dimensional rectangular cavity cooled from the sidewall. The mixture continuity equation is Dρ + ρ∇ ⋅ V = 0 (6.276) Dt where the mixture density is defined as ρ = ε s ρ s + ε A ρA (6.277) and the volume fractions of solid and liquid ε s and ε A satisfy εs + εA = 1 (6.278) The mass-averaged velocity in eq. (6.276) is defined as V = f s Vs + f A VA (6.279) where Vs and VA are velocities of solid and liquid phases, respectively; the mass fractions of solid and liquid satisfy f s + fA = 1 (6.280) The volume fraction ε and the mass fraction f are related by ρ sε s ε s ρ s + ε A ρA ρ Aε A fA = ε s ρ s + ε A ρA fs = (6.281) (6.282) The momentum equations in x- and y-directions are Chapter 6 Melting and Solidification 471 where us and vs are the components of the solid velocity, Vs, in the x- and ydirections, respectively. X and Y are components of body force in the x- and ydirections. The energy equation is ∂ ( ρ h ) + ∇ ⋅ ( ρ Vh ) = ∇ ⋅ (k ∇T ) − ∇ ⋅ [ ρ ( hA − h )(V − Vs )] (6.285) ∂t where the mixture enthalpy is defined as h = f s hs + f A hA (6.286) The mass species equation is ∂ ( ρω ) + ∇ ⋅ ( ρ Vω ) = ∇ ⋅ ( ρ D∇ω ) (6.287) ∂t +∇ ⋅ [ ρ D∇(ωA − ω )] − ∇ ⋅ [ ρ (ωA − ω )(V − Vs )] where the species mass fraction is defined as ω = f s ω s + f A ωA (6.288) The permeability of the two-phase mushy zone, which is modeled as a porous medium, is K = K0 § ·μ ∂( ρu) ∂p ρ + ∇ ⋅ ( ρ Vu ) = − + ∇ ⋅ ¨ μA ∇u ¸ − A ∂t ∂x © ρA ¹K § ·μ ρ ∂( ρ v) ∂p + ∇ ⋅ ( ρ Vv ) = − + ∇ ⋅ ¨ μA ∇v ¸ − A ∂t ∂y © ρA ¹K ρ (u − us ) + X ρA ρ (v − vs ) + Y ρA (6.283) (6.284) ε A3 (1 − ε A ) 2 (6.289) where K0 is the permeability constant. The thermal conductivity of the mixture is k = ε s k s + ε A kA (6.290) Under constant specific heat assumption, the enthalpy and temperature are related by the specific heats as hs = ³ c ps dT = c psT 0 T (6.291) (6.292) hA = ³ c pA dT + h0 = c pAT + h0 0 T where h0 is the reference enthalpy for the liquid phase. The mass diffusivity, in light of the assumption neglecting diffusion in the solid phase, is D = f A DA (6.293) Substituting eqs. (6.291) and (6.292) into eq. (6.286), the enthalpy becomes h = c pT + f A h0 (6.294) where c p = f s c p s + fA c p A (6.295) 472 Transport Phenomena in Multiphase Systems Substituting eq. (6.294) into eq. (6.285), the energy equation becomes (Zeng and Faghri, 1994a) ∂T A(T , ω ) + ∇ ⋅ ( ρ Vc pAT ) = ∇ ⋅ (k ∇T ) − B (T , ω ) (6.296) ∂t where A(T , ω ) is an effective heat coefficient, ∂f A(T , ω ) = ρ c p + a (T , ω ) A (6.297) ∂T and B (T , ω ) is a heat source term ∂f ∂ω B (T , ω ) = −∇ ⋅ ( ρ h0 V ) + ∇ ⋅ { ρ f s [h0 + (c pA − c ps )T ]Vs } − a (T , ω ) A (6.298) ∂ω ∂t where ( ρ A − ρ s )c pT + ρA h0 º ª (6.299) a (T , ω ) = ρ « (c pA − c ps )T + » ρA − f A ( ρA − ρ s ) ¼ ¬ The temperature-transforming model presented by Zeng and Faghri (1994a) can be considered as a combination of the apparent heat capacity method and the latent heat source term method designed to deal with the latent heat evolution in binary phase-change systems. The latent heat evolution with reference to temperature evolution is accounted for in the energy equation by defining the effective heat coefficient A(T , ω ) , which is, in turn, similar to the apparent heat capacity ρ dh / dT , while the latent heat evolution – with reference to concentration variation – is incorporated into the source term B (T , ω ) , which is similar to the source-based method. The model develops and expands the fixedgrid method to a binary phase-change system and links the source-based and the apparent heat capacity methods. For the case of ∂fA / ∂ω = 0 , that is, A(T , ω ) = A(T ) (6.300) B (T , ω ) = B (T ) (6.301) eq. (6.296) reduces to the apparent heat capacity formulation. For the case of ∂f A / ∂T = 0 , A(T , ω ) = ρ c p (6.302) B (T , ω ) = B (T , ω ) (6.303) eq. (6.296) represents the so-called source-based energy equation. It should be pointed out that, in an iterative solution of the discrete equation, the extra term a (T , ω )(∂f A / ∂T ) in the effective heat coefficient A(T , ω ) , which is positive, will play the role of a damping factor, in that the temperature of the nodes in which phase change is occurring will change more slowly from one iteration to the next. This stability may represent a main advantage of the apparent heat capacity method. A disadvantage, however, is that if the liquid fraction curve f A (T , ω ) is very steep or contains some discontinuities, instabilities may occur in the solution procedure. This tendency is reduced in the present formulation, because the binary phase-change system generally releases or Chapter 6 Melting and Solidification 473 absorbs latent heat over a relatively wide temperature range. Furthermore, the energy equation has been transformed into a nonlinear equation with a single dependent variable T. Over the iterative loop that solves the energy equation, the extra term −a(T , ω )(∂f A / ∂ω )(∂ω / ∂T ) in the source term is seen only as a function of temperature. Zeng and Faghri (1994b) solved the solidification of an NH4Cl solution in a two-dimensional rectangular cavity of dimensions 36×144 mm2 (width × height). Initially, the cavity is insulated and charged with a superheated solution at temperature Ti and concentration ωi = 30% of NH4Cl by weight. Solidification is initiated at t = 0 by holding the left vertical boundary to a temperature Tc < Te, while maintaining the right vertical boundary at the initial temperature. Since the initial concentration of the solution was greater than the eutectic concentration ωe the phase change occurred within mushy zone 2, as shown in Figure 2.11. It is further assumed that: (1) the densities of the two phases are equal and constant, and (2) the velocity of the solid phase is zero. The components of body force in the x- and y- directions are: (6.304) X =0 Y = ρ g[ βT (T − Tref ) + β s (ωref − ω )] (6.305) where βT and β s are thermal and solutal expansion coefficients, respectively. With the additional assumptions, a (T , ω ) and B (T , ω ) are reduced to a(T , ω ) = ρ ª(c pA − c ps )T + hsA º ¬ ¼ B (T , ω ) = − a (T , ω ) (6.306) (6.307) ∂f A ∂ω ∂ω ∂t where Tm − TA ∂f A = ∂T (1 − k p )(Tm − T ) 2 Tm − Te ∂f A = ∂ω (1 − k p )(Tm − T )ωe (6.308) (6.309) where Tm is the fusion temperature for ω = 1 , and k p is the equilibrium partition ratio. The above governing equations are solved using a finite volume method. A 40×40 grid, in which the space step in the y-direction remained constant, while the variation in the x-direction was expressed as a symmetric arithmetic progression, was used to compute this problem with a time step of 3 s. Iterations performed within a time step were terminated when the residual mass within each control volume fell below 10-5, and when the following values fell below 0.001%: (1) changes in the minimum and maximum concentration within each control volume, (2) the average cold wall heat extraction, and (3) total liquid fraction. Figure 6.21 shows temperature, concentration, solid fraction, and stream lines of solidification in the rectangular cavity for t = 8 min. Remelting within 474 Transport Phenomena in Multiphase Systems Figure 6.21 Contours for solidification of NH4Cl-H2O binary solution in a rectangular cavity at t = 8 min: (a) temperature, (b) concentration, (c) liquid fraction, (d) streamlines (Zeng and Faghri, 1994b; Reproduced by permission of Routledge/Taylor Francis Group, LLC). the mushy region adjoining the top layer is evident, and channels have begun to form [Fig. 6.21(c)]. Outflow of cold, water-rich solution from the mushy region is pronounced at the top of the cavity. Conditions within the bulk liquid have clearly become double-diffusive in the sense that a stable density interface separates the colder, water-rich solution in an underlying layer, as shown in Fig. 6.21(d). 6.6.4 Volume-Averaging Model Since both liquid and solid are present in the mushy zone, the governing equations in the mushy zone can be obtained by performing volume averaging to each phase to yield a volume-averaging model. The volume-averaged multi-fluid continuity equation (4.22), momentum equation (4.34), energy equation (4.44), and species equation (4.54) introduced in Section 4.3 are also applicable to alloy k solidification. Considering the intrinsic phase-averaged density ρ k is the same as the density ρ k , assuming that the liquid and solid phases are incompressible, and there are no internal heat generations in either phases, and neglecting viscous dissipation, the volume averaged governing equations become Chapter 6 Melting and Solidification 475 ∂ (ε k ρk ) + ∇ ⋅ ε k ρ k Vk ∂t ∂ ε k ρ k Vk ∂t = ∇ ⋅ εk ∂ ε k ρ k hk ∂t ∂ ε k ρ k ωk ,i ∂t ( k )= k jk j =1( j ≠ k ) ¦ Π m′′′ jk (6.310) ( k ( k k ) + ∇ ⋅ (ε ρ V V ) )+ε ρ X + ¦ ( F k k k k Π k k k j =1( j ≠ k ) k + m′′′ jk Vk , I k ) (6.311) ( k ) + ∇ ⋅ ε k ρ k Vk hk ( ) = −∇ ⋅ q′′ + k ª q′′′ + m′′′ h jk k ,I « jk ¬ j =1( j ≠ k ) ¦ Π k º » ¼ ( k ) + ∇ ⋅ ε k ρ k ωk ,i ( k Vk k ) (6.312) ′′′ = −∇ ⋅ m′′,i + mk ,i + k k k j =1( j ≠ k ) ¦ Π m′′′ ,i jk (6.313) where k can be either (s) or liquid ( A ), Vk , I th and hk , I are the intrinsic phase-averaged velocity and enthalpy of the k phase at the interface. The above equations are generalized and applicable to any two-phase system. In the alloy solidification process, the solid phase can be assumed to be rigid and stationary as in purely columnar growth (Beckermann and Viskanta, 1993), i.e., s Vs = 0 . The continuity equations in the solid and liquid phases become ∂ (ε s ρ s ) = mA′′′ s ∂t (6.314) A ∂ (ε A ρA ) + ∇ ⋅ ε A ρA VA ∂t ( ) = m′′′ As (6.315) where (6.316) ′′′ ′′′ mAs + msA = 0 (6.317) The momentum equation for the solid phase is not needed, and the momentum equation for the liquid phase is ∂ A A A ε A ρA VA + ∇ ⋅ ε A ρA VA VA ∂t (6.318) εs + εA = 1 ( = ∇ ⋅ εA ( A A )( )+ε ρ g + F A A ) sA ′′′ + msA VA , I A ˆˆ − ∇ ⋅ ε A ρA VA VA A ( A where VA VA A in eq. (6.311) is replaced by VA VA A ) ˆˆ + VA VA and XA is replaced by g, i.e., gravity is the only body force. Since the liquid velocity in the mushy zone is usually laminar because of the characteristic diameter, the product A ˆˆ of the velocity deviations, V V , is also small and can be neglected. V is A A A,I the intrinsic phase-averaged velocity of the k phase at the interface. The difference between two adjacent phases results solely from the density difference between the two phases. Since the density change is usually very small during th 476 Transport Phenomena in Multiphase Systems solid-liquid phase change ( ρA ≈ ρ s ), the intrinsic phase-averaged velocities for both phases at the solid-liquid liquid interface are approximately equal to each other, i.e., VA , I A ≈ Vs , I s . The solid phase is rigid and stationary, thus the s velocity of the solid phase at the interface, Vs , I = 0 . The fourth term on the right hand side of eq. (6.318) can therefore be neglected. The momentum equation for the liquid phase becomes ∂ A A A A ε A ρA VA + ∇ ⋅ ε A ρA VA VA = ∇ ⋅ ε A A + ε A ρA g + FsA (6.319) ∂t Applying the stress-strain relationship for Newton’s fluid for the first term on the right-hand side and replacing the third term by the Darcy-Brinkman’s relation (see Section 4.6.2), the momentum equation becomes ª 1 ∂ VA A º VA + 2 ∇ ⋅ VA » = −∇p + ρA « ε « ε ∂t » ¬ ¼ (6.320) μA 2 μA ρ A CA ∇ VA + ρA g − VA − 1 2 VA VA ε K K A where VA = ε A VA is extrinsic phase-averaged velocity of the liquid phase. The energy equations for the solid and liquid phases are s ∂ s ′′′ ε s ρ s hs = −∇ ⋅ q′′ + qs′′′ + mAs hs , I (6.321) s A ∂t A ∂ A A ε A ρA hA + ∇ ⋅ ε A ρA VA hA = −∇ ⋅ q′′ + qs′′′ + ms′′′ hA , I (6.322) A A A ∂t A s Adding eqs. (6.321) and (6.322) together, considering hA − hs = hsA , and assuming the solid and liquid phases are in thermal equilibrium, s A ( Ts = TA = T ), the energy equation becomes ( ) ( ) ( ) ( ) ( ) ( ) ∂ε ∂T A + ε A ( ρ c p )A VA ⋅ ∇T = ∇ ⋅ ( k ∇T ) + ρ s hsA s (6.323) ∂t ∂t where eqs. (6.314) and (6.317) are used to determine the mass production rate of the solid and liquid phases per unit volume. The volume-averaged heat capacity and thermal conductivity are defined by ρ c p = ε A ( ρ c p )A + ε s ( ρ c p ) s (6.324) ρcp k = ε A kA + ε s k s (6.325) th Assuming there is no source/sink of the i component in the solid and liquid phases, the conservations of species of the ith component for the solid and liquid phases are s ∂ ′′′ ε s ρ s ωs ,i = −∇ ⋅ m′′,i + mAs ,i (6.326) s ∂t ( ) Chapter 6 Melting and Solidification 477 A A ∂ A ε A ρA ωA ,i + ∇ ⋅ ε A ρA ωA ,i VA = −∇ ⋅ m′′,i + ms′′′,i (6.327) A A ∂t ′′′ ′′′ Adding eqs. (6.326) and (6.327) together and considering mAs ,i + msA ,i = 0 , the conservation of species for the ith component becomes A s A ∂ A ε A ρA ωA ,i + ε s ρ s ωs ,i + ∇ ⋅ ε A ρA ωA ,i VA = −∇ ⋅ m′′,i + m′′,i k s ∂t (6.328) ( ) ( ) ( ) ( ) ( ) CmixAcet height (m) height (m) height (m) distance from heat exchanger (m) distance from heat exchanger (m) distance from heat exchanger (m) (a) (b) (c) height (m) distance from heat exchanger (m) (d) Figure 6.22 Solidification of SNC-acetone solution (t = 20246 s, before deformation): (a) mesh, (b) temperature field, (c) solid fraction, and (d) mixture concentration (Rizzo et al., 2003). CmixAcet height (m) height (m) height (m) distance from heat exchanger (m) distance from heat exchanger (m) distance from heat exchanger (m) (a) (b) (c) height (m) distance from heat exchanger (m) (d) Figure 6.23 Solidification of SNC-acetone solution (t = 26301 s, after deformation): (a) mesh, (b) temperature field, (c) solid fraction, and (d) mixture concentration (Rizzo et al., 2003). 478 Transport Phenomena in Multiphase Systems Rizzo et al. (2003) presented a numerical simulation of solidification and mushy zone deformation of an organic alloy (SNC-acetone) in a two-dimensional rectangular cavity cooled from its left sidewall and all other three walls were kept adiabatic. Succinonitrile (SCN) is a transparent organic substance in the liquid state, and acetone was added to obtain an alloy behavior. The liquidus temperature for the pure SCN is 58.1 °C. The initial concentration of acetone was 11.3 wt% and the initial temperature of the melt alloy was at 42.35 °C. The cold wall temperature was set to 5 °C in numerical simulation. The rectangular cavity was deformed from its initial size of 50×170 mm to a final size of 40×210 mm. Figure 6.22 shows the results at the time relative to the deformation beginning, i.e., at 26246 s from cooling start, while Fig. 6.23 shows results at the end of deformation (26301 s). 6.7 Contact Melting 6.7.1 Fixed Melting and Contact Melting Two types of melting patterns may be observed during the experimental investigation of melting in a two-dimensional rectangular cavity: fixed melting (Fig. 6.24) and contact melting (Fig. 6.25). The solid core is fixed during the fixed melting process, and heat transfer is controlled mainly by natural convection. During the contact melting process, on the other hand, gravity causes the solid core to fall so that it maintains constant contact with the bottom surface of the cavity. Heat transfer in the contact melting process is controlled primarily by conduction between the heated bottom surface and the solid PCM. Since the unmelted solid core always maintains contact with the bottom surface of the cavity, the melting proceeds more quickly in contact melting than it does in fixed melting. A complete review of contact melting in various geometries and their applications can be found in Bejan (1994, 2004). Figure 6.24 Fixed melting in a rectangular cavity heated from two sides and bottom (t1<t2<t3). Chapter 6 Melting and Solidification 479 Figure 6.25 Contact melting in a rectangular cavity heated from two sides and bottom (t1<t2<t3). 6.7.2 Contact Melting in a Rectangular Cavity The contact melting in a rectangular cavity shown in Fig. 6.26 was modeled by Hirata et al. (1991). The initial temperature of the solid is assumed to be at the melting point of the PCM, Tm At time t = 0 , the temperatures of the sidewall and bottom wall suddenly increase to Tw, which is higher than Tm. Melting occurs at the sidewalls and bottom wall simultaneously, and the solid falls due to gravity, always keeping in good contact with the bottom surface of the cavity. As a result, the liquid layer between the solid and the bottom surface of the cavity is very thin. The liquid layer thickness in Fig. 6.26 is exaggerated for clarity of presentation. The following assumptions are needed in order to solve the problem: 1. Compared with melting at the bottom wall, melting at the two sides and above the solid can be neglected. 2. The liquid is Newtonian. 3. Compared with the viscous term, the inertia term of the momentum equation is negligible. 4. Compared with heat conduction, heat transfer by convection is negligible. 5. The thermophysical properties of the liquid are constant. 6. The flow of the liquid layer is symmetric at x = 0. Nusselt liquid film theory is applicable to this melting process. The momentum and energy equations for the liquid film are ∂ 2u dp μA 2 = (6.329) dx ∂y ∂ 2T =0 ∂y 2 (6.330) 480 Transport Phenomena in Multiphase Systems Figure 6.26 Physical model of contact melting. The corresponding boundary conditions are y = 0: u = 0, T = Tw (6.331) y =δ : u = 0, T = Tm (6.332) Integrating eqs. (6.329) and (6.330) and considering boundary conditions (6.331) and (6.332), the velocity and temperature distribution in the liquid layer are 1 dp 2 ( y − δ y) u= (6.333) 2 μA dx y T = Tw + (Tm − Tw ) (6.334) δ The heat balance equation over the range (0, x) in the liquid layer can be given by x δ ∂T dx = ρ hsA ³ udy (6.335) − ³ kA 0 0 ∂y y = 0 The transport of sensible heat has been neglected because the Stefan number is assumed to be much less than unity. Substituting eqs. (6.333) and (6.334) into eq. (6.335) yields x dx ρ h dp 3 −kA (Tm − Tw ) ³ = − A sA δ (6.336) 0δ 12 μA dx Experimental observation (Dong et al., 1991) indicates that the variation of the liquid layer thickness with respect to x is negligible, and therefore δ can be treated as a constant across x. Chapter 6 Melting and Solidification 481 0 δ where m′′ is the molten mass per unit area and unit time, which is independent of x . Integrating eq. (6.337) with respect to x, one obtains the liquid layer thickness as k (T − T ) δ= A w m (6.338) m′′ hsA Substituting eq. (6.338) into eq. (6.336) and integrating it with boundary condition p = 0 at x = W / 2 yields the pressure distribution in the liquid layer: 0 Neglecting the effect of sensible heat in the liquid layer, the heat balance equation for the range of (0, x) can be given as x dx x − kA (Tm − Tw ) ³ = ³ hsA m′′ dx (6.337) 2 ª º ª 2 §W · º hsA p = −6m′′ ν A « » «x − ¨ ¸ » ©2¹» ¬ kA (Tw − Tm ) ¼ « ¬ ¼ The force balance of the buoyancy and pressure for the solid is 4 3 (6.339) 2³ 4 W/2 0 pdx = g ( ρ s − ρA )W ( H − s ) 3 (6.340) Substituting eq. (6.339) into eq. (6.340) yields ª hsAW º (6.341) m′′ ν A « » = g ( ρ s − ρ A )W ( H − s ) ¬ kA (Tw − Tm ) ¼ The relationship between m′′ and s is (see Fig. 6.23) ds = ρ sVs m′′ = ρ s (6.342) dt where Vs is the downward velocity of the solid. Substituting eq. (6.342) into (6.341), the following ordinary differential equation for s is obtained: º § ds · ª hsAW » = g ( ρ s − ρA )W ( H − s ) ¨ ρs ¸ ν A « © dt ¹ ¬ kA (Tw − Tm ) ¼ The downward velocity of the solid is then 4 3 1 3 (6.343) ª k (T − T ) º 4 ª g ( ρ s − ρA )( H − s ) º 4 (6.344) Vs = « A w m » « » ρ sW 2ν A ¬ ρ s hsA ¼ ¬ ¼ The order of magnitude of the downward velocity of the solid was obtained by scale analysis in Chapter 1 and the result was [see eq. (1.208)] § k ΔT · 4 § p · 4 Vs  ¨ A (6.345) ¸ ¨ 2¸ © ρ s hsA ¹ © μ L ¹ We will show below that the result of analytical solution, eq. (6.344), is consistent with the results of the scale analysis, eq. (6.345). It follows from the 3 1 482 Transport Phenomena in Multiphase Systems force balance of the solid, eq. (6.340), that the order of magnitude of the pressure in the liquid layer is p  g ( ρ s − ρA )( H − s ) (6.346) The order of magnitude of the second bracket on the right-hand side of eq. (6.344) is therefore ª g ( ρ s − ρA )( H − s ) º 4 ª p º 4 « » « » 2 ρ sW 2ν A ¬ ¼ ¬ ρ sW ν A ¼ Since ρ s  ρA and W  L , eq. (6.347) becomes 1 1 1 1 (6.347) ª g ( ρ s − ρ A )( H − s ) º 4 ª p º 4 (6.348) « » « 2» ρ sW 2ν A ¬ μA L ¼ ¬ ¼ Substituting eq. (6.348) into eq. (6.344) and considering ΔT = Tw − Tm , eq. (6.345) can be obtained. In general, it is important to find the time it takes to completely melt the solid. Introducing the following dimensionless variables ρ − ρ A gW 3 ρ s H δ½ S= B= Rρ = A Δ = ° Ar = s 2 W W W° ρA νA ρs (6.349) ¾ c pA (Tw − Tm ) α At sS ° Ste = τ = Ste 2 A(τ ) = = ° hsA W HB ¿ where Ar is the Archimedes number and A(τ ) is the melting rate, eq. (6.343) becomes dτ § Ste · 4 § S · 4 =¨ (6.350) ¸ ¨1 − ¸ dS © Ar Pr B ¹ © B ¹ Substituting eq. (6.342) into eq. (6.338) yields k (T − T ) dt δ= A w m (6.351) ρ s hsA ds Substituting eq. (6.349) into eq. (6.351), the dimensionless liquid film thickness becomes δ dτ Δ= = (6.352) W dS Thus, the dimensionless liquid layer thickness can be obtained by substituting eq. (6.350) into eq. (6.352): ª º4 Ste Δ=« » ¬ Ar Pr B[1 − A(τ )] ¼ Equation (6.350) can be rewritten as −1 § Ste · 4 4 dτ = ¨ ¸ [1 − A(τ ) ] BdA(τ ) © Ar Pr B ¹ 1 1 1 −1 (6.353) (6.354) Chapter 6 Melting and Solidification 483 Integrating eq. (6.354) with respect to τ , and considering the initial condition of A(0) = 0 , gives 3 4 § Ste · 4 4 τ= ¨ (6.355) ¸ B 1 − [1 − A(τ ) ] 3 © Ar Pr B ¹ Thus, the melting rate at any time can be obtained by solving eq. (6.355), and then the dimensionless liquid layer thickness at any time can be obtained by using eq. (6.353). When all of the solid has melted, the melting rate A(τ ) = s / H = 1 . The corresponding dimensionless time is 1 { } 4 § Ste · 4 3 4 (6.356) τf = ¨ ¸B 3 © ArPr ¹ which gives the time it takes to completely melt the solid; it can be used to estimate the working time of the latent heat thermal energy storage system with contact melting. 1 6.8 Melting and Solidification in Porous Media Melting and solidification in porous media can find application in a wide range of problems, including freezing and thawing of soil in a cold region, latent heat thermal energy storage, and selective laser sintering of metal powders. Applications in life science include cryopreservation of tissue, human embryos, and organs. While the early studies about phase change in porous media treated melting and solidification as conduction-controlled, the effect of natural convection is often very important. Convection-controlled melting and solidification of PCMs saturated in porous media will be discussed in this section. 6.8.1 Convection-Controlled One-Region Melting Problem Melting in an enclosure filled with a porous medium and solid PCM, studied by Jany and Bejan (1987) will be presented here (see Fig. 6.27). It is assumed that all of the walls of the two-dimensional rectangular cavity are impermeable, and that the initial temperature of the porous medium and PCM, Ti, is at the melting point of the PCM, Tm. At time t = 0, the temperature of the left side of the cavity is increased to Tw , while the right side, top, and bottom of the cavity remain insulated. The following assumptions are made: 1. 2. 3. 4. Liquid flow is very slow, so Darcy’s law is applicable. The PCM and porous medium are in thermal equilibrium. The porous medium is isotropic. Constant thermal properties for both PCM and porous medium. 484 Transport Phenomena in Multiphase Systems Figure 6.27 Convection-controlled melting in a porous medium. Since the temperature in the solid region is equal to the melting point, this is a one-region melting problem. The governing equations for the liquid region are ∂ uA ∂ vA + =0 (6.357) ∂x ∂y where uA and vA are the components of extrinsic phase-averaged velocity in the x- and y-directions [see eq. (4.7)] and they can be obtained by Darcy’s law: A K ∂ pA uA = − (6.358) μA ∂x A ½ K ­ ∂ pA ° ° + ρA g ª1 − β ( T − Tm ) º ¾ ® ¬ ¼ μA ° ∂y ° ¯ ¿ The energy equation is § ∂2 T ∂T ∂T ∂T ∂2 T ϒ + uA + vA = αm ¨ + ¨ ∂x 2 ∂t ∂x ∂y ∂y 2 © vA = − (6.359) · ¸ (6.360) ¸ ¹ where ϒ = ε + (1 − ε ) ( ρ c p ) sm ( ρ c p )A (6.361) is the ratio of the heat capacity of the porous medium saturated with PCM over the heat capacity of the PCM. is the volume fraction of PCM, ( ρ c p )A and ( ρ c p ) sm are heat capacities of the liquid PCM and porous matrix. The thermal diffusivity of the porous medium saturated with the PCM is α m = keff /( ρ c p )A . Chapter 6 Melting and Solidification 485 Equations (6.357) – (6.360) are subjected to the following boundary conditions: uA = 0 , T = Tw , x = 0 (6.362) uA = 0 , vA = 0 , T = Tm , x = s (6.363) = 0, y = 0, H (6.364) ∂y The energy balance at the solid-liquid interface is (Jany and Bejan, 1987) α m c pA § ∂ T ∂s ∂s ∂ T · (6.365) =− − ¨ ¸, x = s ∂t ∂y ∂y ¹ hsA © ∂x where the second term in the parentheses at the right-hand side represents the effect of the inclination of the interface. To change uA and vA into one dependent function, a stream function, ψ , is defined as ∂ψ ∂ψ uA = , vA = − (6.366) ∂y ∂x Introducing the following dimensionless variables uH vH αt x y s X = , Y = , Fo = m2 , S = , U A = A , VA = A αm αm H H H H ∂T θ= T − Tm Tw − Tm , Ψ= c pA (Tw − Tm ) g β KH (Tw − Tm ) ψ , Ra = , Ste = αm ν Aα m hsA (6.367) where Ra and Ste are the Rayleigh and Stefan numbers, respectively, the governing equations and boundary conditions become ∂2Ψ ∂2Ψ ∂θ + 2 = − Ra (6.368) 2 ∂x ∂y ∂X ∂θ ∂θ ∂θ ∂ 2θ ∂ 2θ +U +V = + (6.369) ∂Fo ∂X ∂Y ∂x 2 ∂y 2 ∂Ψ U= (6.370) = 0 , θ = 1, X = 0 ∂Y ∂Ψ U= = 0 , θ = 0, X = S (6.371) ∂Y ∂Ψ ∂θ V =− = 0, = 0, Y = 0,1 (6.372) ∂X ∂Y ∂S § ∂θ ∂S ∂θ · (6.373) = − Ste ¨ − ¸, x = S ∂Fo © ∂X ∂Y ∂Y ¹ for phase-averaging has been dropped in the dimensionless variables where for ease of notation. ϒ 486 Transport Phenomena in Multiphase Systems (a) Ra=12.5 Figure 6.28 Propagation of the melting fronts (Jany and Bejan, 1987). (b) Ra=800 Figure 6.29 Nusselt number for melting in a square cavity (Jany and Bejan, 1988; Reprinted with permission from Elsevier). The dimensionless governing equations were solved numerically using a finite difference method by Jany and Bejan (1987). Propagations of the melting fronts in a square cavity ( L / H = 1 ) for different Rayleigh numbers are shown in Fig. 6.28. It can be seen that both the melting velocity and the shape of the melting front are strongly affected by natural convection. Heat transfer from the heating wall is measured by overall Nusselt number, defined as H ∂T 1 § ∂θ · hH 1 Nu = =− (6.374) ³0 ∂x dy = −³0 ¨ ∂X ¸ X =0 dY keff (Tw − Tm ) © ¹ Figure 6.29 shows Nusselt number as a function of dimensionless time, SteFo, at different Rayleigh numbers for a square cavity. The result for Ra = 0 represents Nusselt number for pure conduction solution. It can be seen that the effect of natural convection is pronounced when the Rayleigh number is greater than 50. For all cases where Ra > 0 , there is a knee point beyond which the Chapter 6 Melting and Solidification 487 Nusselt number drops significantly. The knee point represents the time at which the top of the melting front arrives at the right wall (X = 1) of the cavity. The Nusselt number decreases sharply after the knee point because the height of the PCM being molten decreases after the top of the melting front reaches to the right wall. 6.8.2 An Enthalpy Model for Two-Region Melting/Solidification A more generalized model for solid-liquid phase change in porous media was proposed by Beckermann and Viskanta (1988). Instead of using Darcy’s law, which is valid only at very low velocity, the effect of inertia terms was included in the model. This treatment allows for nonslip conditions at the wall of the cavity, e.g., v = 0 at x = 0 (under Darcy’s law, slip occurs at the wall). The generalized governing equations that are applicable to the entire computational domain, including both liquid and solid phases, are obtained by using an enthalpy model that can be applied to a two-region problem (melting of a subcooled solid or solidification of a superheated liquid). Beckermann and Viskanta (1988) developed their model based on volumeaveraged governing equations (see Section 4.6.6). Figure 6.30 shows the illustration of a phase change process in a porous medium. The volume fraction of pore space in the porous medium is defined as ε p = ΔV p / ΔV (6.375) where Vp is the volume of porous space and V is the total volume. The porous space can be occupied by either liquid ( A ) or solid (s) phase of the PCM. The liquid volume fraction in the porous space is defined as γ = ΔVA / ΔV p (6.376) The volume fraction of the liquid in the porous media is therefore ε A = ΔVA / ΔV = ε pγ (6.377) Porous matrix Solid Liquid V Figure 6.30 Control volume for volume averaging. 488 Transport Phenomena in Multiphase Systems Although the PCM under consideration is single-component and its melting point is well defined, the melting front can have a finite thickness because the phase change can be simultaneously inhabited by solid and liquid in the pore. A two-phase region, similar to a mushy zone for a binary system, exists between the solid and liquid phases. The liquid fraction, γ , in this two phase region varies between 0 to 1, while the average temperature in this two phase region is equal to the melting point, Tm. It was assumed that the liquid fraction and the temperature in the two-phase region were related by T − Tm + ΔT γ= , Tm − ΔT ≤ T ≤ Tm + ΔT (6.378) 2ΔT Viskanta (1988) assumed that: (1) the flow and heat transfer is twodimensional and laminar, (2) the properties of both porous matrix and PCM are homogeneous and isotropic, (3) the porous matrix and the PCM are at thermal equilibrium, (4) the velocities of the porous matrix and solid phase are zero, (5) the flow is incompressible and the Boussinesq approximation (density change is considered only in the buoyancy term) applies, (6) the thermal physical properties are constants, (7) the dispersion flux due to velocity fluctuation is negligible, and (8) the densities of the liquid and solid phases are identical. The continuity, momentum, and energy equations are ∇⋅V = 0 (6.379) ª ∂V 1 º ρA μ + ( V ⋅ ∇)V » = −∇p + A ∇ 2 V « ε A ¬ ∂t ε A εA ¼ (6.380) § μA ρ A C f · − ¨ + 1 2 V ¸ V − ρA gβ (T − Tref ) K ©K ¹ ( ρ c p )eff ∂T ∂γ (6.381) + ρ A c pA (V ⋅ ∇T ) = −∇ ⋅ (keff ∇T ) − ε p ρ A hsA ∂t ∂t A where V = VA is intrinsic phase-averaged velocity. For the porous medium consisting of spherical beads with diameter dp, the permeability can be obtained using the following Kozeny-Carman relation: 2 d pε A3 K= (6.382) 175(1 − ε A ) 2 which indicates that the permeability is zero in the solid region. The Forchheimer inertia coefficient Cf is taken as 0.55 (Ward, 1964). The effective heat capacity ( ρ c p )eff in eq. (6.382) is ( ρ c p )eff = ε p ρ ªγ c pA + (1 − γ )c ps º + (1 − ε p )( ρ c p ) sm ¬ ¼ (6.383) where ( cp)sm is heat capacity of the porous matrix. The effective thermal conductivity of the porous medium saturated with PCM is estimated by (Veinberg, 1967) −k 1/ k (6.384) keff + ε keff3 sm 1/ 3 As − k sm = 0 kAs Chapter 6 Melting and Solidification 489 where kAs is the average thermal conductivity of the PCM: kAs = γ kA + (1 − γ )ks (6.385) Beckermann and Viskanta (1988) applied the above model to solve solidliquid phase change problems as shown in Fig. 6.31. A rectangular cavity is filled with a rigid porous matrix that is saturated with PCM. The left wall temperature is held at Th, which is above the melting point of the PCM, Tm, while the right wall is kept at a temperature Tc, below the melting point. The top and bottom of the cavity are adiabatic. The initial temperature of the porous medium and PCM, Ti, is equal to Tc for a melting problem or Th for a solidification problem. The governing equations were nondimensionalized and numerical simulations for both melting and solidification were performed by Beckermann and Viskanta (1988). They also performed experiments on melting and solidification in a system containing spherical glass beads ( d p = 6 mm , p = 0.385) and gallium as a PCM (Tm = 29.78 °C) in a 3.81×4.76cm2 (L × H) rectangular cavity. Figure 6.32 shows comparisons between numerical and experimental results for melting and solidification. It can be seen that the agreement between numerical and experimental results is very good for both melting and solidification. Under the conditions studied by Beckermann and Viskanta (1988), the effect of natural convection on the melting and solidification was insignificant because the interfaces shown in Fig. 6.32 are nearly planar and Adiabatic H Two-phase region Liquid Th Tc Solid y 0 x L Adiabatic Figure 6.31 Solid-liquid phase change in a rectangular cavity filled with porous media (Beckermann and Viskanta, 1988; Reproduced with permission from Elsevier). 490 Transport Phenomena in Multiphase Systems (a) Melting (Th=45 °C, Tc=20 °C) (b) Solidification (Th=40 °C, Tc=20 °C) Figure 6.32 Comparison between numerical and experimental results (Beckermann and Viskanta, 1988; Reproduced with permission from Elsevier). vertical. The reason for the weak effect of natural convection is that the Rayleigh numbers – as defined in eq. (6.367) – for the cases shown in Fig. 6.32(a) and (b) are respectively 9.22 and 11.52 (Nield and Bejan, 1999), which are comparable to Fig. 6.28(a). Oosthuizen (1988) investigated natural convection and heat transfer in an enclosure containing a porous medium saturated with the PCM [see Fig. 6.33(a)]. The left wall of the cavity is maintained at a temperature Th, which is above the melting point of the PCM, while the right wall temperature is kept at Tc, a temperature below the melting point. The problem under the configuration shown in Fig. 6.33(a) can reach steady-state, at which point the solid-liquid interface ceases to move. At steady-state, the location of the interface is not affected by the latent heat of the PCM, hsA , because there is no phase change in the steady state. The location and shape of the interface are affected by natural convection in the liquid region as well as conduction in the solid region. Oosthuizen (1988) solved the phase change problem shown in Fig. 6.33(a) by using a finite element method. The effects of Rayleigh number (Ra = g KH(Th –Tc)/ (ν Aα m ) =0 ~ 500), aspect ratio (H/L = 0.5 ~ 2), and thermal conductivity ratio of the liquid and solid region ( keff ,A / keff , s = 1  3 ) on the fluid flow and heat transfer were investigated. The effect of the Rayleigh number on the solid-liquid interface for keff ,A / keff , s = 1  3 and θ c = (Tm − Tc ) /(Th − Tc ) = 0.5 is shown in Fig. 6.33(b). In the absence of natural convection (Ra = 0), the phase change is conductioncontrolled and the location of the interface is located at L/2. After the Rayleigh number is greater than 30, the interface becomes inclined and the total volume of the liquid region increases with increasing Rayleigh number. Chapter 6 Melting and Solidification 491 Adiabatic Th Tc Solid y Liquid 0 x Adiabatic L (b) Effect of Ra on the freezing front (a) Physical model Figure 6.33 Steady-state convection and heat transfer in a porous medium with differentially heated side walls (Oosthuizen, 1988). 6.9 Applications of Solid-Liquid Phase Change 6.9.1 Latent Heat Thermal Energy Storage Phase-change thermal energy storage systems can store thermal energy while being subjected to heat input, and then release it to the environment over a long period of time. Therefore, they are especially suitable for space applications involving pulsed power loads, such as a large amount of heat rejection from a power cycle in a short period of time. Heat transfer and its enhancement during melting and solidification in shell-and-tube latent heat thermal energy storage systems, similar to that shown in Fig. 1.17, have been investigated numerically (Cao et al., 1991; Cao and Faghri, 1991; 1992) and analytically (Zhang and Faghri, 1995a, b; 1996a, b; 1997). In this subsection, heat transfer during solidification around a horizontal tube with internal convective cooling (Zhang et al., 1997) will be discussed. The theoretical model employed in this study is shown in Fig. 6.34. At the very beginning of the process (t = 0), the tube, which has a radius of R0, is surrounded by a liquid phase change material with uniform temperature Tf > Tm. The temperature of the working fluid inside the tube is Ti, and the convective heat transfer coefficient between the working fluid and the internal tube wall is hi. Both hi and Ti are kept constant throughout the process. The thickness of the tube is assumed to be very thin, so the thermal resistance of the tube wall can be neglected. The phase change material can be treated as if it were directly in contact with the coolant inside the tube. Liquid adjacent to the cooled surface will be solidified, and the temperature difference between the solid-liquid 492 Transport Phenomena in Multiphase Systems R0 R θ hi, Ti δ y x g Figure 6.34 Solidification around a horizontal tube. interface and the otherwise quiescent liquid will drive natural convection in the liquid region. The liquid region is assumed to be sufficiently large so that the convection can be treated as natural convection between the surface and an extensive fluid medium. Conduction is the only heat transfer mode in the solidified region. The analysis is based on several additional assumptions: 1. The liquid is Newtonian and Boussinesq as well as incompressible. The Prandtl number of the liquid phase change material is greater than unity. 2. The solid is homogeneous and isotropic. 3. The liquid motion induced by volumetric variation during solidification is neglected, i.e., the density of the solid is equal to the density of the liquid. In addition, the phase change material’s properties are constant in the liquid and solid regions. 4. The solid-liquid interface is assumed to be a smooth cylinder concentric with the cooled tube. 5. The effect of natural convection is restricted within the boundary layer and the bulk liquid has a uniform temperature, Tf. Compared to the thickness of the solidified layer, the thickness of the natural convective boundary layer on the solid-liquid interface is very thin, except at the onset of freezing (the boundary layer thickness in Fig 6.34 is significantly exaggerated so as to be clearly visible). However, since the onset of freezing is a very short period compared to the whole solidification process, it is reasonable to neglect the curvature effect in the equations of the boundary layer. The boundary layer equations of the problem can be written as follows: ∂u ∂v + =0 (6.386) ∂x ∂y Chapter 6 Melting and Solidification 493 ν ∂ 2u + g β (T − Tm )sin θ = 0 ∂y 2 (6.387) ∂T ∂ (uT ) ∂ (vT ) ∂ 2T + + =αf 2 (6.388) ∂t ∂x ∂y ∂y Since the Prandtl number of the liquid is much larger than unity, the inertia terms in the momentum equation have been neglected. These equations were solved by an integral method. Assuming a polynomial profile, the temperature and velocity profiles inside the boundary layer of thickness δ are expressed as y· § T = T f − (T f − Tm ) ¨ 1 − ¸ © δ¹ 2 2 (6.389) y§ y· 1− (6.390) δ¨ δ¸ © ¹ where U is a characteristic velocity, δ is the boundary layer thickness, and both of them are functions of time t and angle θ. Integrating both sides of eq. (6.387) with respect to y within the interval (0, δ ), the expression for U in eq. (6.390) should be g β sin θ U= (6.391) (T f − Tm )δ 2 3ν Integrating eq. (6.388) in the same manner, one obtains 1 ∂δ dR 1 g β sin θ ∂ 3δ 2α f − + − (T f − Tm ) + =0 (6.392) 3 ∂t dt 90 ν R ∂θ 3 δ Heat transfer in the solidified layer is dominated by conduction. The governing equation of the solid layer and corresponding boundary conditions are as follows: 1 ∂ § ∂T · 1 ∂T R0 < r < R t > 0 (6.393) ¨r ¸= r ∂r © ∂r ¹ α s ∂t ∂T ks = hi (T − Ti ) r = R0 t > 0 (6.394) ∂r Ts (r , t ) = Tm r = R t > 0 (6.395) u =U ks ∂Ts ∂r − kf r=R ∂T ∂y = ρ hsA y =0 dR dt r =R t >0 (6.396) Equations (6.392) – (6.396) can be nondimensionalized and solved using integral approximate method. When Bi → ∞ , it corresponds to boundary conditions of the first kind, i.e., the tube wall temperature is Ti and is kept steady throughout the process. Wang et al. (1991) experimentally investigated the solidification process around a horizontal cooled tube. A comparison of the predicted solidification rate, V / V0 = ( R / R0 ) 2 and the experimental results is shown in Fig. 6.35. When Ra = 0, i.e., no superheat exists in the liquid region or the solidification process is dominated by conduction, the predicted value is 18% 494 Transport Phenomena in Multiphase Systems Figure 6.35 Comparison of predicted solidification rate with experiments (Zhang et al., 1997). lower than the experimental data. During the conduction-dominated freezing process, the front of the freezing layer is a dendritic layer. Therefore, it is believed that the solid-liquid interface is extended by the dentric layer. As the Rayleigh number increases, natural convection occurs in the liquid region, and the solid-liquid interface becomes smooth because of the natural motion of the liquid. When Ra = 1.8 × 105 , the predicted value is only 8% lower than the experimental data, so the agreement is satisfactory. 6.9.2 Heat Pipe Startup from Frozen State When a high-temperature heat pipe starts from room temperature, the working fluid within the wick structure is in the frozen state. During startup from the frozen state, large thermal gradients generate significant internal stresses within the pipe wall. These stresses may severely shorten the life of the heat pipe. As a result, information on the stresses is needed for design purposes, and it is necessary to determine the temperature distribution during the transient frozen startup period. During frozen startup of cylindrical heat pipes, the wall temperature drops sharply from the active hot region to the inactive cold region. The hot zone length increases over time until it reaches the condenser end cap, at which time the frozen startup period is completed. Therefore, the startup process is characterized by the movement of a high-temperature front down the length of the pipe. If the input heat flux at the evaporator is sufficient, the heat pipe working temperature will continue to rise after the front reaches the condenser end cap until a steady-state condition is attained. In the hot zone, the vapor is in the continuum state, while the vapor in the cold region is in the rarefied Chapter 6 Melting and Solidification 495 Figure 6.36 Evolution of heat pipe startup process from the frozen state (not to scale; Cao and Faghri, 1993a). condition. Therefore, both types of vapor flows must be considered simultaneously and linked at the interface between the hot and cold regions. The startup process of a liquid metal heat pipe from the frozen state may be divided into several periods for convenience of analysis (Cao and Faghri, 1993a): 1. The heat pipe wick is initially in the frozen state. The heat pipe core can be considered to be completely evacuated before startup. The heat pipe in its initial state is illustrated in Fig. 6.36 (period 1). 2. Power is applied to the evaporator section of the heat pipe. Heat conduction through the wall and melting of the working fluid in the porous wick take place. However, the liquid-solid melting interface has not yet reached the wick/vapor boundary, and the heat pipe core is still evacuated because no evaporation has taken place in the evaporator section (Fig. 6.36, period 2). 496 Transport Phenomena in Multiphase Systems 3. In the evaporator section, the working fluid in the wick is completely melted, and evaporation takes place at the liquid-vapor interface. The vapor pressure is so low that the vapor flow is in a rarefied or free molecular condition. Also, some working fluid in the wick is still in the solid state in the adiabatic and condenser sections of the heat pipe (Fig. 6.36, period 3). 4. As evaporation continues, the amount of vapor accumulated in the core is large enough in the evaporator that a continuum flow is established there. Near the condenser end of the heat pipe, however, the vapor flow is still in the rarefied or free molecular regime. This period may be called the intermediate period (Fig. 6.36, period 4). 5. The working fluid in the wick is completely melted and continuum flow is established over the entire heat pipe length. The heat pipe continues to operate and gradually reaches steady state (Fig. 6.36, period 5). Since most heat pipes have a long, slender configuration, and during startup the heat pipe wall temperature drops sharply in the axial direction from the hot zone to the cold zone, the heat pipe startup process is actually the movement of the hot temperature front from the evaporator region to the condenser cap. Therefore, most of the startup time is during the intermediate period (period 4). In addition, when the heat input at the evaporator is high enough – after the hot zone reaches the condenser end cap and the continuum vapor flow is established over the entire heat pipe length – the heat pipe working temperature will continue to rise until a steady-state condition is reached. During startup, continuum and rarefied vapor flows coexist in different portions of the heat pipe, so these two different flows should be considered simultaneously, and be coupled to simulate the entire startup process. Cao and Faghri (1993b) proposed a two-region model to link the continuum and rarefied vapor flows in the heat pipe. In the continuum flow region, the continuum momentum and energy equations are solved; in the rarefied flow region, the self-diffusion model proposed by Cao and Faghri (1993a) is applied. For the porous wick structure, the change of phase of the frozen working fluid must be considered. The wick structure is assumed to be isotropic and homogeneous. The temperature transforming model (Section 6.5.4) is employed to study melting of working fluid during the heat pipe startup. The energy equation in the wick structure (a porous medium) is ∂ ∂§ ∂T · 1 ∂ § ∂T · ∂ ( ρ eff cT ) = ¨ keff ¸+ ¨ keff r ¸ − ( ρ eff b) (6.397) ∂t ∂z © ∂z ¹ r ∂r © ∂r ¹ ∂t The coefficients c, b, and keff in eq. (6.397) are T < Tm − ΔT ­cse ° c(T ) = ®cme + hsA /(2ΔT ) Tm − ΔT < T < Tm + ΔT (6.398) °c T > Tm + ΔT ¯ Ae Chapter 6 Melting and Solidification 497 ­cse (ΔT − Tm ) ° b(T ) = ®cme ΔT + hsA / 2 − ( cme + hsA /(2ΔT ) ) Tm °c ΔT + h − c T sA Ae m ¯ se T < Tm − ΔT Tm − ΔT < T < Tm + ΔT T > Tm + ΔT (6.399) k se T < Tm − ΔT ­ ° keff = ®k se + (kAe − k se )(T − Tm + ΔT ) /(2ΔT ) Tm − ΔT < T < Tm + ΔT (6.400) °k T > Tm + ΔT ¯ Ae where ρeff = ερA + (1 − ε ) ρ ws , cse = ε cs + (1 − ε )cws and cAe = ε cA + (1 − ε )cws are effective density, and effective specific heats of wick saturated with a working fluid in frozen or liquid states. The specific heat in the mushy zone is cme = 0.5(cse + cAe ) . The values of kse and kAe depend on the particular heat pipe wick structure. Figure 6.37 compares the numerical and experimental outer wall temperatures for high temperature heat pipes with multiple heat sources and sinks during all of the frozen startup periods for Cases 11a-11f (Faghri et al., 1991) with evaporator 1 active, which is closest to the evaporator end cap. The agreement between the numerical solutions and the experimental data is generally good. The total heat inputs in the different cases were maintained at different levels in the active evaporator and incremented from 311 W (Case 11a) to 598 W (Case 11f). The corresponding power outputs at the condenser were 119 W (Case 11a) and 298 W (Case 11f). The hot zone moved from the active evaporator down the length of the heat pipe with time. At t = 185 min., which Figure 6.37 Outer wall temperature along the heat pipe compared with experimental data for case 11a-11f (Cao and Faghri, 1993b). 498 Transport Phenomena in Multiphase Systems corresponds to the steady state for Case 11a, the hot zone reached the condenser section, but the low heat input prevented the hot zone from reaching the condenser end cap. At t > 185 min., the heat input was increased to the higher level in Case 11b, and the wall temperature increased. At the heat input corresponding to Case 11c, the hot zone filled the heat pipe completely, and continuum vapor flow was established along the entire heat pipe length. Cases 11a to 11c covered the frozen startup process. After Case 11c, the performance is the transient heat pipe operation without phase change from the frozen state. For the numerical results at the steady state, the total heat output at the condenser was within 1.5% of the input at the evaporator in all cases. 6.9.3 Thermal Protection from Intense Localized Heating In some applications, a surface is hit by a moving high intensity heat source, as shown in Fig. 6.38. The major concern here is how to protect the surface from being burned out by the moving heat flux. This is indeed a concern in laser thermal threats and re-entry situations. Ablation and heat pipe technologies are usually sources of protection for surfaces in danger of being burned out by a high heat flux. Phase-change materials (PCMs) have a large melting heat, so they offer an efficient means of absorbing the heat energy while the materials are subjected to heat input, and releasing it afterward at a relatively constant temperature. Figure 6.38 Schematic of localized heating. Figure 6.39 Configuration of thermal protection from localized heating using PCM (Cao and Faghri, 1990b; Reprinted with permission from Elsevier). Chapter 6 Melting and Solidification 499 Another alternative, as shown in Fig. 6.39 (Cao and Faghri, 1990b), advantageously incorporates the merits of the above technologies as protection for surfaces attacked by a high heat flux. There, the incoming heat input moves along the surface with speed U and the heat is conducted through the outside wall to the PCM. The PCM beneath the surface melts and absorbs a large amount of the incoming heat. Because of the large melting latent heat of the PCM and the constant melting temperature Tm, the peak wall temperature will be maintained at a temperature moderately higher than Tm. With a low or moderate Tm, the reduction of the peak wall temperature is evident. The dividing sheet, the soft insulating material and the supporting plate may be used to prevent the PCM from separating from the surface wall during the melting process, and to prevent the incoming heat from being conducted into the cabin. An analysis in a fixed-coordinate system is difficult for this problem. It is convenient to study it in a moving coordinate system where the origin is fixed at the heated spot. For an imaginary observer riding along with the incoming heat beam, the outside wall and PCM will travel at the same speed, –U. The energy equation in the Cartesian coordinate system for this problem is ∂h ∂h ∂ § ∂T · ∂ § ∂T · ∂ § ∂T · = ¨k ρ − ρU (6.401) ¸ + ¨k ¸ + ¨k ¸ ∂t ∂x ∂x © ∂x ¹ ∂y © ∂y ¹ ∂z © ∂z ¹ where the second term on the left-hand side is a convection term, while the terms on the right-hand side are diffusive terms. Since the heat source is circular, Cao and Faghri (1990b) transformed eq. (6.401) into the form for cylindrical coordinate system and further rewrote the energy equation using the enthalpy transforming model (see Section 6.5.2), i.e., ∂ ( ρ h) 1 ∂ (rvr ρ h) 1 ∂ (vθ ρ h) + + ∂t r ∂r r ∂θ 1 ∂ ª ∂ ( Γ h ) º 1 ∂ ª 1 ∂ ( Γh ) º ∂ ª ∂ ( Γh ) º (6.402) = + + r r ∂x « ∂r » r ∂θ « r ∂θ » ∂z « ∂z » ¬ ¼ ¬ ¼ ¬ ¼ 1 ∂ ª ∂S º 1 ∂ ª 1 ∂S º ∂ 2 S + + + r r ∂x « ∂r » r ∂θ « r ∂θ » ∂z 2 ¬ ¼ ¬ ¼ where vr = −U cos θ and vθ = U sin θ are the velocity components in the cylindrical coordinate system, and Γ and S for the PCM can be obtained from eqs. (6.231) and (6.232). For the wall, no phase change occurs, and therefore h = cwT , Γ = k w / cw and S = 0 . Equation (6.402) with corresponding boundary conditions is solved using a finite volume method. The resulting algebraic equations are solved with the Gauss-Seidel iterative method. Figure 6.40 shows the effect of the Stefan number [ Ste = cs (Tm − Ti ) / hsA ] on the steady-state wall temperature, (T–Ti)/(Tm–Ti), at steady state along the x-z * plane. Other parameters are defined as δ * = δ / Rh , δ p = δ p / Rh , Nv = Ucw(Tm– ′′ ′′ Ti)/ qh , N Q = qh Rh /(Tm − Ti )k w , Csw=( c)s/( c)w, CAw = ( ρ c)A /( ρ c) w , Ksw=ks/kw ′′ and K Aw = kA / kw where δ and δ p are the thicknesses of the wall and PCM, qh 500 Transport Phenomena in Multiphase Systems Figure 6.40 Steady-state wall temperature (Cao and Faghri, 1990b; Reprinted with permission from Elsevier). is heat flux of the heat source (flat hat distribution), and Rh is the radius of the heat source. Also shown in the figure is a pure wall (without PCM) having the same weight per unit surface area as that of the wall-PCM module. As can be seen, the reduction of the peak wall temperature is significant. In the figure, * δ m = δ m / Rh is the dimensionless maximum melting front depth, where δ m is the maximum melting front depth from the wall surface corresponding to the steady state. 6.9.4 Microwave Thawing of Food and Biological Materials Food thawing processes are becoming more important, because the demand for frozen food products continuously increases. The interest in microwave thawing is stimulated by its ability to avoid some of the disadvantages associated with conventional thawing: long processing times, large space requirements, microbial problems, chemical deterioration, drip loss, and high fresh water consumption (Zeng and Faghri, 1994c). During microwave heating, the transient temperature distribution within the material is determined by the internal heat generation attributed to dissipation of electrical energy from microwave radiation, and the heat transfer by conduction and convection. The moisture transfer within the material and evaporation at the surface can also influence the temperature profile. Modeling microwave thawing heat transfer is difficult due to the effects of highly nonlinear phenomena, such as the rate of energy dissipation and the energy distribution within the material. These phenomena are governed by the thermal, electrical, and physical properties of the material and vary with temperature during the microwave thawing process. In addition, for most food and biological Chapter 6 Melting and Solidification 501 materials, phase change occurs over a range of temperature, which typically results in mushy region phase change effects. Moreover, the thermal and electrical properties of moderate and high moisture content food materials subjected to microwave radiation vary considerably within the temperature range of the mushy phase. Zeng and Faghri (1994c) presented a two-dimensional mathematical model to deal with the complicated thawing process, which is presented here. The model includes a complete description of the frozen, mushy, and thawed phases and the evaporation of water. Microwaves are electromagnetic waves with wavelengths ranging from 1 cm to 1 m. The heat generation within the microwave-irradiated frozen materials results from dipole excitation and ion migration. There may also be other mechanisms of interaction between the material and the electromagnetic field which cause the dissipation and heating effects of microwave energy (Datta, 1990). The equations governing the absorption of microwave radiation by a conducting material are the Maxwell equations of electromagnetic waves. In the time-varying case, the different field equations are coupled, since a changing magnetic flux induces an electric field and a time-varying electric flux induces a magnetic field. If it is assumed that the currents in the material obey Ohm's law, these equations can be written as ∂B (6.403) +∇×E = 0 ∂t ∂E ε (6.404) +σE = ∇× ∂t ∇ ⋅ (ε ) = 0 (6.405) ∇⋅B = 0 (6.406) B = μH (6.407) where E is the electric field strength, B is the magnetic flux density, H is the magnetic field intensity, and , , and are magnetic permeability, electric permissivity and conductivity of the material, respectively. Maxwell’s first equation is Faraday’s law of induction, which states that a time variation of flux density is accompanied by the curl of the electrical field. Maxwell’s second equation refers to Ampere’s law, whose integral form states that the magnetic field over a closed path is equal to the enclosed current. Equations (6.405) and (6.406) are known respectively as Gauss’ magnetic law and Gauss’ electric law, while eq. (6.407) is the constitutive relation for a simple medium. For the microwave thawing process, the energy equation, including a heat generation term, is given by ∂h = ∇ ⋅ (k ∇T ) + q′′′ ∂t (6.408) where h is enthalpy and k is the thermal conductivity. The heat generation term q′′′ accounts for the conversion of microwave energy to heat energy. Its relationship to the electrical field intensity E at that location can be derived from eqs. (6.403) – (6.407) (Metaxas and Meridith, 1983): 502 Transport Phenomena in Multiphase Systems q′′′ = 2π f ε 0ε ′′ E (6.409) where the magnetic losses in the food material have been ignored. The microwave frequency f is 2450 MHz for most domestic ovens and 915 MHz for some industrial equipment. ε 0 is the dielectric constant of free space, and ε ′′ is the loss factor for the dielectric being heated. The basic assumptions concerning the heat generation term include the following: 2 1. The microwaves are planar and propagate perpendicularly to the material. 2. The microwave field at the material surface is uniform. 3. Microwave energy entering from different sides is considered separately, and decreases exponentially from the surface into the material. 4. The total amount of heat absorbed by a sample during a heating cycle is constant. Generally, it is difficult to either measure or calculate the strength of the electrical field inside the material from Maxwell’s equations, so an exponential variation in the field intensity from the surface of the material into its interior is generally assumed (Datta, 1990): 2 Ex = E02 exp(−2α x) (6.410) Substituting eq. (6.410) into eq. (6.409), we have q′′′ = q′′′0 exp(−2α x) (6.411) x x ′′′ where qx 0 is the rate of heat generation corresponding to an electric field E0 at the surface of the material, and x is the distance from the surface into the material. The heat generation term can be expressed in general as ′′′ q′′′ = qx 0 exp( −2α x) + q′′′0 exp( −2α y ) + q′′′ exp(−2α z ) (6.412) y z0 ′′′ y where qx 0 , q′′′0 and q′′′ are the rates of heat generation corresponding to an z0 electric field E0 at the surfaces of the material in the different coordinate directions. The attenuation coefficient α is calculated from the dielectric constant ε ′ , the dielectric loss factor ε ′′ , and the wavelength λ0 in vacuum, which determines the energy distribution within the material (Mudgett, 1986): 2 λ0 where the dielectric constant ε ′ is a measure of the ability of a material to store electrical energy; it varies significantly with composition and temperature. Since phase change occurs within a range of temperature, mushy zone exists between the solid and liquid regions. The thermal and dielectric properties in the mushy zone can be interpolated according to the liquid fraction (Zeng and Faghri, 1994c). The temperature transforming model for a binary system (see Section 6.6.3) can be employed to solve the thawing process. α= 2π ε ′{[1 + (ε ′′ / ε ′) 2 ]1/ 2 − 1} (6.413) Chapter 6 Melting and Solidification 503 Figure 6.41 Comparison of numerical and experimental temperatures during microwave thawing of a short Tylose cylinder (D =4.8 cm): (a) H/D =1.0, P =351.0 W, (b) H/D =1.5, P =367.0 W, and (c) H/D =2.0, P =400.0 W (Zeng and Faghri, 1994c). Figure 6.41 presents the thawing history for the 48-mm-diameter sample with H / D = 1, 1.5, and 2. In Fig. 6.41(a), H / D = 1 , for the time interval of 0 < t ≤ 10 s , predicted temperatures vary rapidly with time. This phenomenon is attributed to the constant absorbed power assumption in the mathematical model. For the time interval of 10 < t ≤ 80 s , the temperature data at the three locations are nearly the same, which suggests the temperatures at r = 0, R / 2, and 3R / 4 range from –8 to –0.6°C, which is defined as the mushy region temperature range. As the thawing process is completed in the outer layer of the sample, the predicted and measured temperatures differ more significantly. Figure 6.41(b) compares the experimental and numerical results for D×H = 48×72 mm at selected locations. From t = 20 to 100 s, the temperatures indicate the existence of the mushy region at r = 0, R/2, and 3R/4. Achievement of better accuracy is limited by the local completion of thawing. In the case of 2450 MHz, the penetration depth of microwaves in the frozen sample is approximately 5.9 cm. However, in thawed Tylose, the penetration depth is approximately 1.45 cm. As 504 Transport Phenomena in Multiphase Systems time progresses and temperatures increase, the thawed layer absorbs larger amounts of energy, while less energy reaches the mushy or frozen regions. Similar trends in the thawing history can be found in Fig. 6.41(c) for D×H = 48×96 mm. A comparison of the predicted and experimental thawing times shows good agreement. At t = 60 s, the mushy region extends over the entire domain. 6.9.5 Laser Drilling Laser drilling is very important in many industries, including automotive, aerospace, electronics, and materials processing. It is a very complex process, as both melting and vaporization are accomplished. Considerable research has been carried out to develop a theoretical laser drilling model, and many detailed investigations are in the existing literature (Ganesh et al., 1997a). A very detailed numerical model and complete physical model for the laser drilling process was developed by Ganesh et al. (1997a, b) employing a free surface and phase change simulation. Zhang and Faghri (1999) modeled melting and vaporization, conjugated with conduction in the workpiece. Figure 6.42 shows the physical model of the laser drilling process. A laser beam with an intensity of I(r,t) – usually Gaussian – is produced and directed toward a solid target material at an initial temperature of Ti < Tm, which absorbs a fraction of the incident light energy. The laser beam produces a very intensive Figure 6.42 Physical model of laser drilling process. Chapter 6 Melting and Solidification 505 heat flux into the solid, which brings the temperature of the solid to its melting point, melts the solid, and then vaporizes the resulting liquid. The following assumptions will be needed to analyze the laser drilling process: 1. The target material is opaque. 2. No plasma is generated in the laser-drilled hole. 3. Liquid metal flow inside the hole is neglected and the temperature distribution in the liquid metal is assumed to be linear; to account for the unsteady state effect, the sensible heat required to bring the liquid layer temperature to the average value of saturation temperature and melting point is added to the latent heat of vaporization, i.e., hA′v = hAv + c pA (Tsat − Tm ) / 2 . 4. The liquid and solid properties are independent of temperature and the densities of the liquid and solid phases are assumed to be identical. 5. Heat loss to the environment due to convection and surface radiation is neglected. The laser-material interaction can be divided into three stages: (1) preheating, (2) melting, and (3) vaporization. Assuming vaporization has started and the liquid-vapor interface is formed, the geometric shapes of the liquid-vapor and solid-liquid interfaces are respectively expressed as: f1 ( z , r , t ) = z − s1 (r , t ) = 0 (6.414) f 2 ( z , r , t ) = z − s2 ( r , t ) = 0 (6.415) where s1 (r , t ) and s2 (r , t ) are locations of interfaces between solid-liquid and liquid vapor, respectively. The temperatures at the two interfaces satisfy the following conditions: T = Tm , z = s1 ( r , t ) (6.416) T = Tsat , z = s2 (r , t ) (6.417) The energy balance at the liquid-vapor and solid-liquid interfaces can be used to obtain the locations of the two interfaces (Zhang and Faghri, 1999): 2 r2 −2 T − Tm ª § ∂ s1 · º ∂s ′v 1 = α a I 0 e R − kA sat (6.418) ρ hA «1 + ¨ ¸ » , z = s1 (r , t ) ∂t s2 − s1 « © ∂ r ¹ » ¬ ¼ 2 T −T ºª § ∂ s · º ∂ s ª ∂T ρ hsA 2 = « ks s + kA sat m » «1 + ¨ 2 ¸ » , z = s2 (r , t ) (6.419) ∂t ¬ ∂ z s2 − s1 ¼ « © ∂ r ¹ » ¬ ¼ where the vaporization temperature, Tsat, equals the saturation temperature corresponding to the recoil pressure. Ganesh et al. (1997a) derived an equation of saturation temperature by a using gas dynamic model, i.e., ªh § 1 § ∂T · γ +1 1 ·º (6.420) γ Rg Tsat ¨ I abs − kA A ¸ = p0 exp « Av ¨ − ¸» ¨ ¸ ∂n1 ¹ γ hAv « » © ¬ Rg © Tsat ,0 Tsat ¹ ¼ 506 Transport Phenomena in Multiphase Systems where Iabs is the rate of energy absorption along the normal direction of the liquid-vapor interface, n1, and γ = c p ,v / cc ,v is the ratio of specific heats for the metal gas. In order to find a simple approximate solution of conduction heat loss, it is assumed that conduction takes place only in the local surface normal, i.e., the conduction in the solid is assumed to be locally one-dimensional. The conduction equation in the solid phase was given in a transformed coordinate system, n2, by Zhang and Faghri (1999) and solved using an integral approximate solution. The laser drilling on a workpiece made of Hastelloy-X was simulated and the results were compared with the experimental data (Bellantone and Ganesh, 1991) in Fig. 6.41. In order to obtain dimensional results for comparison, the value of absorptivity of the target materials was needed; it was taken as 0.85, as recommended by Kar et al. (1992). Other dimensionless variables in Fig. 6.43 are defined as Ste = c pA (Tsat − Tm ) / hsA , Rh = hAv / hsA , H Av = hAv /( Rg Tsat ,0 ) , θ m = Tm / Tsat ,0 , and N p 0 = γ p0 c pA R /[(γ + 1)kA γ Rg Tsat ,0 . Tsat ,0 is the saturation temperature corresponding to the ambient pressure, p0 , and the τ = α A t / R 2 is the dimensionless time. The average material removal rate in the laser drilling process, which was obtained by the following integration 2πρ ∞ MR = s1 (r , t p )rdr (6.421) t p ³0 was compared with experimental data. Figure 6.43 shows the comparison of Figure 6.43 Comparison of predicted and experimental material removal rate (Zhang and Faghri, 1999). Chapter 6 Melting and Solidification 507 calculated material removal rate, for the cases with and without conduction heat loss, and that of experimental data. The experimental material removal rate was obtained by scaling micrographs of single-shot drilled holes for a pulse-on time of tp = 700 s and radius of 0.254 mm at the Pratt & Whitney drilling facility at North Haven, CT (Bellantone and Ganesh, 1991). It can be seen that the predicted materials removal rate with conduction heat loss is slightly lower than that without conduction heat loss, but the difference is less than 2%. This suggests that the overall effect of conduction heat loss on the material removal rates is not significant. As can be seen from Fig. 6.43, the material removal rate predicted by the present model is higher than that of experimental data for most cases. The possible causes of the overprediction include uncertainty of absorptivity and possible partial laser beam blockage due to plasma, which is not taken into account in the model. Considering these complicated phenomena, the agreement between the calculated results and the experimental data is very good. 6.9.6 Selective Laser Sintering (SLS) of Metal Powders SLS of metal powder involves fabrication of near full density objects from powder via melting induced by a directed laser beam (generally CO2 or YAG) and resolidification. For sintering of metal powder, the latent heat of fusion is usually very large, and therefore melting and resolidification phenomena have a significant effect on the temperature distribution in the parts and powder, the residual stress in the part, local sintering rates, and the final quality of the parts. A significant change of density accompanies the melting process, because after melting, the volume fraction of gas(es) in the powder decreases from a value as large as 0.6 to nearly zero. In addition, the liquid metal infiltrates the unsintered region due to capillary and gravity forces. In order to overcome balling phenomena caused by surface tension, a powder mixture containing two powders with significantly different melting points can be used in the SLS process. Only the low melting point powder will be molten and resolidified during the SLS process. The low melting point powder melts and infiltrates the unsintered region due to capillary and gravity forces. The solid particles of the high melting point powder may also move downward, because the high melting point powder cannot sustain the powder bed alone. It is very clear that both liquid and solid in the powder bed have their own velocities, and these velocities may have a significant effect on the energy transport in the powder bed. Zhang et al. (2000) developed a very detailed three-dimensional thermal model of SLS of a metal powder bed that contains a mixture of two powders with significantly different melting points. The predicted results were compared with experimental results obtained with nickel braze and AISI 1018 steel powder. The physical model of the problem is shown in Fig. 6.44. A powder bed, which contains two powders with significantly different melting points, with a uniform initial temperature, Ti, below the melting point of the low melting point powder, Tm, is in a cavity with a size of xD×2yD×zD (length × width × height). 508 Transport Phenomena in Multiphase Systems Figure 6.44 Physical model of three-dimensional sintering of two-component metal powder (Zhang et al., 2000). The coordinate system is also shown in Fig. 6.44. A round Gaussian laser beam scans the surface of the cavity starting from a point x = x0 toward the positive direction of the x-axis with a constant velocity, ub. As the laser beam interacts with the powder, the temperature of the powder increases to Tm, which induces melting. The molten metal infiltrates the unsintered region of the powder bed due to the driving forces of capillarity and gravity. The top surface of the melted powder bed recedes, since the low melting point material is molten and the high melting point material alone cannot sustain the powder bed. After the laser beam moves away, the liquid pool cools and resolidifies to form a densified heat affected zone (HAZ). These zones, formed by multiple laser scans, are subsequently interwoven to form a layer part. Since the melting and nonmelting powder particles have sizes of the same order of magnitude, the nonmelting powder’s skeletal structure collapses after the low melting point powder is liquefied. In this case, the conservation of mass principles can be used to show that the porosity of the powder bed remains constant (while the volume of the powder bed shrinks; Pak and Plumb, 1997). The system is symmetric about the (x,z) plane; thus, only half of the cavity ( 0 < y < yD ) needs to be studied. The intrinsic phase-averaged liquid velocity, VA , must satisfy the continuity equation: ∂ε A 0 (6.422) + ∇ ⋅ (ε A VA ) = Φ A ∂t Chapter 6 Melting and Solidification 509 where VA = VA A is the intrinsic phase-averaged velocity of the liquid phase, ε A is volume fraction of the liquid low melting metal, and Φ 0 is the volumetric A production rate of the liquid due to melting. Likewise, the solid low melting point powder vanishes at the same rate, i.e., ∂ε s ∂ (ε s ws ) + = −Φ 0 (6.423) A ∂t ∂z where ε s is volume fraction of the solid low melting metal, and ws is the shrinkage velocity in the z-direction. The continuity equation for the high melting point material is ∂ε H ∂ (ε H ws ) + =0 (6.424) ∂t ∂z where ε H is volume fraction of the high melting point powder. The volume fraction of each species satisfies ε p + ε s + ε H = 1 , where the porosity of the powder bed, ε p , is defined as the total volume of void – including the volumes of gas and liquid – relative to the total volume of the powder bed ( ε p = ε A + ε g ). The solid velocity can be determined by integrating eq. (6.424), and the result is z>s ­0 ° (6.425) ws = ® ε s ,i ∂s z<s ° ¯ 1 − ε ∂t where ε s ,i is initial porosity of the loose powder. The liquid flow occurs in three directions, and the velocities can be determined using Darcy’s law: KK rA VA − ws k = (∇pc + ρ A gk ) (6.426) εAμ where k is unit vector in the z-direction, K is permeability of the porous medium, K rA is relative permeability, and pc is capillary pressure. The temperature transforming model using a fixed grid method is employed to describe melting and resolidification in the powder bed. Assuming that the thermal equilibrium exists between two different powders, the energy equation is ∂ ªε H ρ H c pH + (ε A + ε s ) ρ L c pL º T + ∇ ⋅ (ε A VA ρ L c pLT ) ¼ ∂t ¬ ∂ (6.427) + ª ws (ε H ρ H c pH + ε s ρ L c pL )T º = ∇ ⋅ (k ∇T ) ¼ ∂z ¬ ∂ ­∂ ½ − ρ L ® [ (ε A + ε s )b ] + ∇ ⋅ (ε A VA b) + ( ε s ws b ) ¾ ∂z ¯ ∂t ¿ where the effective specific heat for the low melting point metal, cpL and the source term b can be obtained by using eqs. (6.241) and (6.242). The thermal conductivity for the unsintered powder can be calculated using the empirical correlation proposed by Hadley (1986), while the thermal conductivity of a liquid { } 510 Transport Phenomena in Multiphase Systems or resolidified part is obtained based on the weighted average of the thermal conductivities of both powders. The difference between the liquid and sintered regions is that the temperature of the latter is lower than Tm, and it exists in the form of a solid. Obviously, the low melting point material in the sintered region cannot flow because it exists in the solid state. For the convenience of programming, the volume fraction of the low melting point metal is still represented by ε A , and the summation of ε A and ε g holds constant, ε p , even in the sintered region. The viscosity of the low melting point metal in the sintered region can be set as a very large value so that a low melting point metal velocity of zero can be achieved. The governing equations and the corresponding boundary and initial conditions can be nondimensionalized and solved numerically. The predicted results were compared with experimental results obtained using nickel braze as the low melting point powder and AISI 1018 steel as the high melting point powder. Figure 6.45a) shows the predicted shape of the sintered region in three different cross sections. The dimensionless laser intensity is defined as N i = α a I 0 R /[k H (Tm − Ti )] , where α a is absorptivity, I0 is the laser intensity at the center of the Gaussian laser beam, R is the radius of the laser beam, kH is the thermal conductivity of the high melting point powder, and Ti is the initial Figure 6.45 Comparison of cross-section area obtained by numerical simulation and experiment (Ni = 0.0749, Ub= 0.124; Zhang et al., 2000). Chapter 6 Melting and Solidification 511 temperature of the powder bed. The dimensionless scanning velocity is U b = ub R / α H where α H is the thermal diffusivity of the high melting point powder. The laser beam starts to scan the powder bed from x = 28.1 mm and beam irradiation is curtailed at x = 48.1 mm. Figure 6.45(b) shows three cross sections corresponding to the numerical results in Fig. 6.45(a). It can be seen that the shapes of the numerical results agree fairly well with the experimental results. As can be seen from Fig. 6.45(b), some voids can be observed in the cross section of the sintered part. Again, this phenomenon was not observed in the numerical results, because the porosity of the powder bed was assumed to remain constant in the process. However, the predicted and the measured total volumes of powder sintered in the process are very close. 6.10 Microscale Phase Change 6.10.1 Overview Microscale heat transfer has drawn the attention of many researchers because it is important in many advanced manufacturing and materials processes. The traditional phenomenological laws, such as Fourier’s law of heat conduction, are challenged in the microscale regime because (1) the characteristic lengths of the various heat carriers are comparable to each other and to the characteristic length of the system considered, and/or (2) the characteristic times of the various heat carriers are comparable to the characteristic energy excitation time (Wang and Prasad, 2000). Thus, microscale heat transfer can be referred to as heat transfer occurring on both the micro-length and micro-time scales. Microscale heat transfer finds applications in thin film (micro- length scale) as well as ultra-short pulsed laser processing (micro- time scale). Compared to long pulsed laser processing, short-pulsed laser processing enables users to precisely control the size of the heat-affected zone, the heat rate, and the solidification speed. The carriers of heat in solids include electrons, phonons, and photons. Shortpulsed laser melting of thin film involves the following three steps (Qiu and Tien, 1993; Kuo and Qiu, 1996): (1) absorption of photon energy by free electrons, (2) energy transfer between the free electrons and the lattice, and (3) phase change of the lattice due to the propagation of energy. The electron-lattice thermalization time is on the order of picoseconds (10-12 sec), which is significantly shorter than the time scale for most heat transfer problems. It is common practice to assume that the electrons and the lattice are always in thermal equilibrium and have the same temperature. However, when picosecond or femtosecond (10-15 sec) lasers are used in the material processing, the time scale of the laser pulse is comparable to or shorter than the electron-lattice thermalization time, and so the electrons and lattice are not in thermal equilibrium during short-pulse laser materials processing. Furthermore, rapid phase-change phenomena do not experience 512 Transport Phenomena in Multiphase Systems controlled heat transfer at the solid-liquid interface as suggested by eq. (6.8). The solid-liquid interface can be heated well above the melting point during a rapid melting process, in which case the solid becomes superheated. Similarly, the solid-liquid interface can be cooled far below the melting point in the rapid solidification process, in which case the liquid becomes undercooled. Both superheated solid and undercooled liquid are thermodynamically unstable and cannot be described using equilibrium thermodynamics. Once phase change is triggered in a superheated solid or undercooled liquid, the solid-liquid interface can move at an extremely high velocity (on the order of 10 to 102 m/s). In this case, the solid-liquid interface velocity is dominated by nucleation dynamics, which deviates significantly from equilibrium thermodynamics. 6.10.2 Two-Step Model Short-pulsed laser melting of thin film involves two nonequilibrium processes: (1) deposition of laser energy, and (2) initiation of melting. Therefore, successful modeling of short-pulsed laser melting requires modeling of these two fundamental processes. When the photons from the laser beam interact with the solid, the photon is first absorbed by the free electrons and then transferred to the lattice by interaction between the free electrons and lattice. Since the time scales of the laser pulse and the electron-lattice thermalization are comparable, the free electrons and the lattice have their own temperatures. The picoseconds laser melting of free-standing gold film is solved by Kuo and Qiu (1996). The problem can be approximated to be one-dimensional because the radius of the laser beam is significantly larger than the film thickness of the gold. The physical model of the picosecond melting is shown in Fig. 6.46, from which the energy equations of the free electrons and the lattice can be written separately as ∂T e ∂ ª e e A ∂T e º e A ′′′ C e (T e ) = « k (T , T ) » − G (T − T ) + qlaser (6.428) ∂t ∂x ¬ ∂x ¼ Figure 6.46 Laser melting of thin film. Chapter 6 Melting and Solidification 513 ∂T A ∂ § A ∂T A · A e = ¨k (6.429) ¸ + G (T − T ) ∂t ∂x © ∂x ¹ where the superscripts e and A represent the free electrons and the lattice, respectively. G is the electron-lattice coupling factor. The heat capacity of the electrons, Ce, is proportional to the temperature of only the electrons, but the thermal conductivity of the electrons depends on the temperatures of both electrons and lattice. Qiu and Tien (1993) suggest that the thermal conductivity of the electron be written as Te k e (T e , T A ) = A keq (T ) (6.430) T where keq(T) is the thermal conductivity of the electron when the electrons and lattice are in thermal equilibrium. The source terms in eq. (6.428) can be described by the following equation: 2 ªx §t·º 1− R ′′′ (6.431) qlaser = 0.94 J exp « − − 2.77 ¨ ¸ » ¨ tp ¸ » t pδ «δ © ¹¼ ¬ where R is reflectivity of the thin film, tp is laser pulse duration (s), is laser penetration depth (m), and J is laser pulse fluence (J/m2). For the conventional melting process, the velocity of the solid-liquid interface is obtained by energy balance as specified by eq. (6.15). However, this is not the case for rapid melting/solidification processes, because the velocity of the interface is dominated by nucleation dynamics. For short-pulsed laser melting of gold, Kuo and Qiu (1996) recommend that the velocity of the solid-liquid interface be obtained by § h ·ª § h T − T A ·º ds (6.432) = Vs exp ¨ − sA , a ¸ «1 − exp ¨ − sA ,a m A I ¸ » dt © kbTm ¹ « © kbTm TI ¹» ¬ ¼ where Vs is the maximum interface velocity that can be approximated as the speed of sound in the liquid phase. hsA , a is the latent heat of fusion per atom, kb is CA the Boltzmann constant, and TIA is the lattice temperature at the interface. The time t = 0 is defined as the time when the peak of a laser pulse reaches the film surface. Therefore, the initial conditions of the problem are (6.433) T e ( x, −2t p ) = T A ( x, −2t p ) = Ti The boundary conditions of the problem can be specified by assuming that the heat loss from the film surface can be neglected, i.e., ∂T e ∂T e ∂T A ∂T A = = = =0 (6.434) ∂x x = 0 ∂x x = L ∂x x = 0 ∂x x = L Equations (6.428) and (6.429) apply to short-pulsed laser melting problems, cases in which the electrons and lattice are not in thermal equilibrium. On the other hand, when the electrons and lattice are in thermal equilibrium, 514 Transport Phenomena in Multiphase Systems T e = T A = T , one can add eqs. (6.428) and (6.429) to obtain the following onestep model: ∂T ∂ § ∂T · = ¨k (6.435) C aser ¸ + ql′′′ ∂t ∂x © ∂x ¹ Figure 6.47 Comparison between two-step and one-step models ( t p = 20 ps , L = 100 nm , and R = 0.6 ; Kuo and Qiu, 1996). Chapter 6 Melting and Solidification 515 where C = C e + C A and k = k e + k A are the heat capacity and thermal conductivity of the thin film, respectively. Equation (6.435) is clearly very similar to the energy equation in Section 6.3 except for the heat source term S. In arriving at eq. (6.435), it is assumed that the electron-lattice thermalization time is significantly smaller than the laser pulse; this is the case for long-pulse laser melting. The corresponding initial and boundary conditions for eq. (6.435) are T ( x, −2t p ) = Ti (6.436) ∂T ∂T = =0 (6.437) ∂x x = 0 ∂x x = L The one-step model represented by eqs. (6.435) – (6.437) is presented here in comparison with the two-step model. Numerical solutions for both two-step and one-step models were obtained by Kuo and Qiu (1996). Figure 6.47 shows the transient melting depth, interface temperature, and interface velocity during laser melting of a gold film with a thickness of L = 1000 nm and laser pulse duration of tp = 20 ps. The results obtained by both models were presented for comparison. The two-step model predicts that the melting process is initiated near or after the end of laser irradiation, which is consistent with experimental observation. The one-step model, on the other hand, predicts that the film surface always starts melting during laser pulse irradiation, which is not true for picosecond laser melting. The superheating and undercooling predicted by the two-step model are significantly smaller than those predicted by the one-step model. Consequently, the maximum melting and resolidification velocities based on the two-step model are much smaller than those predicted by the one-step model. These results show that microscopic energy transfer and phase change are very important for picosecond laser processing of materials. References Alexiades, A., and Solomon, A.D., 1993, Mathematical Modeling of Melting and Freezing Processes, Hemisphere, Washington, DC. Bardsley, W., Hurle, D.T.J., and Mullin, T.B., 1979, Crystal Growth: A Tutorial Approach, North-Holland, Amsterdam. Beckermann, C., and Viskanta, R., 1988, “Natural Convection Solid/Liquid Phase Change in Porous Media,” International Journal of Heat and Mass Transfer, Vol. 31, pp. 35-46. Beckermann, C., and Viskanta, R., 1993, “Mathematical Modeling of Transport Phenomena during Alloy Solidification,” Applied Mechanics Reviews, Vol. 46, pp. 1-27. 516 Transport Phenomena in Multiphase Systems Beckermann, C., and Wang, C.Y., 1995, “Multi-Phase/-Scale Modeling of Transport Phenomena in Alloy Solidification," Annual Review of Heat Transfer, Vol. 6, pp. 115-198. Bellantone, R., and Ganesh, R.K., 1991, Analysis Model for Laser Hole Drilling: Final Report, Contract II, report submitted to Pratt and Whitney Aircraft, East Hartford, CT, 1991 Bejan, A., 1994, “Contact Melting Heat Transfer and Lubrication,” Advances in Heat Transfer, Vol. 24, pp. 1-38. Bejan, A., 2004, Convection Heat Transfer, 2nd ed., John Wiley & Sons, Inc, New York. Bennon, W.D., and Incropera, F.P., 1987a, “A Continuum Model for Momentum, Heat and Species Transport in Binary Solid-Liquid Phase Change Systems – I: Model Formulation,” International Journal of Heat and Mass Transfer, Vol., 30, pp. 2161-2170. Bennon, W.D., and Incropera, F.P., 1987b, “A Continuum Model for Momentum, Heat and Species Transport in Binary Solid-Liquid Phase Change Systems – II: Application to Solidification in a Rectangular Cavity,” International Journal of Heat and Mass Transfer, Vol. 30, pp. 2171-2178. Bonacina, C., Comini, G., Fasano, A., and Primicerio, A., 1973, “Numerical Solution of Phase Change Problems,” International Journal of Heat and Mass Transfer, Vol. 16, pp. 1285-1832. Braga, S. L., and Viskanta, R., 1990, “Solidification of a Binary Solution on a Cold Isothermal Surface,” International Journal of Heat and Mass Transfer, Vol. 33, pp. 745-754. Budhia, H., and Kreith, F., 1973, “Heat Transfer with Melting or Freezing in a Wedge,” International Journal of Heat and Mass Transfer, Vol. 16, pp. 195-211. Cao, W.Z., and Poulikakos, D., 1991, “Freezing of a Binary Alloy Saturating a Packed Bed of Spheres,” AIAA Journal of Thermophysics and Heat Transfer, Vol. 5, pp. 46-53. Cao, Y., and Faghri, A., 1989, “A Numerical Analysis of Stefan Problem of Generalized Multi-Dimensional Phase-Change Structures Using the Enthalpy Transforming Model,” International Journal of Heat and Mass Transfer, Vol. 32, pp. 1289-1298. Cao, Y., and Faghri, A., 1990a, “A Numerical Analysis of Phase Change Problem including Natural Convection,” ASME Journal of Heat Transfer, Vol. 112, pp. 812-815. Cao, Y., and Faghri, A., 1990b, “Thermal Protection from Intense Localized Moving Heat Fluxes Using Phase-Change Materials,” International Journal of Heat and Mass Transfer, Vol. 33, No. 1, pp. 127-138. Chapter 6 Melting and Solidification 517 Cao, Y., and Faghri, A., 1991, “Performance Characteristics of a Thermal Energy Storage Module: A Transient PCM/Forced Convection Conjugate Analysis,” International Journal of Heat and Mass Transfer, Vol. 34, pp. 93-101. Cao, Y., and Faghri, A., 1992, “A Study of Thermal Energy Storage System with Conjugate Turbulent Forced Convection,” ASME Journal of Heat Transfer, Vol. 114, pp. 1019-1027. Cao, Y., and Faghri, A., 1993a, “Simulation of the Early Startup Period of High Temperature Heat Pipes from the Frozen State by a Rarefied Vapor SelfDiffusion Model,” ASME Journal of Heat Transfer, Vol. 115, pp. 239-246. Cao, Y., and Faghri, A., 1993b, “A Numerical Analysis of High Temperature Heat Pipe Startup from the Frozen State,” ASME Journal of Heat Transfer, Vol. 115, pp. 247-254. Cao, Y., Faghri, A., and Juhasz, A., 1991, “A PCM/Forced Convective Conjugate Transient Analysis of Energy Storage Systems with Annular and Counter-Current Flows,” ASME Journal of Heat Transfer, Vol. 113, No. 1, pp. 37-42. Cho, S.H., and Sunderland, J.E., 1981, “Approximate Temperature Distribution for Phase Change of Semi-Infinite Body,” ASME Journal of Heat Transfer, Vol. 103, pp. 401-403. Crank, J., 1984, Free and Moving Boundary Problems, Clarendon Press, Oxford. Datta, A.K., 1990, “Heat and Mass Transfer in the Microwave Processing of Food,” Chemical Engineering Process, Vol. 86, pp. 47-53. Dong, Z.F., Chen, Z.Q., Wang, Q.J., and Ebadian, M.A., 1991, “Experimental and Analytical Study of Contact Melting in a Rectangular Cavity,” AIAA Journal of Thermophysics and Heat Transfer, Vol. 5, pp. 347-354. Eckert, E.R.G., and Drake, R.M., 1987, Analysis of Heat and Mass Transfer, Hemisphere, Washington, DC. El-Genk, M.S., and Cronenberg, A.W., 1979, “Solidification in a Semi-Infinite Region with Boundary Condition of the Second Kind: An Exact Solution,” Letters in Heat and Mass Transfer, Vol. 6, pp. 321-327. Faghri, A., Buchko, M., and Cao, Y., 1991, “A Study of High Temperature Heat Pipes with Multiple Heat Sources and Sinks, Part I: Experimental Methodology and Frozen Startup Profiles,” ASME Journal of Heat Transfer, Vol. 113, pp. 1003-1009. Fang, L.J., Cheung, F.B., Linehan, J.H., and Pedersen, D.R., 1984, “Selective Freezing of a Dilute Salt Solution on a Cold Ice Surface,” ASME Journal of Heat Transfer, Vol. 106, pp. 385-393. Ganesh, R.K., Faghri, A., and Hahn, Y.A, 1997a, “Generalized Thermal Modeling for Laser Drilling Process, Part I - Mathematical Modeling & 518 Transport Phenomena in Multiphase Systems Numerical Methodology,” International Journal of Heat and Mass Transfer, Vol. 40, pp. 3351-3360. Ganesh, R.K., Faghri, A., and Hahn, Y.A, 1997b, “Generalized Thermal Modeling for Laser Drilling Process, Part II - Numerical Simulation & Results,” International Journal of Heat and Mass Transfer, Vol. 40, pp. 3361-3373. Goodman, T.R. 1958, “The Heat-Balance Integral and its Application to Problems Involving a Change of Phase,” Transactions of ASME, Vol. 80, pp. 335-342. Hadley, G.R., 1986, “Thermal Conductivity of Packed Metal Powders,” International Journal of Heat and Mass Transfer, Vol. 29, pp. 909-920. Hirata, T., Makino Y., and Kaneko, Y., 1991, “Analysis of Close-Contact Melting for Octadecane and Ice inside Isothermally Heated Horizontal Rectangular Capsule,” International Journal of Heat and Mass Transfer, Vol. 34, pp. 3097-3106. Hobbs, P.V., 1974, Ice Physics, Clarendon Press, Oxford. Kar, A., Rockstroh, T., and Mazumder, J., 1992, “Two-Dimensional Model for Laser Induced Material Damage: Effect of Assist Gas and Multiple Reflection inside the Cavity,” Journal of Applied Physics, Vol. 71, pp. 2560-2569. Kuo, L.S., and Qiu, T.Q., 1996, “Microscale Energy Transfer During Picosecond Laser Melting of Metal Films,” ASME HTD-Vol. 323, Vol. 1, pp. 149-157. Jany, P., and Bejan, 1987, Melting in the Presence of Natural Convection in a Rectangular Cavity Filled with Porous Medium, Report DU-AB-6, Department of Mechanical Engineering and Materials Sciece, Duke University. Jany, P., and Bejan, 1988, “Scaling Theory of Melting with Natural Convection in an Enclosure,” International Journal of Heat and Mass Transfer, Vol. 31, pp. 1221-1235. Metaxas, A.C., and Meredith, R.J., 1983, Industrial Microwave Heating, I.E.E. Power Engineering Series, 4, Peregrinus, London. McDonough, M. W., and Faghri, A., 1993, “Ultrasonic Measurement of Interface Positions During the Solidification of an Aqueous Sodium Carbonate Solution around a Vertical Cylinder,” Experimental Heat Transfer Journal, Vol. 6, pp. 215-230. Mudgett, RE., 1986, “Microwave Properties and Heating Characteristics of Foods,” Food Technology, Vol. 40, pp. 84-93. Nield, D.A., Bejan, A., 1999, Convection in Porous Media, 2nd ed., SpringerVerlag, New York. Chapter 6 Melting and Solidification 519 Oosthuizen P.H., 1988, “The Effects of Free Convection on Steady State Freezing in a Porous Medium Filled Cavity,” ASME HTD-Vol. 96, Vol. 1, pp. 321–327. Ozisik, M.N., 1993, Heat Conduction, 2nd ed., Wiley-Interscience, New York. Pak, J., Plumb, O. A., 1997, “Melting in A Two-component Packed Bed”, ASME Journal of Heat Transfer, Vol. 119, pp. 553-559. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York. Poots, G., 1962, “An Approximate Treatment of Heat Conduction Problem Involving a Two-Dimensional Solidification Front,” International Journal of Heat and Mass Transfer, Vol. 5, pp. 339-348. Qiu, T.Q., and Tien, C.L., 1993, “Heat Transfer Mechanism During Short-Pulsed Laser Heating of Metals,” ASME Journal of Heat Transfer, Vol. 115, pp. 835841. Rathjen, K.A., and Jiji, L.M., 1971, “Heat Conduction with Melting or Freezing in a Corner,” ASME Journal of Heat Transfer, Vol. 93, pp. 101-109. Rizzo, E. M. da S., Santos, R. G. and Beckermann, C., 2003, “Modeling Solidification and Mushy Zone Deformation of Alloys,” Journal of Brazil Society of the Mechanical Science and Engineering, Vol. 25, pp. 180-184. Shamsunder, N., and Sparrow, E.M., 1975, “Analysis of Multidimensional Conduction Phase Change via the Enthalpy Model,” ASME Journal of Heat Transfer, Vol. 97, pp. 333-340. Tien, R.H., Geiger, G.E., 1967, “A Heat-Transfer Analysis of the Solidification of A Binary Eutectic System,” ASME Journal of Heat Transfer, Vol. 91, pp. 230-234. Veinberg, A.K., 1967, “Permeability, Electrical Conductivity, Dielectric Constant and Thermal Conductivity of a Medium with Spherical and Ellipsoidal Inclusions,” Soviet Phys. Dokl., Vol. 11, pp. 593-595. Viskanta, R., 1988, “Heat Transfer During Melting and Solidification of Metals,” ASME Journal of Heat Transfer, Vol. 110, 1205-1219. Wang, G.X., and Prasad, V., 2000, “Microscale Heat and Mass Transfer and nonEquilibrium Phase Change in Rapid Solidification,” Materials Science and Engineering, A., Vol. 292, pp. 142-148. Wang, Q.J., Wang, C., and Chen, Z.Q., 1991, “Heat Transfer during Superheated Solidification around Bare Tube and Finned Tubes,” Experimental Heat Transfer, Fluid Mechanics and Thermodynamics 1991, Elsevier Science Publishing Co., Inc, pp. 1171-1176. 520 Transport Phenomena in Multiphase Systems Ward, J. C., 1964, “Turbulent Flow in Porous Media,” ASCE J. Hyd. Div., Vol. 90, HY5, pp. 1-12. Zeng, X., and Faghri, A, 1994a, “Temperature-Transforming Model for Binary Solid-Liquid Phase-Change Problems Part I: Physical and Numerical Scheme,” Numerical Heat Transfer, Part B, Vol. 25, pp. 467-480. Zeng, X., and Faghri, A., 1994b, “Temperature-Transforming Model for Binary Solid-Liquid Phase-Change Problems Part II: Numerical Simulation,” Numerical Heat Transfer, Part B, Vol. 25, pp. 481-500. Zeng, X., and Faghri, A., 1994c, “Experimental and Numerical Study of Microwave Thawing Heat Transfer for Food Materials,” ASME Journal of Heat Transfer, Vol. 116, No. 2, pp. 446-455. Zhang, Y., Chen, Z., and Wang Q., 1993, “Analytical Solution of Melting in a Subcooled Semi-Infinite Sold with Boundary Conditions of the Second Kind,” Journal of Thermal Science, Vol. 2, pp. 111-115, Science Press, Beijing. Zhang, Y., Chen, Z.Q., and Faghri, A., 1997, “Heat Transfer During Solidification around a Horizontal Tube with Internal Convective Cooling,” ASME Journal of Solar Energy Engineering, Vol. 119, pp. 44-47. Zhang, Y., and Faghri, A., 1995a, “Analysis of Thermal Energy Storage System with Conjugate Turbulent Forced Convection,” AIAA Journal of Thermophysics and Heat Transfer, Vol. 9, No. 4, pp. 722-726. Zhang, Y., and Faghri, A., 1995b, “Analysis of Forced Convection Heat Transfer in Microencapsulated Phase Change Material Suspensions,” AIAA Journal of Thermophysics and Heat Transfer, Vol. 9, No. 4, pp. 727-732. Zhang, Y., and Faghri, A., 1996a, “Semi-Analytical Solution of Thermal Energy Storage System with Conjugate Laminar Forced Convection,” International Journal of Heat and Mass Transfer, Vol. 39, pp. 717-724. Zhang, Y., and Faghri, A., 1996b, “Heat Transfer Enhancement in Latent Heat Thermal Energy Storage System by Using an External Radial Finned Tube,” Journal of Enhanced Heat Transfer, Vol. 3, pp. 119-127. Zhang, Y., and Faghri, A., 1997, “Analysis of Freezing in an Eccentric Annulus,” ASME Journal of Solar Energy Engineering, Vol. 119, No. 3, pp. 237-241. Zhang, Y., and Faghri, A., 1998, “A Thermal Model for Mushy Zone Formation in Binary Solutions,” ASME Journal of Solar Energy Engineering, Vol. 120, pp. 144-147. Zhang, Y., and Faghri, A., 1999, “Vaporization, Melting and Heat Conduction in the Laser Drilling Process,” Int. J. Heat Mass Transfer, Vol. 42, pp. 1775-1790. Zhang, Y., Faghri, A., Buckley, C.W., and Bergman, T.L., 2000, “ThreeDimensional Sintering of Two-Component Metal Powders with Stationary and Moving Laser Beams,” ASME Journal of Heat Transfer, Vol. 122, pp. 150-158. Chapter 6 Melting and Solidification 521 Problems 6.1. A horizontal cylinder is embedded in infinite solid paraffin. The initial temperature of the paraffin is equal to its melting temperature and the surface temperature of the cylinder is greater than the melting point of paraffin. The melting occurs around the cylinder. The three possible shapes of solid-liquid interfaces are given in Fig. P6.1. Discuss which one should form and for what reason. If the PCM becomes ice (the surface temperature of the cylinder is less than 4 ºC), what is your answer? Figure P6.1 6.2. A horizontal cylinder is embedded in liquid paraffin with a temperature above its melting point. From time t = 0 , the surface temperature of the cylinder is suddenly decreased to a temperature below the melting point, and then the solidification process occurs around the cylinder. Discuss the effect of natural convection on the solid-liquid interface shape. 6.3. In the thermal energy storage system that is shown in Fig. 1.17, the PCM is n-Octadecane (the thermal conductivity of solid is 0.358 W/m-K) and the transfer fluid is water. The liquid with a temperature below the melting point flows through the tube. The solidification process occurs in the annulus. To enhance the heat transfer in the system, one engineer suggests installing longitudinal fins on the PCM side. Others suggest installing longitudinal fins on the fluid side. In your opinion, who is right? 6.4. The energy balance at the solid-liquid interface for a melting problem given by eq. (6.8) is obtained by assuming that Fourier’s conduction model is valid. When the temperature gradient is very high or the duration of heating is extremely short, Fourier’s conduction model is no longer valid. Instead, the heat flux and temperature gradient are related by 522 Transport Phenomena in Multiphase Systems ∂q′′( x, t ) ∂T ( x, t ) = −k ∂t ∂x where τ is the thermal relaxation time. Apply the above equation for heat flux in both solid and liquid phases and derive the energy balance at the interface. The thermal relaxation times for both phases can be assumed to be the same. 6.5. When heat flux and temperature gradient are related by the equation of the previous problem, there may be a discontinuity of temperature at the solidliquid interface (see Fig. P6.2). Modify the energy balance that you obtained in Problem 6.4 to account for the contribution of sensible heat to the interfacial energy balance. q′′( x, t ) + τ Figure P6.2 6.6. For a two-dimensional melting/solidification problem, the energy balance equation at the interface is ks ∂Ts ( x, y, t ) ∂T ( x, y, t ) − kA A = ρ hsA vn ∂n ∂n x = s ( y, t ) The above equation is not convenient for the analytical solution or numerical solution because it contains derivatives of the temperature along n -direction. Prove that the above equation can be rewritten as ª § ∂s · 2 º ª ∂T ∂T º ∂s «1 + ¨ ¸ » « ks s − kA A » = ρ hsA ∂x ¼ ∂t « © ∂y ¹ » ¬ ∂x ¬ ¼ x = s ( y, t ) 6.7. A solid PCM with a uniform initial temperature Ti below the melting point Tm is in a half-space x > 0. At time t = 0, the temperature at the boundary Chapter 6 Melting and Solidification 523 x = 0 is suddenly increased to a temperature T0 above the melting point of the PCM. The melting process occurs from the time t = 0. The densities of the PCM for both phases are assumed to be the same and the natural convection in the liquid phase can be neglected. Show that eqs. (6.35) – (6.41), which are nondimensional governing equations for solidification in a semi-infinite region, are also valid for this problem provided that subscripts 1 and 2 represent liquid and solid, respectively. 6.8. The governing equation of the one-region phase change problem discussed in Section 6.3.3, eq. (6.42), can also be reduced to an ordinary differential equation by introducing a similarity variable: η = X /(2 τ ) , which is referred to as the similarity solution. Obtain the exact solution of the oneregion problem using the similarity solution. 6.9. A solid PCM with a uniform initial temperature equal to the melting point Tm is in a half-space x > 0. At time t = 0, the surface of the semi-infinite body is exposed to a fluid at temperature T∞ . The convective heat transfer coefficient between the PCM and the fluid is h. Assume that the densities of the PCM for both phases are the same and that natural convection in the liquid phase is negligible. The temperature distribution in the liquid phase can be constructed based on the solution of heat conduction in a semiinfinite region subject to the boundary condition of the third kind. Find the transient location of the solid-liquid interface. 6.10. To solve heat conduction in a semi-infinite body with a convective boundary condition by using the integral approximate solution, the temperature distribution in the thermal penetration depth must be determined. If it can be assumed to be a second-order polynomial function, determine the unknown constants in the polynomial function. 6.11. Find thermal penetration depth for heat conduction in a semi-infinite body with convective boundary condition. The temperature distribution in the thermal penetration depth can be assumed to be a second-order polynomial function. 6.12. Solve for melting in a semi-infinite body with convective heating. The temperature distribution in the liquid phase can be assumed to be a second order polynomial function. 6.13. Solve for melting in a subcooled semi-infinite solid with constant wall temperature by using the integral approximate method. The temperature distribution in the liquid phase can be assumed to be a linear distribution. The temperature distribution in the solid phase can be assumed to be a second-order polynomial function. 6.14. The dimensionless governing equations of melting in a subcooled semiinfinite solid with constant heat flux heating were given in eqs. (6.145) – 524 Transport Phenomena in Multiphase Systems (6.151). Solve the problem using a semi-exact solution: exact solution in the liquid (see Example 6.1) and integral approximate solution in the solid (see Section 6.4.3). 6.15. You are given a one-dimensional heat sink to analyze. It works by melting a block of encased wax. The block of wax (PCM) is one centimeter thick. The heat sink operates by keeping one wall at 10 °C above the wax melt temperature. Assume that the solid PCM is uniformly at the melting point at time t = 0. The heat sink is used up when all of the PCM is liquid. The wax has the following properties: Tm = 0o C, cp = 1 J/kg-K, hsA = 100 J/kg, In calculating thermal capacity, what error would result if only the latent heat – and not also the sensible heat rise – is accounted for? 6.16. Since the surface temperature of an ablation material is very high, the surface is subjected to convective and radiation cooling. Modify the energy balance at the surface of the ablating material – eq. (6.175) – and obtain the ablation velocity with the convective and radiation cooling effects considered. 6.17. The surface heat flux in the ablation problem discussed in Section 6.4.4 was assumed as constant, and only steady-state solution was discussed. The heat flux during ablation is usually a function of time. A finite slab with thickness L and an initial temperature of Ti ( Ti < Tm ) is subject to a variable ′′ heat flux q0 (t ) = Aet at the surface x = 0. The surface x = L can be assumed to be adiabatic. Since the initial temperature of the slab is below the melting point, it must be preheated before ablation begins. Analyze the preheating and ablation problem using integral approximate method. It is assumed that the thermal penetration depth in the ablating material never reaches to x = L so the ablating material can be treated as a semi-infinite solid. The temperature distribution in the solid can be assumed as the second-degree polynomial function. 6.18. Derive the quasi-steady formulation for the temperature of the ablative ′′ material and substrate when a constant heat flux qw is applied to the surface of the ablative material, as shown in Fig. P6.3. Develop the governing equations, including the boundary conditions and initial conditions needed. The two materials in this problem are assumed to be homogenous, isotropic and in good contact with each other. The surface temperature is assumed to remain constant during the process at a temperature that is lower than the melting temperature of the substrate. The inner surface is assumed to be adiabatic. The temperature distribution is initially uniform and then exposed to a constant external heat flux on the outer surface of the ablative material. A phase change occurs at the outer surface and the new material is immediately removed upon formation. Find the ablation rate for the limiting case of infinite thickness of the ablation material at steady state. Chapter 6 Melting and Solidification 525 Ablative Material Substrate L Adiabatic ′′ qw T1 k1 cp,1 1 T q T2 k2 cp,2 2 y y=0 x=0 y=s Figure P6.3 x 6.19. A horizontal cylinder is embedded in solid paraffin with a temperature equal to the melting point of paraffin. From time t = 0 , a constant heat flux is applied to the inner surface of the tube so that melting occurs around the tube. Assume that the thickness of the tube is negligible. Solve the melting problem using the integral approximate method. 6.20. The physical model of melting in an eccentric annulus is shown in Figure P6.4. The PCM fills the eccentric annulus with outer radius ro, inner radius ri, and eccentricity e. The initial temperature of the PCM is assumed to be at its freezing temperature, Tm. At the beginning, t = 0, the temperatures of the outer and inner walls of the eccentric annulus are suddenly reduced to a temperature Tc, below the freezing temperature. The freezing process will occur simultaneously at the outer and inner walls of the eccentric annulus. It is assumed that the freezing at the inside of the outer wall and at the outside of the inner wall can be solved independently. Before the freezing fronts meet, the freezing volume is the sum of the volumes at the inside of the outer wall and the outside of the inner wall. After the freezing fronts meet, the freezing rate will decrease since the freezing process stopped at the locations where the freezing fronts met. Therefore, the freezing rate needs to be calculated for two different stages: before and after the freezing fronts meet. The radius of the inner and outer freezing fronts can be determined by an integral approximate solution. Obtain the rate of freezing as a function of time. 526 Transport Phenomena in Multiphase Systems ro PCM Tc e ri φ Tc r Figure P6.4 6.21. In a shell-and-tube latent heat thermal energy storage system, as shown in Fig. P6.5, the solid PCM at melting point fills in the shell side, while the transfer fluid flows inside the tube. It is assumed that (1) the inlet velocity of the transfer fluid is fully developed, but the heat transfer on the tube side is in the thermal entry region, (2) the quasi-steady assumption is applicable to convective heat transfer inside the tube, (3) the axial conductions in both the tube side and shell side are assumed to be negligible, and (4) the tube wall is assumed to be very thin so that its heat capacity and axial conduction can be neglected. Find the transient location of the melting front before it reaches the adiabatic wall of the thermal energy storage device by analyzing forced convective heat transfer inside the tube and melting on the shell side. Figure P6.5 6.22. A line heat source with intensity of q′ (W/m) is located at r = 0 in an infinite solid as shown in Fig. P6.6. The initial temperature of the solid, Ti, Chapter 6 Melting and Solidification 527 is below the melting point, Tm, so that this is a two-region melting problem. Develop the governing equations of the problem and obtain the transient location of the melting front. The temperature distributions in the liquid and solid phases can be constructed using an exponential function in the −z ∞e form of Ei[ − r 2 /(4α t )] , where Ei ( z ) = ³ dz . z z Line heat source Tm Ti s(t) r Figure P6.6 0 6.23. One-region melting in a semi-infinite body with constant surface temperature can be described by eqs. (6.42) – (6.45). In a quasi-stationary approximate solution, the derivative of the temperature in the liquid phase with respect to time is assumed to be zero ( ∂θ/∂τ = 0 ). Consequently, the energy equation (6.42) is simplified as d 2θ =0 dX 2 0 < X < S (τ ), τ > 0 Solve the melting problem using quasi-stationary approximate solution and compare your results with integral solution for Ste = 0.1, 0.2, and 0.4. 6.24. A very long cylindrical biological material with a radius of r0 is placed in a freezer at temperature T0 , which is assumed to be the surface temperature of the cylinder during freezing. Assume that the biological material can be treated as single-component PCM with a well defined melting point, Tm, and that its initial temperature is at Tm. The freezing process is so slow that heat transfer in a frozen layer can be regarded as a quasi-steady-state process. Estimate the time it takes to freeze the entire cylinder. 6.25. Solve melting around a hollow sphere with boundary condition of the first kind by using the quasi-stationary method. The initial temperature of the PCM is assumed to be equal to the melting point. 528 Transport Phenomena in Multiphase Systems 6.26. A liquid PCM with initial temperature Ti ( Ti > Tm ) enters a circular tube with diameter D and wall temperature Tw ( Tw < Tm , see Fig. P6.7). The liquid PCM is cooled in the tube and solidification begins at x = Ls . While the diameter of the liquid channel shrinks as a result of solidification, the Nusselt number is assumed to be a constant before and after solidification begins. Assuming that axial conduction in the liquid and solid PCM are negligible, and the solidification inside the tube is steady-state, obtain the length at which solidification begins and the thickness of the solid layer. The properties of the liquid and solid PCM can be assumed to be the same. Ls Tw δ Solid D Solid Figure P6.7 6.27. Show that the energy equation for the enthalpy model, eq. (6.211),satisfies the energy equations for both the liquid and solid phases, eqs. (6.20) and (6.22), as well as, energy balance, eq. (6.26). 6.28. In the enthalpy method, we obtained an equation for the numerical solution: h n +1 = j Δt Ti ui ρ ( Δx ) § ¨ ¨ 2© k j + 1 T jn+1 + k j − 1 T 2 2 n· ¸ j −1 ¸ ¹ + hn − j § ¨ ¨ ¨ ¨ ¨ ¨ © Δt ¨ k j + 1 + k j − 1 ¸ ¨ ¸ © 2 2 § · ρ ( Δx ) ¹ 2 T · ¸ n¸ ¸ j¸ ¸ ¸ ¹ j = 2,, N − 1 Prove that the above equation is stable only if the following equation is satisfied: max(α s ,α A )Δt ( Δx ) 2 ≤ 1 2 6.29. Phase change of a PCM occurs over a range of temperatures ( Tm − ΔT , Tm + ΔT ). A finite slab of this PCM with thickness L has an initial temperature Ti , which is below Tm − ΔT . At time t = 0, the temperature of the left surface suddenly increases to T0 > Tm + ΔT , while the right surface remains insulated. Give the governing equations using enthalpy method and design the solution procedure. Chapter 6 Melting and Solidification 529 6.30. You are given a one-dimensional melting problem that you need to solve with an explicit numerical scheme. The phase change material is nOctadecane. a) b) c) Given a time step of 0.1 seconds, find the minimum x necessary for stable numerical solution. Given a x of 1mm, find the minimum time step necessary for stable numerical solution. What is the requirement for numerical simulation using the implicit method? 6.31. Show that the enthalpy-temperarture relation in Fig. 6.18 can be represented analytically by eq. (6.237). 6.32. Most materials that are used as solid-to-liquid PCMs have poor thermal conductivity. Describe at least three ways of increasing thermal conductivity to melt the material located “far” from the heat source. 530 Transport Phenomena in Multiphase Systems