3
HEAT CONDUCTION
3.1 Introduction
Heat conduction is the process of heat transfer across a stationary medium, either solid or fluid. The distinction between conduction and convection is that there is no bulk motion of heat transfer media in conduction. Heat conduction is an atomic or molecular activity that transfers thermal energy from a region with higher temperature to a region with lower temperature. The mechanism of heat conduction is different for different substances. For an electrically nonconducting solid, conduction is attributed to atomic activity in the form of lattice vibrational waves (or phonons). The mechanism of conduction in an electrically conducting solid is a combination of lattice vibration and translational motion of electrons. Heat conduction in a liquid or gas is due to the random motion and interaction of the molecules. Since heat conduction always occurs in an incompressible substance and the velocity is zero in the conducting media, the energy equation (2.92) becomes ∂T ρcp = −∇ ⋅ q′′ + q′′′ (3.1) ∂t where q′′′ is the internal heat source (W/m3) that can be caused by electrical heating, chemical reaction heat, blood perfusion (for bioheat transfer), or coupling between electron and phonon (for ultrafast heat transfer in metal). For most engineering problems, the heat transfer rate is related to the temperature gradient by Fourier’s law: q′′ = −k ⋅ ∇T (3.2) where the negative sign on the right-hand side means that the heat flux vector, q′′ , is in the opposite direction of the temperature gradient, ∇T . For anisotropic material, the thermal conductivity is a tensor of the second order that has nine components: ⎡ k xx k xy k xz ⎤ ⎢ ⎥ (3.3) k = ⎢ k yx k yy k yz ⎥ ⎢ k zx k zy k zz ⎥ ⎣ ⎦ where k xy = k yx , k xz = k zx , and k yz = k zy according to the reciprocity relation derived from the Onsagar’s principle of thermodynamics of irreversible processes (Ozisk, 1993). The dot product of a second order tensor, k, and a vector, ∇T ,
Chapter 3 Heat Conduction 213
results in a vector, q′′ . The energy equation for an anisotropic material can be obtained by substituting eq. (3.2) into eq. (3.1), i.e., ∂T ρcp = ∇ ⋅ (k ⋅ ∇T ) + q′′′ (3.4) ∂t In a Cartesian coordinate system, eq. (3.4) becomes ∂T ∂ ⎛ ∂T ∂T ∂T ⎞ ∂ ⎛ ∂T ∂T ∂T ⎞ = ⎜ k xx + k xy + k xz + k yy + k yz ρcp ⎟ + ⎜ k yx ⎟ ∂t ∂x ⎝ ∂x ∂y ∂z ⎠ ∂y ⎝ ∂x ∂y ∂z ⎠
∂ ⎛ ∂T ∂T ∂T ⎞ + k zy + k zz (3.5) ⎜ k zx ⎟ + q′′′ ∂z ⎝ ∂x ∂y ∂z ⎠ For the case that all non-diagonal components in eq. (3.3) are absent, the anisotropic material becomes orthotropic material and the energy equation for heat conduction becomes ∂T ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ρcp = ⎜ k xx (3.6) ⎟ + ⎜ k zz ⎟ + ⎜ k yy ⎟ + q′′′ ∂t ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ Heat conduction in an orthotropic material is much easier to compute than that for the anisotropic material because cross-derivatives of the space variables are absent. If the thermal conductivities in the three directions are the same ( k xx = k yy = k zz = k ), the material becomes isotropic and the thermal conductivity +
becomes a scalar. Fourier’s law for an isotropic material is q′′ = −k ∇T (3.7) and the energy equation (3.1) and (3.6) are simplified as ∂T ρcp = ∇ ⋅ (k ∇T ) + q′′′ (3.8) ∂t ∂T ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ρcp = ⎜k (3.9) ⎟+ ⎜k ⎟ + ⎜k ⎟ + q′′′ ∂t ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ While eqs. (3.5), (3.6) and (3.9) are given in the Cartesian coordinate system, the corresponding equations in the cylindrical and spherical coordinate systems can be obtained by using the ∇ operator in these coordinate systems (see Appendix G, and Problem 3.2 and 3.3). Thermal conductivity is a very important thermophysical property of a material because the heat flux increases with increasing thermal conductivity under a prescribed temperature gradient. Since the interaction between molecules in the solid is strongest among the three phases, the thermal conductivity for the solid is usually higher than that in the liquid. For the same reason, the thermal conductivity for the gas is usually lowest among the three phases. The thermal conductivity for the metal is higher than that for a non-conducting solid because of the additional carrier – free electrons – of heat in metal. The typical ranges of thermal conductivities for various materials are shown in Fig. 3.1. It can be seen that the thermal conductivities for different materials can cross four orders of
214 Advanced Heat and Mass Transfer
Figure 3.1 Typical range of thermal conductivities (Ozisik, 1993)
magnitude. Thermal conductivity is usually a function of temperature. The thermal conductivity for metal increases with increasing temperature, but it decreases with increasing temperature for gases. The thermal conductivities for some materials can be found in Appendix B-D. The energy equation or the boundary condition is said to be linear if there are no products of unknown variable (temperature) or its derivative in the energy equation or boundary condition. For example, eq. (3.8) or (3.9) is nonlinear because product of k – a function of temperature – and the derivative of temperature appeared. When the thermal conductivity is independent from temperature, eq. (3.8) and (3.9) become 1 ∂T q′′′ (3.10) = ∇ 2T + α ∂t k 1 ∂T ∂ 2T ∂ 2T ∂ 2T q′′′ (3.11) = + + + α ∂t ∂x 2 ∂y 2 ∂z 2 k where k α= (3.12) ρcp is the thermal diffusivity – another thermophysical property. Equations (3.10) and (3.11) are linear, because they do not contain products of T or its derivative. Furthermore, a linear energy equation or boundary condition is said to be homogeneous if the temperature T can be replaced by CT, where C is a nonzero constant, and the equation is still satisfied. Thus, eq. (3.10) or (3.11) is nonhomogeneous due to the presence of the internal heat generation term. They become linear homogeneous when there is no internal heat generation.
Chapter 3 Heat Conduction 215
This chapter presents a number of prototypical problems of heat conduction based on different models; its goal is to establish the physics of heat conduction and to demonstrate the variety of tools available for its solution. Sections 3.2 and 3.3 respectively discuss the classical heat conduction under steady and unsteady conditions, while numerical solution of heat conduction will be addressed in Section 3.4. Heat conduction with phase change – melting and solidification – will be discussed in detail in Section 3.5. This chapter closes with a discussion on microscale heat conduction (Section 3.6), which includes pure conduction as well as melting and solidification.
3.2 Steady State Heat Conduction
Under steady state conditions, the energy equation for classical heat conduction, eq. (3.8) is reduced to ∇ ⋅ (k ∇T ) + q′′′ = 0 (3.13) If the thermal conductivity is independent from the temperature, eq. (3.13) is further simplified to q′′′ ∇ 2T + =0 (3.14) k Analytical solutions of one-dimensional and multidimensional heat conduction under various physical conditions and geometrical configurations will be discussed.
3.2.1 One Dimensional Systems
One-dimensional heat conduction for the cases without heat generation, with heat generation, and heat conduction from extended surfaces will be discussed in this subsection.
Without Heat Generation
One-dimensional heat conduction can occur in different geometric configurations. Figure 3.2 shows heat conduction in a plane wall and in a hollow cylinder or sphere. The energy equation (3.13) in different geometric configurations can be expressed as (see Appendix G): d ⎛ dT ⎞ (3.15) plane wall (x1 < x < x2 ) ⎜k ⎟ = 0, dx ⎝ dx ⎠ 1 d ⎛ dT ⎞ (3.16) ⎜ kr ⎟ = 0, in hollow cylinder (r1 < r < r2 ) r dr ⎝ dr ⎠ 1 d ⎛ 2 dT ⎞ (3.17) ⎜ kr ⎟ = 0, in hollow sphere (r1 < r < r2 ) r 2 dr ⎝ dr ⎠
216 Advanced Heat and Mass Transfer
T1 T2 T1 T2 0 r1 r2 r
0
x1
x2 (a) Plane wall
x (b) Hollow cylinder or sphere
Figure 3.2 Heat conduction in different geometric configurations
Equations (3.15) – (3.17) can be represented using the following general form: d⎛ dT ⎞ (3.18) ⎜ kA( s ) ⎟ = 0, s1 < s < s2 ds ⎝ ds ⎠ where s is space variable and A(s) can be viewed as the heat transfer area. The heat transfer area in a plane wall, hollow cylinder and sphere are A( s) = LW , 2πsL, and 2πs2, respectively (where W is the width of the plane wall, and L is length of the plane wall or cylinder). Figure 3.3 represents the general onedimensional heat conduction problem. Equation (3.18) can also be rewritten as dq = 0, s1 < s < s2 (3.19) ds or dT q = − kA( s ) = const, s1 < s < s2 (3.20) ds which means the heat transfer rate at different cross-sections is a constant. It should be pointed out that the heat flux q′′( s ) = − kdT / ds is not a constant unless the cross-sectional area A(s) is constant. If the thermal conductivity is independent from temperature, the temperature distribution can be obtained by integrating eq. (3.18) twice: s 1 T =a ds + b (3.21) s1 A( s ) where a and b are unspecified integral constants. If the temperatures of both inner and outer surfaces are specified, i.e., T = T1 , s = s1 (3.22) T = T2 , s = s2 (3.23)
∫
Chapter 3 Heat Conduction 217
Figure 3.3 Generalized 1-D heat conduction problem (Arpaci, 1966)
The integral constants a and b in eq. (3.21) can be determined from eqs. (3.22) and (3.23), and the temperature becomes s T −T 1 T= s 2 1 ds + T1 (3.24) 2 1 s1 A( s ) ds s1 A( s ) The rate of heat transfer can be obtained using Fourier’s law: T −T dT q = − kA( s ) = 1s 2 (3.25) ds 1 2 ds k s1 A( s ) which can be rewritten as T −T (3.26) q= 1 2 Rc where 1 s2 ds Rc = (3.27) k s1 A( s ) is the thermal resistance for heat conduction. Equation (3.26) is similar to Ohm’s law V −V I= 1 2 (3.28) R where V1 and V2 are, respectively, the electrical potential at two ends of the electrical resistance, R, and I is the current through the electrical resistance (see Fig. 3.4).
∫
∫
∫
∫
218 Advanced Heat and Mass Transfer
I V1 R q T1 Rc T2 V2
Figure 3.4 Analogy between conductions of electricity and heat
If the boundary conditions at the inner and outer surfaces are convective conditions, i.e., dT −k = hi (Ti − T ), s = s1 (3.29) ds dT −k = ho (T − To ), s = s2 (3.30) ds where hi and ho are heat transfer coefficients at the inner and outer surfaces, respectively. The temperatures of the fluids that are in contact with the inner and outer surfaces are Ti and To , respectively. If the temperatures of the inner and outer walls are represented by T1 and T2, (both of them are unknown) eq. (3.26) is still valid. By multiplying eqs. (3.29) and (3.30) by the heat transfer areas at the inner and outer surfaces, A1 and A2, one obtains dT q = − kA( s1 ) = hi A1 (Ti − T1 ), s = s1 (3.31) ds dT q = − kA( s2 ) = ho A2 (T2 − To ), s = s2 (3.32) ds Under steady state condition, the heat transfer rate obtained from eqs. (3.26), (3.31) and (3.32) are identical. By eliminating T1 and T2 from these equations, the following expression for heat transfer rate is obtained: Ti − To (3.33) q= Ri + Rc + Ro where 1 Ri = (3.34) hi A1 1 Ro = (3.35) ho A2
Chapter 3 Heat Conduction 219
are convective thermal resistances between inner fluid and inner wall and between outer fluid and outer wall, respectively. Equation (3.33) can be viewed as a system with three thermal resistances connected in series. Equation (3.33) can also be written as q = U o A2 (Ti − To ) (3.36) where Uo is the overall- coefficient of heat transfer based on the outer surface area. Comparing eqs. (3.33) and (3.36) yields 1 = Ri + Rc + Ro (3.37) U o A2 or A A s2 ds 1 1 = 2+ 2 + (3.38) s1 A( s ) U o hi A1 k ho For the heat conduction in three different coordinates as shown in eqs. (3.15) – (3.17), the overall coefficients for heat transfer are 1 1L1 plane wall (x2 − x1 = L) =++ (3.39) U o hi k ho
∫
r r ⎛r ⎞ 1 1 = 2 + 2 ln ⎜ 2 ⎟ + , hollow cylinder U o hi r1 k ⎝ r1 ⎠ ho
(3.40)
⎞1 r2 r ⎛r 1 = 2 2 + 2 ⎜ 2 − 1⎟ + , hollow sphere (3.41) U o hi r1 k ⎝ r1 ⎠ ho If the conducting wall shown in Fig. 3.2 has multiple layers and each layer has different thermal conductivity, there will be multiple conduction thermal resistances between two convection thermal resistances. If the number of layers is represented by N, the overall coefficient of heat transfer will be expressed as N A 1 1 si+1 ds 1 = N +1 + AN +1 + (3.42) si U o hi A1 A( s ) ho i =1 ki where A1 and AN+1 are the areas of the inner and outer surfaces, respectively; and ki is the thermal conductivity of the ith layer ( si < s < si +1 ). The heat transfer rate becomes q = U o AN +1 (Ti − To ) (3.43) or Ti − To (3.44) q= N 1 1 si+1 ds 1 + + hi A1 i =1 ki si A( s ) ho AN +1
∑∫
∑∫
Example 3.1 Steam flows through a 10-m-long carbon steel (k = 50 W/m-K) pipe with an outer diameter of 100 mm and a thickness of 5mm. The paper faced glass fiber with thermal conductivity of 0.038 W/m-K is wrapped on the steam pipe as insulation material (see Fig. 3.5). The
220 Advanced Heat and Mass Transfer
T1 T2 T3
0 r1 r2 r3
Figure 3.5 Heat conduction in a steam pipe
steam temperature is 100 oC and the inner surface temperature of the steam pipe is assumed to be equal to the steam temperature. The outer surface temperature of the insulation material is 20 oC. The heat loss per unit length of the pipe must be less than 1000 W. Find the required thickness of the insulation material and the temperature at the interface between the carbon steel and glass fiber.
Solution: The thermal conductivities of the carbon steel and paper faced glass fiber are k1 = 50 W/m-K and k2 = 0.038 W/m-K , respectively. The inner and outer diameters of the carbon steel tube are r1 = 0.045 m and r2 = 0.050 m . Since the temperatures at the inner and outer surfaces are known, the convective thermal resistances in eq. (3.44) do not need to be included and the heat transfer rate becomes T1 − T3 q= ⎛r ⎞ ⎛r ⎞ 1 1 ln ⎜ 2 ⎟ + ln ⎜ 3 ⎟ 2π Lk1 ⎝ r1 ⎠ 2π Lk2 ⎝ r2 ⎠ which can be rearranged as ⎛ r ⎞ 2π k2 L(T1 − T3 ) k2 ⎛ r2 ⎞ − ln ⎜ ⎟ ln ⎜ 3 ⎟ = q k1 ⎝ r1 ⎠ ⎝ r2 ⎠
= i.e.,
2π × 0.038 × 10 × (100 − 20) 0.038 ⎛ 0.050 ⎞ − ln ⎜ ⎟ = 0.19093 1000 50 ⎝ 0.045 ⎠
r3 = 1.2104 r2 Thus, the outer radius of the glass fiber is r3 = 0.0605m , and the thickness of the glass fiber is δ = r3 − r2 = 0.0105 m = 10.5 mm . If we take the steel tube only, the heat transfer rate is
Chapter 3 Heat Conduction 221
q=
T1 − T2 ⎛r ⎞ 1 ln ⎜ 2 ⎟ 2π Lk1 ⎝ r1 ⎠
which can be rearranged to ⎛r ⎞ q 1000 ⎛ 0.050 ⎞ o T2 = T1 − ln ⎜ 2 ⎟ = 100 − ln ⎜ ⎟ = 99.97 C 2π Lk1 ⎝ r1 ⎠ 2π × 10 × 50 ⎝ 0.045 ⎠ In arriving at eq. (3.42), it is assumed that different layers are in perfect contact so that the temperatures across the interfaces between different layers are continuous. When two rough surfaces are in contact, there may be a temperature drop, ΔT, across the interface between different materials (see Fig. 3.6). The heat flux and this temperature drop is related by ΔT (3.45) q′′ = i ′′ Rct ,i where Rc′′t is contact thermal resistance (m2-K/W). Heat transfer across an imperfect contact interface is due to combined effects of conduction at the actual contact area, convection of the air entrapped in the gap, and radiation between the two surfaces that are not in direct contact. The contact thermal resistance can be reduced by applying pressure in the direction perpendicular to the interface, or applying conduction grease at the interface. More information about contact thermal resistance can be found in Fletcher (1988). When contact thermal resistances are present, eq. (3.44) can be modified as Ti − To (3.46) q= N N ′′ Rtc ,i 1 1 si+1 ds 1 + + + hi A1 i =1 ki si A( s ) i = 2 Ai ho AN +1
∑∫
∑
Ti
si
q′′
Δ Ti
Ti+1 ith layer
Figure 3.6 Contact thermal resistance
i+1th layer
222 Advanced Heat and Mass Transfer
With Internal Heat Generation
The 1-D steady-state heat conductions that we discussed so far are limited to the case without internal heat generation. Heat conduction with internal heat generation can be encountered in many applications such as electrical heating, chemical reaction, or nuclear reaction in the conduction medium. Let us consider the generalized 1-D heat conduction problem shown in Fig. 3.3. With uniform internal heat source, q′′′, the energy equation is (see Problem 3.5) 1 d⎛ dT ⎞ q′′′ (3.47) = 0, s1 < s < s2 ⎜ A( s ) ⎟+ A( s ) ds ⎝ ds ⎠ k where the thermal conductivity is assumed to be independent from temperature. Equation (3.47) is subject to boundary conditions specified by eqs. (3.22) and (3.23). Multiplying eq. (3.47) by A(s) and integrating the resultant equation in the interval of (s1, s) yields dT q′′′ s dT =− (3.48) A( s ) A( s )ds + A1 ds k s1 ds s = s1
∫
Dividing eq. (3.48) by A(s) and integrating the resultant equation in the interval of (s1, s), we have s s1 ⎤ q′′′ s ⎡ 1 dT T =− ds + T1 (3.49) ⎢ A( s ) s A( s )ds ⎥ ds + A1 ds s1 A( s ) k s1 ⎣ 1 ⎦ s = s1
∫
∫
∫
where the temperature gradient at s = s1 is still unknown. Equation (3.49) already satisfies eq. (3.22) but not eq. (3.23). Substituting eq. (3.49) into eq. (3.23), one obtains s s2 ⎤ q′′′ s2 ⎡ 1 dT 1 T2 = − A( s ) ds ⎥ ds + A1 ds + T1 ⎢ A( s ) s k s1 ⎣ ds s = s1 s1 A( s ) 1 ⎦
∫
∫
∫
which can be rearranged to get A1 dT ds (T2 − T1 ) + =
s = s1
q′′′ k
∫
s2
s1
∫
s2
s1
⎡1 ⎢ A( s ) ⎣ 1 ds A( s )
∫
⎤ A( s )ds ⎥ ds s1 ⎦
s
Therefore, the temperature profile becomes s ⎤ q′′′ s ⎡ 1 T =− ⎢ A( s ) s A( s )ds ⎥ ds k s1 ⎣ 1 ⎦
∫
∫
s ⎡1 ⎤ ⎢ A( s ) s A( s )ds ⎥ ds s1 1 ⎣ ⎦ + s2 1 ds s1 A( s ) The heat transfer rate can be obtained by Fourier’s law
(T2 − T1 ) +
q′′′ k
∫
s2
∫
∫
∫
s
s1
1 ds + T1 A( s )
(3.50)
Chapter 3 Heat Conduction 223
dT ds Substituting eq. (3.48) into eq. (3.51), the heat transfer rate becomes q = − kA( s ) q = q′′′ where
(3.51)
∫
s
s1
A( s )ds + q( s1 )
(3.52)
s ⎡1 ⎤ ⎢ A( s ) s A( s ) ds ⎥ ds s1 dT 1 ⎣ ⎦ = (3.53) q( s1 ) = − kA1 s2 1 ds s = s1 ds s1 A( s ) is the heat transfer rate at s = s1 . It is evident from eq. (3.52) that the heat transfer rate is no longer independent from s when there is internal heat generation. It should be pointed out that eqs. (3.50) and (3.52) are valid for any coordinate system. For a Cartesian coordinate system with the origin of the coordinate on the left surface (i.e., x1 = 0 in Fig. 3.2), A( s ) = A is a constant and eqs. (3.50) and (3.52) become (T − T ) q′′′ (3.54) T ( x) = [ ( L − x) x ] + 2 1 x + T1 2k L L ⎞ kA(T1 − T2 ) ⎛ (3.55) q = q′′′AL ⎜ x − ⎟ + 2⎠ L ⎝ For the heat conduction in a cylindrical and spherical coordinate system, the general solution, eqs. (3.50) and (3.52), can be simplified by considering the variation of conduction area (see Problem 3.7 and 3.8).
k (T1 − T2 ) − q′′′
∫
s2
∫
∫
Heat Transfer from Extended Surfaces
The total thermal resistance in eq. (3.37) includes two convective thermal resistances and one conduction thermal resistance. For the cases that one of the convection thermal resistances is dominant (i.e., significantly greater than the conduction thermal resistance and the other convective thermal resistance), one can increase the heat transfer coefficient (hi or ho) or the heat transfer area (A1 or A2) to enhance the overall heat transfer. Since increasing the heat transfer coefficient is constrained by the type of working fluid and the power required driving the flow, increasing the heat transfer area becomes a natural choice. The increase of surface area can be done by using fins that extend into the fluid. Figure 3.7 shows some examples of fin configurations. It can be seen that the fins can have either uniform or variable cross-sectional areas. For cases where the fluid is cooling the fin, the fin temperature is the highest at the base ( x = 0 ) and decreases with increasing x as convection takes place on the fin surface. The degree of heat transfer enhancement can be maximized by minimizing the temperature variation, which can be achieved by using fin materials with large
224 Advanced Heat and Mass Transfer
Figure 3.7 Fin configurations: (a) straight fin of uniform cross-section on plane wall, (b) straight fin of uniform cross-section on circular tube, (c) annular fin, and (d) straight pin fin (Holman, 2002)
T∞, h
0 T0 x x+dx
A(x)
x
Figure 3.8 Heat conduction in a fin with variable cross-section
thermal conductivity. The solution of heat transfer from extended surfaces (fins) will provide (a) temperature distribution in the fin, and (b) total heat transfer from the finned surface. Figure 3.8 shows the generalized physical model for heat transfer form an extended surface. It is a steady-state conduction problem in variable crosssectional area. It is different from the preceding subsection in that convection occurs on the extended surface. The following assumptions are made to simplify the problem: 1. Heat transfer is steady-state, 2. The thermal conductivity of the fin is independent from the temperature, 3. There is no internal heat generation in the fin, 4. Both convective heat transfer coefficient and fluid temperature are constants, 5. The heat transfer is one-dimensional, i.e., the temperature is uniform in the same cross-sectional area. Heat transfer in a fin can be modeled as a steady-state one-dimensional conduction with internal heat generation, described by eq. (3.47), if the convection from the extended surface is treated as an equivalent internal heat source. For the differential control volume ( dV = Adx ) shown in Fig. 3.8, the convective heat transfer from the side is
Chapter 3 Heat Conduction 225
dqconv = hPdx [T∞ − T ( x) ] where P is the perimeter of the fin. The equivalent intensity of the internal heat source in the control volume due to the convection on the extended surface is dq hP q′′′ = conv = [T∞ − T ( x)] dV A Therefore, the energy equation becomes 1 d ⎛ dT ⎞ hP (3.56) (T − T∞ ) = 0 ⎜A ⎟− A dx ⎝ dx ⎠ Ak or d 2T ⎛ 1 dA ⎞ dT hP +⎜ − (3.57) (T − T∞ ) = 0 ⎟ dx 2 ⎝ A dx ⎠ dx Ak which requires two boundary conditions. At the base of the fin, the temperature is T = T0 , x = 0 (3.58) There are commonly three kinds of different boundary conditions at the tip of the fin: 1. When the fin is sufficiently long, the temperature at the tip can be assumed to be equal to the fluid temperature, T = T∞ , x = L (3.59) 2. It is assumed that the heat transfer at the tip of the fin is negligible, dT = 0, x = L (3.60) dx 3. When the fin is not long enough, the tip temperature is higher than the fluid temperature and convection occurs at the fin tip, dT −k = h(T − T∞ ), x = L (3.61) dx Heat transfer from an extended surface is a non-homogeneous problem because eq. (3.56) or (3.57) is not homogonous. If convection at the fin tip is considered, eq. (3.61) is also non-homogeneous. Defining excess temperature ϑ ( x) = T ( x) − T∞ Equations (3.57) – (3.61) become d 2ϑ ⎛ 1 dA ⎞ dϑ hP (3.62) +⎜ − ϑ =0 ⎟ dx 2 ⎝ A dx ⎠ dx Ak ϑ = ϑ0 , x = 0 (3.63) (3.64) ϑ = 0, x = L dϑ = 0, x = L (3.65) dx dϑ −k = hϑ , x = L (3.66) dx which are all homogeneous. The total amount of heat transfer from a fin can be obtained by using Fourier’s law at the base of the fin
226 Advanced Heat and Mass Transfer
⎛ dϑ ⎞ (3.67) q f = −k ⎜ A ⎟ ⎝ dx ⎠ x = 0 or accounting for convective heat transfer throughout the fin surface by qf = h
∫
L
0
Pϑ dx
(3.68)
Under steady state condition, the heat conduction from the base of the fin obtained from eq. (3.67) is identical to the total convective heat transfer obtained from eq. (3.68), however, eq. (3.67) is much easier to apply than eq. (3.68). The performance of a fin is usually measured by the fin efficiency, η f , which is defined as the ratio of the actual heat transfer rate of a fin to the heat which would be transferred if the entire fin area is at base temperature, i.e., q qf (3.69) ηf = f = q f ,max hϑ L Pdx 0
∫
0
While the general solution of the heat conduction in a fin with arbitrarily variable cross-sectional area cannot be obtained, analytical solutions for some special cases are possible. For fins with uniform cross-sections, the second term in eq. (3.62) is dropped and the energy equation becomes d 2ϑ − m 2ϑ = 0 (3.70) dx 2 where hP (3.71) m2 = Ak The general solution of eq. (3.70) is ϑ = C1e mx + C2 e − mx (3.72) where C1 and C2 are integral constants that need to be determined by boundary conditions at the base and the tip of the fin. The temperature distributions, heat transfer rate, and fin efficiency for the three boundary conditions represented by eqs. (3.59) – (3.61) are: 1. For the case that the fin’s top temperature equals the fluid temperature:
ϑ = e − mx ϑ0
(3.73) (3.74) (3.75)
q f = hPkAϑ0
1 mL 2. For the case that the fin tip is adiabatic ϑ cosh[m( L − x)] = cosh( mL) ϑ0
ηf =
(3.76)
(3.77)
q f = hPkAϑ0 tanh(mL)
Chapter 3 Heat Conduction 227
tanh( mL) (3.78) mL 3. For the case of a convective fin tip ϑ cosh[m( L − x)] + (h / mk )sinh[m( L − x)] = (3.79) ϑ0 cosh(mL) + (h / mk )sinh(mL) sinh(mL) + (h / mk ) cosh(mL) (3.80) q f = hPkAϑ0 cosh(mL ) + (h / mk )sinh(mL) 1 sinh(mL) + (h / mk )cosh( mL) ηf = (3.81) mL cosh(mL) + (h / mk )sinh( mL) Among the three possible boundary conditions, the first case is valid only if the fin is sufficiently long. For this reason it cannot be applied in most cases. While the convective fin tip is most accurate for many applications, its results are very complicated. Practically, one can use the adiabatic fin tip for the case with convection on the fin tip by using a corrected fin length, Lc, to replace L in eqs. (3.76) – (3.78). The corrected length is Lc = L + t / 2 for a rectangular fin [see Fig. 3.7 (a) and (b), t is fin thickness], and Lc = L + D / 4 for a pin fin [see Fig. 3.7(d), D is pin fin diameter]. This treatment is based on the assumption that the amount of heat transfer from the fin tip can be considered by an equivalent amount of heat transfer from the increased fin surface area. The above solution is valid for the case that the cross-sectional area is uniform. If the cross-sectional area is not uniform [e.g., Fig. 3.7(c)], the second term in eq. (3.62) is retained and the solution becomes more complicated. Figure 3.9 shows the physical model of heat transfer in an annular fin with uniform thickness. The perimeter and the cross-sectional area of the annular fin are P = 2π r A = 2π rt
ηf =
t
r1
r2
Figure 3.9 Heat transfer in an annular fin with uniform thickness
228 Advanced Heat and Mass Transfer
The energy equation (3.62) becomes d 2ϑ 1 dϑ − − m 2ϑ = 0, r1 < r < r2 (3.82) dr 2 r dr where 2h m= (3.83) kt The boundary conditions for eq. (3.82) are θ = ϑ0 , r = r1 (3.84) dϑ = 0, r = r2 (3.85) dr Equation (3.82) is a modified Bessel’s equation of zero order and has the following general solution ϑ (r ) = C1 I 0 (mr ) + C2 K 0 (mr ) (3.86) where I0 and K0 are zero order modified Bessel function of the first and second kinds, respectively (See Appendix G). Substituting eq. (3.86) into eqs. (3.84) and (3.85), the two integral constants in eq. (3.86) can be determined and the temperature distribution in the annular fin becomes K (mr ) I (mr ) + I1 (mr2 ) K 0 (mr ) θ (r ) = 1 2 0 (3.87) K1 (mr2 ) I 0 (mr1 ) + I1 (mr2 ) K 0 (mr1 ) The heat transfer rate can be determined using Fourier’s law at the base ( r = r1 ), i.e., q f = − k 2π r1t dϑ dx (3.88)
r = r1
Substituting eq. (3.87) into eq. (3.88) yields I (mr ) K (mr ) − K1 (mr2 ) I1 (mr1 ) q f = k 2π r1tϑ0 m 1 2 1 1 (3.89) I1 (mr2 ) K 0 ( mr1 ) + K1 (mr2 ) I 0 (mr1 ) If the entire fin surface is at the base temperature, the heat transfer amount is q f ,max = 2π (r22 − r12 )hϑ0 (3.90)
The fin efficiency for a straight annular fin is q I1 (mr2 ) K1 (mr1 ) − K1 (mr2 ) I1 ( mr1 ) 2r1 ηf = f = (3.91) 2 2 q f ,max m(r2 − r1 ) I1 (mr2 ) K 0 ( mr1 ) + K1 (mr2 ) I 0 (mr1 ) The fin efficiencies for other configurations can be found in Incropera et al. (2006).
Bioheat equation
Heat transfer in living tissue can be considered as heat conduction with internal heat generation from two sources: metabolic heat generation and blood perfusion. The former results from a series of chemical reactions that occur in the
Chapter 3 Heat Conduction 229
living cells, while the latter is due to energy exchange between the living tissue and blood flowing in the small capillaries. If the heat transfer in a tissue can be considered as one-dimensional and the thermal conductivity of the tissue is assumed to be independent from temperature, the energy equation for living tissue can be obtained by modifying eq. (3.47), i.e., ′′′ ′′′ 1 d⎛ dT ⎞ qm + qb (3.92) = 0, s1 < s < s2 ⎜ A( s ) ⎟+ A( s ) ds ⎝ ds ⎠ k ′′′ ′ where qm and qb′′ are the rate of metabolic heat generation and blood perfusion, respectively. Equation (3.92) is generalized for Cartesian ( s = x, A( s ) = const ), cylindrical ( s = r , A( s) = r ), and spherical ( s = r , A( s ) = r 2 ) coordinate systems. The internal heat generation due to blood perfusion can be obtained by ′ qb′′ = ωρb c pb (Ta − T ) (3.93) where ω is the blood perfusion rate – defined as volume flow rate of blood flow per unit volume of the tissue, ρb is the density of blood, c pb is the specific heat of blood under constant pressure, and Ta is the blood temperature. Depending on whether the blood temperature is higher than the tissue temperature, the internal heat generation obtained from eq. (3.93) can be either positive ( Ta > T ) or negative ( Ta < T ). Substituting eq. (3.93) into eq. (3.92), the energy equation for the tissue becomes ′′′ 1 d⎛ dT ⎞ qm + ωρb c pb (Ta − T ) = 0, s1 < s < s2 ⎜ A( s ) ⎟+ A( s ) ds ⎝ ds ⎠ k which is nonhomogeneous due to the presence of metabolic heating and blood perfusion. The above equation can be rearranged into ′′′ 1 d⎛ dT ⎞ ωρb c pb [T − Ta − qm /(ωρb c pb )] = 0, s1 < s < s2 ⎜ A( s ) ⎟− A( s ) ds ⎝ ds ⎠ k (3.94) If the blood temperature and the internal heat generation due to metabolic heating are constant, eq. (3.94) can be homogenized by introducing the following excess temperature ′′′ (3.95) ϑ = T − Ta − qm /(ωρb c pb ) i.e., 1 d⎛ dϑ ⎞ 2 ⎜ A( s ) ⎟ − m ϑ = 0, A( s ) ds ⎝ ds ⎠ where m2 = s1 < s < s2 (3.96)
ωρb c pb
k
(3.97)
230 Advanced Heat and Mass Transfer
Equation (3.96) is identical to eq. (3.62) and therefore, the solution for heat conduction from an extended surface can be readily applied to the solution of the bioheat equation. For different parts of the body, the generalized bioheat equation (3.96) can be changed to different forms. For the tissue with thickness much less than its curvature radius, the tissue can be treated as a plane wall and eq. (3.96) becomes d 2ϑ − m 2ϑ = 0, 0 < x < L (3.98) dx 2 If the effects of curvature need to be considered, eq. (3.96) can be rewritten for cylindrical or spherical coordinate systems, i.e., d 2ϑ 1 dϑ + − m 2ϑ = 0, r1 < r < r2 (Cylindrical) (3.99) dr 2 r dr d 2ϑ 2 dϑ + − m 2ϑ = 0, r1 < r < r2 (Spherical) (3.100) dr 2 r dr The bioheat equation in different coordinate systems can be solved with appropriate boundary conditions.
Example 3.2 The cylindrical biological tissue with a radius of ro, length L, and an initial temperature Ti is immersed into a fluid with a temperature of T∞. The convective heat transfer coefficient between the tissue and fluid is h. Find the temperature distribution in the cylindrical tissue and heat loss at the surface of the tissue. Solution: Assuming the heat conduction in the cylinder is onedimensional along the radial direction, the energy equation for this problem is d 2ϑ 1 dϑ + − m 2ϑ = 0, 0 < r < ro (3.101) dr 2 r dr and the boundary conditions for the biological tissue are dT = 0, r = 0 dr dT −k = h(T − T∞ ), r = ro dr which can be written in term of excess temperature defined in eq. (3.95), i.e., dϑ = 0, r = 0 (3.102) dr dϑ −k = h(ϑ − ϑ∞ ), r = ro (3.103) dr where ′′′ (3.104) ϑ∞ = T∞ − Ta − qm /(ωρb c pb )
Chapter 3 Heat Conduction 231
It is apparent that eq. (3.103) is not homogeneous even though it is in term of excess temperature. While this is different from heat conduction in an annular fin, eq. (3.101) is still a modified Bessel’s equation of zero order and has the following general solution ϑ (r ) = C1 I 0 (mr ) + C2 K 0 (mr ) (3.105) where C1 and C2 are unspecified constants that need to be determined by the boundary conditions. The gradient of excess temperature is (3.106) ϑ ′(r ) = C1mI1 (mr ) − C2 mK1 (mr ) Substituting eq. (3.106) into eq. (3.102), one obtains C1mI1 (0) − C2 mK1 (0) = 0 The characteristics of Bessel functions give us I1 (0) = 0 and K1 (0) = ∞ , thus, C2 must be zero and the excess temperature and its gradient become ϑ (r ) = C1 I 0 (mr ) and ϑ ′(r ) = C1mI1 (mr ) (3.107) Substituting eq. (3.107) into eq. (3.103) yields −kC1mI1 ( mro ) = h[C1 I 0 (mro ) − ϑ∞ ] which can be rearranged to obtain hϑ∞ C1 = hI 0 (mro )+kmI1 (mro ) The excess temperature thus becomes hϑ∞ I 0 (mr ) ϑ (r ) = (3.108) hI 0 (mro )+kmI1 (mro ) The heat loss from the surface of the biological tissue is then 2π ro Lhϑ∞ mI1 (mro ) (3.109) q = −2π ro Lkϑ ′(ro ) = − hI 0 (mro )+kmI1 (mro )
3.2.2 Multidimensional Systems
The heat transfer problems discussed in the preceding subsection are for the case that the temperature varies in one dimension only. For most applications, this is oversimplified and it is necessary to consider heat conduction in multidimensional systems. Analytical solution of multidimensional heat conduction using different methods will be discussed in this subsection. The numerical solution to be discussed in Section 3.4 can also be employed to solve steady-state heat conduction problems.
Two-Dimensional Heat Conduction without Internal Heat Generation
The problem under consideration is shown in Fig. 3.10. The left, right and bottom of the two-dimensional domain are kept at T1 while the top boundary has
232 Advanced Heat and Mass Transfer
y H T=f(x)
T=T1
T=T1
0
T=T1
L
x
Figure 3.10 Two-dimensional steady-state heat conduction
variable temperature f(x). It is assumed that the thermal conductivity is independent from the temperature. The energy equation for this problem is ∂ 2T ∂ 2T + = 0, 0 < x < L, 0 < y < H (3.110) ∂x 2 ∂y 2 with the following boundary conditions: T = T1 , x = 0, L, 0 < y < H (3.111) T = T1 , y = 0, 0 < x < L (3.112) T = f ( x), y = H , 0 < x < L (3.113) Although eq. (3.110) is linear homogeneous, eqs. (3.111) – (3.113) are not homogeneous. Introducing excess temperature ϑ = T − T1 eqs. (3.110) – (3.113) become ∂ 2ϑ ∂ 2ϑ + = 0, 0 < x < L, 0 < y < H (3.114) ∂x 2 ∂y 2 ϑ = 0, x = 0, L, 0 < y < H (3.115) ϑ = 0, y = 0, 0 < x < L (3.116) ϑ = f ( x) − T1 , y = H , 0 < x < L (3.117) It can be seen that there is only one non-homogeneous boundary condition, eq. (3.117). This problem can be solved using the method of separation of variables. In this methodology, it is assumed that the temperature can be expressed as ϑ ( x, y ) = X ( x)Y ( y ) (3.118) where X(x) and Y(y) are functions of x and y, respectively. Substituting eq. (3.118) into eq. (3.110), one obtains 1 d 2 X 1 d 2Y − = (3.119) X dx 2 Y dy 2
Chapter 3 Heat Conduction 233
Since the left hand side of (3.119) is a function of x only and the right-hand side is function of y only, the only way that eq. (3.119) is valid for any x and y is that both sides equal to a constant – termed separation constant. If this separation constant is represented by λ 2 , which is unknown at this point, eq. (3.119) will become two ordinary differential equations. d2X + λ2 X = 0 (3.120) dx 2 d 2Y − λ 2Y = 0 (3.121) dy 2 The general solutions of eqs. (3.120) and (3.121) are X = C1 cos λ x + C2 sin λ x (3.122) Y = C3e − λ y + C4 eλ y (3.123) where C1, C2, C3 and C4 are integral constants to be determined by appropriate boundary conditions. Substituting eq. (3.118) into eqs. (3.115) and (3.116), the following boundary conditions of eqs. (3.120) and (3.121) are obtained (3.124) X (0) = 0 X ( L) = 0 (3.125) Y (0) = 0 (3.126) Substituting eq. (3.122) into eq. (3.124) yields C1 = 0 and eq. (3.122) becomes X = C2 sin λ x (3.127) Applying eq. (3.125) to (3.127), one obtains C2 sin λ L = 0 (3.128) While eq. (3.128) cannot be used to determine C2, the separation constant that satisfies eq. (3.125) can be obtained as nπ λn = (3.129) L where n is an integer. Equation (3.129) indicates that there are many possible values for the separation constants. Substituting eq. (3.123) into eq. (3.126), we get C3 = −C4 and eq. (3.123) becomes Y = C4 (eλ y − e − λ y ) = 2C4 sinh λ y (3.130) Substituting eqs. (3.127) and (3.130) into eq. (3.118) and considering eq. (3.129), a solution for the excess temperature is obtained ⎛ nπ x ⎞ ⎛ nπ y ⎞ (3.131) ϑn = Cn sin ⎜ ⎟ sinh ⎜ ⎟ ⎝L⎠ ⎝L⎠ where Cn = 2C2C4 . For any n, eq. (3.131) satisfies eqs. (3.114) – (3.116) but not eq. (3.117). Since the two-dimensional heat conduction problem under
234 Advanced Heat and Mass Transfer
consideration is a linear problem, the sum of different ϑn for each value of n also satisfies eqs. (3.114) – (3.116). ∞ ⎛ nπ x ⎞ ⎛ nπ y ⎞ ϑ = Cn sin ⎜ (3.132) ⎟ sinh ⎜ ⎟ L⎠ ⎝ ⎝L⎠ n =1 where the term for n = 0 is not included because sinh(0) = 0. Applying the boundary condition at the top of the domain, eq. (3.117), yields ∞ ⎛ nπ x ⎞ ⎛ nπ H ⎞ f ( x) − T1 = Cn sin ⎜ ⎟ sinh ⎜ ⎟ ⎝L⎠ ⎝L⎠ n =1 Multiplying the above equation by sin ( mπ x / L ) and integrating the resulting equation in the interval of (0, L), one obtains L ⎛ mπ x ⎞ [ f ( x) − T1 ]sin ⎜ ⎟ dx 0 ⎝L⎠ (3.133) ∞ ⎛ nπ H ⎞ L ⎛ nπ x ⎞ ⎛ mπ x ⎞ Cn sinh ⎜ = ⎟ sin ⎜ ⎟ sin ⎜ ⎟ dx ⎝ L ⎠0 ⎝L⎠ ⎝L⎠ n =1 The integral in the right-hand side of eq. (3.133) can be evaluated as L ⎡ sin(nπ x / L − mπ x / L) ⎛ nπ x ⎞ ⎛ mπ x ⎞ sin ⎜ ⎟ sin ⎜ ⎟ dx = ⎢ 0 ⎝L⎠ ⎝L⎠ ⎣ 2(nπ x / L − mπ x / L)
∑
∑
∫
∑
∫
∫
sin(nπ x / L + mπ x / L) ⎤ − = 0 (for m ≠ n) 2(nπ x / L + mπ x / L) ⎥ 0 ⎦
L
L ⎡ mπ x 1 ⎛ 2mπ x ⎞ ⎤ L ⎛ mπ x ⎞ sin ⎜ ⎟ dx = ⎢ L − 2 sin ⎜ L ⎟ ⎥ = 2 (for m = n) 0 2mπ ⎣ ⎝L⎠ ⎝ ⎠⎦0 This is referred to as the orthogonal property of the sine function. This property indicates that only one term (m = n) on the right hand of eq. (3.133) is not zero and eq. (3.133) becomes L L ⎛ mπ x ⎞ ⎛ mπ H ⎞ [ f ( x) − T1 ]sin ⎜ ⎟ dx = Cm sinh ⎜ ⎟ 0 2 ⎝L⎠ ⎝L⎠ i.e., 1 2L ⎛ mπ x ⎞ Cm = [ f ( x) − T1 ]sin ⎜ ⎟ dx sinh ( mπ H / L ) L 0 ⎝L⎠ Changing notation from m to n, we get 1 2L ⎛ nπ x ⎞ (3.134) Cn = [ f ( x) − T1 ]sin ⎜ ⎟ dx sinh ( nπ H / L ) L 0 ⎝L⎠ Thus, the temperature profile, eq. (3.132), becomes ∞ sinh ( nπ y / L ) 2 ⎛ nπ x ⎞ L ⎛ mπ x ⎞ ϑ= sin ⎜ ⎟ 0 [ f ( x) − T1 ]sin ⎜ ⎟ dx (3.135) L n =1 sinh ( nπ H / L ) ⎝ L ⎠ ⎝L⎠
∫
L
L
2
∫
∫
∫
∑
∫
Chapter 3 Heat Conduction 235
If the temperature at the top of the domain is constant, i.e,. f ( x) = T2 , the temperature distribution becomes T − T1 2 ∞ (−1) n +1 + 1 sinh ( nπ y / L ) ⎛ nπ x ⎞ (3.136) = sin ⎜ ⎟ T2 − T1 π n =1 2 sinh ( nπ H / L ) ⎝ L ⎠
∑
Two-Dimensional Heat Conduction with Internal Heat Generation
The condition under which the two-dimensional heat conduction can be solved by separation of variables is that the governing equation must be linear homogeneous and no more than one boundary condition is nonhomogeneous. For the case that an internal heat source is present, the energy equation is no longer homogeneous and the method of separation of variables will not work. Let us consider heat conduction in a rectangular domain with a dimension of 2 L × 2 H with uniform heat generation (see Fig. 3.11). All four boundaries are kept at a constant temperature of T1. Since this is an axisymmetric problem, we can study one fourth of the problem. The governing equations and corresponding boundary conditions are ∂ 2T ∂ 2T q′′′ + + = 0, 0 < x < L, 0 < y < H (3.137) ∂x 2 ∂y 2 k with the following boundary conditions: ∂T = 0, x = 0, 0 < y < H (3.138) ∂x T = T1 , x = L, 0 < y < H (3.139) ∂T = 0, y = 0, 0 < x < L (3.140) ∂y T = T1 , y = H , 0 < x < L (3.141) y H T1 –L 0 T1 q′′′ T1 L x
–H
T1
Figure 3.11 Two-dimensional steady-state heat conduction with internal heat generation
236 Advanced Heat and Mass Transfer
Although eq. (3.137) – (3.141) are all linear, eqs. (3.137), (3.139) and, (3.141) are nonhomogeneous. Introducing excess temperature ϑ = T − T1 eqs. (3.137) – (3.141) become ∂ 2ϑ ∂ 2ϑ q′′′ + + = 0, 0 < x < L, 0 < y < H (3.142) k ∂x 2 ∂y 2 ∂ϑ = 0, x = 0, 0 < y < H (3.143) ∂x ϑ = 0, x = L, 0 < y < H (3.144) ∂ϑ = 0, y = 0, 0 < x < L (3.145) ∂y ϑ = 0, y = H , 0 < x < L (3.146) Equation (3.142) is still nonhomogeneous but all boundary conditions are homogeneous. Mathematically, the solution of a nonhomogeneous problem can be obtained by superposing a particular solution of the nonhomogeneous problem and the general solution of the corresponding homogeneous problem. Assuming the solution is ϑ ( x, y ) = ϕ ( x ) + ψ ( x, y ) (3.147) where φ(x) is the solution of the following problem d 2ϕ q′′′ + = 0, 0 < x < L (3.148) dx 2 k dϕ = 0, x = 0 (3.149) dx ϕ = 0, x = L (3.150) and ψ ( x, y ) is the solution of the following problem ∂ 2ψ ∂ 2ψ + = 0, 0 < x < L, 0 < y < H (3.151) ∂x 2 ∂y 2 ∂ψ = 0, x = 0, 0 < y < H (3.152) ∂x ψ = 0, x = L, 0 < y < H (3.153) ∂ψ = 0, y = 0, 0 < x < L (3.154) ∂y ψ = −ϕ ( x), y = H , 0 < x < L (3.155) It can be shown that the summation of the solutions of the above two simpler problems gives the solution of the original problem (see Problem 3.16). Integrating eq. (3.148) twice and determining the integral constants using eqs. (3.149) and (3.150), we obtain q′′′L2 ⎛ x 2 ⎞ ϕ ( x) = (3.156) ⎜1 − ⎟ 2k ⎝ L2 ⎠
Chapter 3 Heat Conduction 237
The solution of eqs. (3.151) – (3.155) can be obtained by using the method of separation of variables and the result is ∞ 2q′′′L2 (−1) n cosh(λn y ) cos(λn x) ψ ( x, y ) = − (3.157) k n = 0 (λn L)3 cosh(λn H ) where the eigenvalue is (2n + 1)π λn = (3.158) , n = 0,1, 2,... 2L Substituting eqs. (3.156) and (3.157) into eq. (3.147), the solution of heat conduction in a rectangular domain with internal heat generation can be obtained. The solution of heat conduction with internal heat generation was obtained by solving a 1-D heat conduction problem in the x-direction with an internal heat source and a 2-D heat conduction problem without internal heat generation. Alternatively, the solution can also be assumed to have the following form ϑ ( x, y ) = ϕ ( y ) + ψ ( x , y ) (3.159) and an equivalent but different solution can be obtained (see Problem 3.17).
∑
Multidimensional Heat Conduction from Buried and Surface Objects
The above two-dimensional heat conduction occurs in a regular shape such that the method of separation of variables can be easily employed. Let us turn our attention to a more complicated problem – an infinitely long cylinder (in the direction perpendicular to the paper) buried in homogeneous medium extending to infinity (see Fig. 3.12). This configuration can find its application in buried large-current cable or steam pipe for residential or industrial heating. The cylinder surface temperature is uniformly T0 and the surface temperature is uniformly Ts. Our objective is to find the temperature distribution, T(x, y) and the total amount of heat transfer from the buried cylinder. Ts y z r0 T0 x
Figure 3.12 Buried long cylinder
238 Advanced Heat and Mass Transfer
Equation (3.110) is still valid for this problem and the boundary conditions for this problem are T = Ts , y = 0 (3.160)
T = T0 , x 2 + ( y − z ) 2 = r02 (3.161) The solution technique discussed previously cannot be applied because the heat transfer domain is very irregular and cannot be easily described using either Cartesian or cylindrical coordinate systems. A different method that uses heat source and sink (see Fig. 3.13) will be used to obtain the solution (Eckert and Drake, 1987). In this methodology, the cylinder is treated as a line heat source, q′ , located at a depth a ( a < z ). If the conduction medium extends infinitely in both the x- and y-direction, the isotherms caused by this line heat source will be a series of circles with its center at y = a but eq. (3.160) cannot be satisfied. In order to fulfill the boundary condition at the surface y = 0, a fictitious heat sink, − q′ , which is symmetric to the heat source, must be used. The combination of the heat source and sink will yield a constant surface temperature, Ts, at y = 0.
r0
z
a
Ts a y x
r2
z
r1 r0
P(x, y)
T0
Figure 3.13 Heat source and sink method
Chapter 3 Heat Conduction 239
The heat generated from the heat source located at (0, a) will keep the surface temperature of the cylinder equal to T0,1. If the temperature at the point P caused by the heat source is T1(x, y), the amount of heat transfer from the surface of the cylinder is 2π k (T0,1 − T1 ) q′ = ln(r1 / r0 ) The temperature at point P caused by the heat source alone is then, ⎛r ⎞ q′ T1 = T0,1 − ln ⎜ 1 ⎟ (3.162) 2π k ⎝ r0 ⎠ Similarly, the heat sink located at ( 0, − a ) will maintain the surface temperature of the fictitious cylinder at T0,2 and result in the temperature at the point P equal to T2(x, y), i.e., 2π k (T0,2 − T2 ) − q′ = ln(r2 / r0 ) Consequently, the temperature at point P caused by the heat sink alone is, ⎛r ⎞ q′ T2 = T0,2 + ln ⎜ 2 ⎟ (3.163) 2π k ⎝ r0 ⎠ If both heat source and sink are present, the temperature at point P is T ( x, y ) = T1 ( x, y ) + T2 ( x, y ) (3.164) Substituting eqs. (3.162) and (3.163) into eq. (3.164) yields ⎛r ⎞ q′ T = (T0,1 + T0,2 ) + ln ⎜ 2 ⎟ (3.165) 2π k ⎝ r1 ⎠ At the surface y = 0 where r1 = r2 , eq. (3.160) requires T0,1 + T0,2 = Ts and eq. (3.165) becomes ⎛r ⎞ q′ ln ⎜ 2 ⎟ T = Ts + 2π k ⎝ r1 ⎠ Defining excess temperature θ = T − Ts , we have
θ=
⎛r ⎞ q′ ln ⎜ 2 ⎟ 2π k ⎝ r1 ⎠
(3.166)
The radii r1 and r2 can be expressed as r1 = x 2 + ( y − z ) 2 x 2 + ( y − a)2 (3.167) r2 = x 2 + ( y + z ) 2 x 2 + ( y + a)2 (3.168) which are valid because a is very close to z as will become evident later. Thus, the excess temperature in the conducting medium can be expressed as ⎡ x 2 + ( y + a)2 ⎤ ⎡ x 2 + ( y + a)2 ⎤ q′ q′ θ= ⎥= ln ⎢ ln ⎢ 2 ⎥ (3.169) 2π k ⎢ x 2 + ( y − a) 2 ⎥ 4π k ⎣ x + ( y − a) 2 ⎦ ⎣ ⎦
240 Advanced Heat and Mass Transfer
which can be rewritten as
⎛ 4π kθ ⎞ x 2 + ( y + a) 2 exp ⎜ (3.170) ⎟= 2 2 ⎝ q′ ⎠ x + ( y − a ) When the left hand side of eq. (3.170) is held at a constant C, the following equation for isotherms is obtained: x 2 + ( y + a)2 =C x 2 + ( y − a)2 which can be rewritten as ⎡ 4Ca 2 ⎛ 1 + C ⎞⎤ (3.171) x + ⎢y + a⎜ = ⎟⎥ (1 − C ) 2 ⎝ 1 − C ⎠⎦ ⎣ It represents a series of circles with the center located at [0, − a(1 + C ) /(1 − C ) ]
2 2
and radius of 2 Ca /(1 − C ) ; both the center of radius of the circular isotherms depend on the constant C, or the temperature. For the cylinder surface, the excess temperature is θ 0 = T0 − Ts and the constant C is ⎛ 4π kθ 0 ⎞ C0 = exp ⎜ (3.172) ⎟ ⎝ q′ ⎠ The depth of the cylinder measured from the horizontal surface to the center of the cylinder is 1 + C0 (3.173) z= a 1 − C0 which can be rearranged to 1 − C0 a= z (3.174) 1 + C0 It is evident from eq. (3.174) that the depth of the line heat source must be shallower than that of the center of the cylinder because (1 − C0 ) /(1 + C0 ) < 1 . The radius of the cylinder is related to the constant C0 by [see eq. (3.171)] 2 C0 a 2 C0 z r0 = = (3.175) (1 − C0 ) (1 + C0 ) Solving C0 in term of r0 yields
⎛z⎞ z (3.176) C0 = + ⎜ ⎟ − 1 r0 ⎝ r0 ⎠ Substituting eq. (3.172) into eq. (3.176) and solving for the heat transfer rate, one obtains 2π k (3.177) q′ = (T0 − Ts ) 2 ln ⎡ z / r0 + ( z / r0 ) − 1 ⎤ ⎢ ⎥ ⎣ ⎦ which can be rewritten as
2
Chapter 3 Heat Conduction 241
q′ =
2π k (T0 − Ts ) cosh −1 ( z / r0 )
(3.178)
For a cylinder with length of L ( L 2r0 ), the total heat transfer rate can be expressed as Q = q′L = kS ΔT (3.179) where 2π L 2π L (3.180) = S= ⎡ z / r + ( z / r )2 − 1 ⎤ cosh −1 ( z / r0 ) ln 0 ⎢0 ⎥ ⎣ ⎦ is referred to as the shape factor – a parameter that is related to the geometric configuration of the problem only. Since the isotherms for the above problem are circles at different centers, heat conduction in a cylinder with an eccentric circular bore (see Fig. 3.14) can be analyzed using eq. (3.177). In this methodology, the inner and outer cylinders can be viewed as two isotherms caused by a line heat source at (0, a) and a heat sink with the same intensity at (0, –a). Under steady-state, the heat transfer rate across both the isotherms T1 and T2 must be identical, i.e., (T1 − Ts ) (T1 − Ts ) (3.181) q′ = = 2 ⎡ z / r + ( z / r ) − 1 ⎤ [1/(2π k )]cosh −1 ( z1 / r1 ) [1/(2π k )]ln 1 1 11 ⎢ ⎥ ⎣ ⎦ (T2 − Ts ) (T2 − Ts ) (3.182) q′ = = ⎡ z / r + ( z / r )2 − 1 ⎤ [1/(2π k )]cosh −1 ( z2 / r2 ) [1/(2π k )]ln 2 2 2 2 ⎢ ⎥ ⎣ ⎦
Ts y z2 z1 a
x
r1 r2 T1 T2
ε
Figure 3.14 Heat conduction in a cylinder with eccentric bore
242 Advanced Heat and Mass Transfer
Equations (3.181) and (3.182) can be combined to eliminate Ts to yield 2π k (T1 − T2 ) 2π k (T1 − T2 ) q′ = (3.183) = −1 ⎡ z / r + ( z / r )2 − 1 ⎤ cosh ( z1 / r1 ) − cosh −1 ( z2 / r2 ) 11 1 1 ⎥ ln ⎢ ⎢ z / r + ( z / r )2 − 1 ⎥ 2 2 ⎣2 2 ⎦ It is desirable to eliminate z1 and z2 from eq. (3.183) so that the heat transfer rate will be related to the size and locations of cylinders only. It follows from Fig. 3.14 that 2 z12 − r12 = z2 − r22 = a 2 z2 − z1 = ε where ε is the distance between the centers of the two cylinders and is referred to as the eccentricity. Solving for z1 and z2 from the above two equations yields r2 − r2 − ε 2 z1 = 2 1 (3.184) 2ε r2 − r2 + ε 2 z2 = 2 1 (3.185) 2ε Substituting eqs. (3.184) and (3.185) into eq. (3.183), one obtains 2π k (T1 − T2 ) q′ = (3.186) −1 cosh ⎡ (8r12 + 8r22 − ε 2 ) /(4r1r2 ) ⎤ ⎣ ⎦ If the length of the cylinder and cylindrical bore is L ( L 2r0 ), the total heat transfer rate can be obtained from eq. (3.179) with ΔT = T1 − T2 and the shape factor of 2π L (3.187) S= −1 2 cosh ⎡ (8r1 + 8r22 − ε 2 ) /(4r1r2 ) ⎤ ⎣ ⎦ Heat transfer rate for many multidimensional conduction problems can be calculated using eq. (3.179) and appropriate shape factors. The shape factors for selected geometric configurations can be found in Incropera et al. (2006). More complete results for shape factors can be found in heat transfer handbooks, e.g., Rohsenow et al. (1998).
3.3 Unsteady Heat Conduction
Our discussions thus far have been limited to the case that the temperature is not a function of time. For many applications, it is necessary to consider the variation of temperature with time. In this case, the energy equation for classical heat conduction, eq. (3.8), should be solved. If the thermal conductivity is independent from the temperature, the energy equation is reduced to eq. (3.10). Analysis of transient heat conduction using lumped (zero-dimensional), onedimensional, and multidimensional models will be presented in this section.
Chapter 3 Heat Conduction 243
3.3.1 Lumped Analysis
Let us consider an arbitrarily shaped object with volume V, surface area As, and a uniform initial temperature of Ti as shown in Fig. 3.15. At time t = 0, the arbitrarily shaped object is exposed to a fluid with temperature of T∞ ( T∞ < Ti ) and the convective heat transfer coefficient between the fluid and the arbitrarily shaped object is h. Since the object is cooled from the surface, it is expected that the temperature at the center will be higher than that at the surface. We will consider a special case in this subsection that the temperature in the arbitrarily shaped object is uniform at any time – referred to as the lumped capacitance method. Qualitatively speaking, this assumption is valid when the thermal conductivity of the arbitrarily shaped object is very large or when its size is very small. The quantitative criterion about the validity of lumped capacitance method will be discussed later in this subsection. The convective heat transfer from the surface is q = hAs (T − T∞ ) Similar to what we did for heat transfer from an extended surface, the heat loss due to surface convection can be treated as an equivalent heat source hA (T − T∞ ) q q′′′ = − = − s V V Since the temperature is assumed to be uniform, the energy equation (3.10) becomes hA (T − T∞ ) ∂T =− s (3.188) ρcp ∂t V which is subject to the following initial condition T = Ti , t = 0 (3.189) ϑ = T − T∞ , eqs. (3.188) and (3.189) become Introducing excess temperature, hAs dϑ =− ϑ (3.190) dt ρVc p
ϑ = ϑi , t = 0
h, T∞ As
(3.191)
T(t) V
Figure 3.15 Lumped capacitance method
244 Advanced Heat and Mass Transfer
Integrating eq. (3.190) and determining the integral constant using eq. (3.191), the solution becomes
ϑ = e−t /τ ϑi
t
(3.192)
where
(3.193) hAs is referred to as the thermal time constant. At time t = τ t , the ratio of excess temperature and initial excess temperature is θ / θ i = e −1 = 0.368 . The value of the thermal time constant is a measure of how fast the temperature of the object reacts to its thermal environment. The cooling process requires transferring heat from the center of the object to the surface and the further transfering heat away from the surface by convection. When the lumped capacitance method is employed, it is assumed that the conduction resistance within the object is negligible compared with the convective thermal resistance at the surface, therefore, the validity of the lumped analysis depends on the relative thermal resistances of conduction and convection. The conduction thermal resistance can be expressed as L Rcond = c kAc where Lc is the characteristic length and Ac is the area of heat conduction. The convection thermal resistance at the surface is 1 Rconv = hAs Assuming As = Ac , one can define the Biot number, Bi, as the ratio of the conduction and convection thermal resistances R hLc Bi= cond (3.194) Rconv k If the characteristic length is chosen as Lc = V / As , the lumped capacitance method is valid when the Biot number is less than 0.1, or the conduction thermal resistance is one order of magnitude smaller than the convection thermal resistance at the surface.
τt =
ρVc p
3.3.2 One Dimensional Transient Systems
For the case that the Biot number is greater than 0.1, the temperature distribution can no longer be treated as uniform and the knowledge about the temperature distribution is also of interest. We will now consider the situation where temperature only varies in one spatial dimension. Both homogeneous and nonhomogeneous problems will be considered.
Chapter 3 Heat Conduction 245
Homogeneous Problems Cartesian Coordinate System Figure 3.16 shows a finite slab with thickness of L and a uniform initial temperature of Ti. At time t = 0, the left side of the slab is insulated while the right side of the slab is exposed to a fluid with temperature of T∞ ( T∞ < Ti ). The convective heat transfer coefficient between the fluid and finite slab is h. In contrast to the lumped capacitance method that assumes uniform temperature, we will present a more generalized model that takes non-uniform temperature distribution in the slab into account. It is assumed that there is no internal heat generation in the slab. The energy equation for this one-dimensional transient conduction problem is ∂ 2T 1 ∂T = , 0 < x < L, t > 0 (3.195) ∂x 2 α ∂t subject to the following boundary and initial conditions ∂T = 0, x = 0 (3.196) ∂x ∂T −k = h(T − T∞ ), x = L (3.197) ∂x T = Ti , 0 < x < L, t = 0 (3.198) It should be pointed out that eqs. (3.195) – (3.198) are also valid for the case that both sides of a finite slab with thickness of 2L are cooled by convection. This is a nonhomogeneous problem because eq. (3.197) is not homogeneous. By introducing the excess temperature, ϑ = T − T∞ , the problem can be homogenized, i.e., ∂ 2ϑ 1 ∂ϑ = , 0 < x < L, t > 0 (3.199) ∂x 2 α ∂t ∂ϑ = 0, x = 0 (3.200) ∂x
adiabatic
h, T∞
0
Figure 3.16 Transient conduction in a finite slab
L
x
246 Advanced Heat and Mass Transfer
∂ϑ = hϑ , x = L (3.201) ∂x ϑ = ϑi = Ti − T∞ , 0 < x < L, t = 0 (3.202) To express our solution in a compact form so that it can be used for all similar problems, one can define the following dimensionless variables ϑ x αt θ = , X = , Fo = 2 (3.203) ϑi L L and eqs. (3.199) – (3.202) will be nondimensionalized as ∂ 2θ ∂θ = (3.204) , 0 < X < 1, Fo > 0 2 ∂X ∂Fo ∂θ = 0, X = 0 (3.205) ∂X ∂θ − = Biθ , X = 1 (3.206) ∂X θ = 1, 0 < X < 1, Fo = 0 (3.207) This problem can be solved using the method of separation of variables. Assuming that the temperature can be expressed as θ ( X , Fo) = Θ( X )Γ(Fo) (3.208) where Θ and Γ are functions of X and Fo, respectively, eq. (3.204) becomes Θ′′( X ) Γ′(Fo) = Θ( X ) Γ(Fo) Since the left hand side is a function of X only and the right-hand side of the above equation is a function of Fo only, both sides must be equal to a separation constant, μ , i.e., Θ′′( X ) Γ′(Fo) = =μ (3.209) Θ( X ) Γ(Fo) The separation constant μ can be either a real or a complex number. The solution of Γ from eq. (3.209) will be Γ = e μ Fo . If μ is a positive real number, we will have Γ → ∞ when Fo → ∞ , which does not make sense, therefore, μ cannot be a positive real number. If μ is zero, we will have Γ = const , and Θ is a linear function of X only. The final solution for θ will also be a linear function of X only, which does not make sense either. It can also be shown that the separation constant cannot be a complex number (see Problem 3.24), therefore, the separation variable has to be a negative real number. If we represent this negative number by μ = −λ 2 , eq. (3.209) can be rewritten as the following two equations Θ′′ + λ 2 Θ = 0 (3.210) 2 Γ′ + λ Γ = 0 (3.211) The general solutions of eqs. (3.210) and (3.211) are Θ = C1 cos λ X + C2 sin λ X (3.212) −k
Chapter 3 Heat Conduction 247
Γ = C3e − λ Fo (3.213) where C1, C2, and C3 are integral constants. Substituting eq. (3.208) into eqs. (3.205) and (3.206), the following boundary conditions of eq. (3.210) are obtained Θ′(0) = 0 (3.214) −Θ′(1) = BiΘ(1) (3.215) Substituting eq. (3.212) into eq. (3.214) yields Θ′(0) = −C1λ sin(0) + C2 λ cos(0) = C2 λ = 0 Since λ cannot be zero, C2 must be zero and eq. (3.212) becomes Θ = C1 cos λ X (3.216) Applying the convection boundary condition, eq. (3.215), one obtains = cot λn (3.217) Bi where n is an integer. Equation (3.217) indicates that many possible values for λ – termed eigenvalues – can satisfy the convection condition. The eigenvalue λn can be obtained by solving eq. (3.217) using an iterative procedure. The dimensionless temperature with eigenvalue λn can be obtained by substituting eqs. (3.216) and (3.213) into eq. (3.208), i.e., 2 θ n = Cn cos ( λn X ) e− λn Fo (3.218) where Cn = C1C3 . Equation (3.218) is a solution that satisfies eqs. (3.204) – (3.206). At this point, the constant Cn is still unspecified and eq. (3.207) is unused, however, if we substitute eq. (3.218) into (3.207), the constant Cn that satisfies the initial condition cannot be found. Since the one-dimensional transient heat conduction problem under consideration is a linear problem, the sum of different θ n for each value of n also satisfies eqs. (3.204) – (3.206).
2
λn
θ=
∑ C cos ( λ X ) e
n n n =1
∞
2 − λn Fo
(3.219)
Substituting eq. (3.219) into eq. (3.207) yields 1=
∑ C cos ( λ X )
n n n =1
∞
Multiplying the above equation by cos ( λm X ) and integrating the resulting equation in the interval of (0, 1), one obtains
∫
1
0
cos(λm X )dX =
∑ ∫ cos(λ
Cn
1 n =1 0
∞
m
X ) cos(λn X )dX
(3.220)
The integral in the right-hand side of eq. (3.220) can be evaluated as
248 Advanced Heat and Mass Transfer
⎧ λm sin λm cos λn − λn cos λm sin λn , m≠n 2 ⎪ 1 λm − λn2 ⎪ (3.221) cos(λm X )cos(λn X ) dX = ⎨ 0 ⎪ 1 ⎛ λ + sin 2λm ⎞ , m=n ⎪ 2λm ⎜ m 2⎟ ⎝ ⎠ ⎩ Equation (3.217) can be rewritten as Bi = λn tan λn Similarly, for eigenvalue λm , we have Bi = λm tan λm Combining the above two equations, we have λm tan λm − λn tan λn = 0 or λm sin λm cos λn − λn cos λm sin λn = 0 therefore, the integral in eq. (3.221) is zero for the case that m ≠ n , and the right hand side of eq. (3.220) becomes ∞ 1 C⎛ sin 2λm ⎞ Cn cos(λm X ) cos(λn X )dX = m ⎜ λm + (3.222) 0 2λm ⎝ 2⎟ ⎠ n =1 Substituting eq. (3.222) into eq. (3.220) and evaluating the integral at the lefthand side of eq. (3.220), we have C⎛ sin 2λm ⎞ 1 sin λm = m ⎜ λm + 2λm ⎝ 2⎟ λm ⎠ i.e., 4 sin λm Cm = 2λm + sin 2λm Changing notation from m to n, we get 4 sin λn Cn = (3.223) 2λn + sin 2λn The dimensionless temperature, therefore, becomes ∞ 2 4sin λn (3.224) θ= cos ( λn X ) e − λn Fo n =1 2λn + sin 2λn If the Biot number becomes infinite, the convection boundary condition becomes θ = 0, X = 1 (3.225) which is an isothermal condition at the right-hand side of the wall. Equation (3.217) becomes cos λn = 0 (3.226) and the eigenvalue is therefore λn = ( nπ − π / 2 ) , n = 1, 2,3,... (3.227) The temperature distribution for this case is then
∫
∑∫
∑
Chapter 3 Heat Conduction 249
θ=
∑ 2(n − π / 2) + sin 2(n − π / 2) cos [(n − π / 2) X ] e
n =1
∞
4sin(n − π / 2)
− ( n −π / 2) 2 Fo
(3.228)
When Fourier’s number is greater than 0.2, only the first term in eq. (3.219) is necessary and the solution becomes
θ = θ1 = C1 cos ( λ1 X ) e− λ Fo
2 1
(3.229)
which is source of Hiesler’s charts (Incropera et al., 2007).
Cylindrical and Spherical Coordinate Systems The method of separation of variables can also be applied to solve heat conduction in cylindrical or spherical coordinate systems. Its application on transient conduction in a cylindrical system will be illustrated in the following example. Example 3.3 A long cylinder with radius of ro and a uniform initial temperature of Ti is exposed to a fluid with temperature of T∞ ( T∞ < Ti ). The convective heat transfer coefficient between the fluid and cylinder is h. Assuming that there is no internal heat generation and constant thermophysical properties, obtain the transient temperature distribution in the cylinder. Solution: Since the temperature changes along the r-direction only, the energy equation is ∂ 2T 1 ∂T 1 ∂T + = , 0 < x < ro , t > 0 (3.230) ∂r 2 r ∂r α ∂t subject to the following boundary and initial conditions ∂T = 0, r = 0 (axisymmetric) (3.231) ∂r ∂T −k = h(T − T∞ ), r = ro (3.232) ∂r T = Ti , 0 < r < ro , t = 0 (3.233) Defining the following dimensionless variables hr T − T∞ αt r θ= , R = , Fo = 2 , Bi = o (3.234) Ti − T∞ ro ro k eqs. (3.230) – (3.233) will be nondimensionalized as ∂ 2θ 1 ∂θ ∂θ + = , 0 < R < 1, Fo > 0 (3.235) ∂R 2 R ∂R ∂Fo ∂θ = 0, R = 0 (3.236) ∂R ∂θ − = Biθ , R = 1 (3.237) ∂R
250 Advanced Heat and Mass Transfer
θ = 1, 0 < R < 1, Fo = 0 (3.238) Assuming that the temperature can be expressed as θ ( R, Fo) = Θ( R )Γ(Fo) (3.239) and substituting eq. (3.239) into eq. (3.235), one obtains 1⎡ 1 ⎤ Γ′ 2 (3.240) ⎢Θ′′ + R Θ′⎥ = Γ = −λ Θ⎣ ⎦ which can be rewritten as the following two equations 1 Θ′′ + Θ′ + λ 2 Θ = 0 (3.241) R Γ′ + λ 2 Γ = 0 (3.242) Equation (3.101) is a Bessel’s equation of zero order and has the following general solution Θ( R ) = C1 J 0 (λ R ) + C2Y0 (λ R ) (3.243) where J0 and Y0 are Bessel functions of the first and second kind, respectively. The general solution of eq. (3.242) is
Γ = C3e − λ Fo (3.244) where C1, C2, and C3 are integral constants. The boundary conditions for eq. (3.241) can be obtained by substituting eq. (3.239) into eqs. (3.236) and (3.237), i.e., Θ′(0) = 0 (3.245) −Θ′(1) = BiΘ(1) (3.246) The derivative of Θ is Θ′( R) = −C1λ J1 (λ R) − C2 λY1 (λ R) (3.247) Since J1 (0) = 0 and Y1 (0) = −∞ , C2 must be zero. Substituting eqs. (3.243) and (3.247) into eq. (3.246) and considering C2 = 0, we have −λn J1 (λn ) + BiJ 0 (λn ) = 0 (3.248) where n is an integer. The eigenvalue λn can be obtained by solving eq. (3.248) using an iterative procedure. The dimensionless temperature with eigenvalue λn is
2
θ n = Cn J 0 ( λn R ) e − λ Fo
2 n
(3.249)
where Cn = C1C3 . Equation (3.249) is a solution that satisfies eqs. (3.235) – (3.237) but not eq. (3.238). For a linear problem, the sum of different θ n for each value of n also satisfies eqs. (3.235) – (3.237).
θ=
∑C J
n n =1
∞
∞
0
( λn R ) e− λ Fo
2 n
(3.250)
Substituting eq. (3.250) into eq. (3.238) yields 1=
∑C J
n n =1
0
( λn R )
Chapter 3 Heat Conduction 251
Multiplying the above equation by RJ 0 ( λm R ) and integrating the resulting equation in the interval of (0, 1), one obtains
∫
1
0
RJ 0 (λm R )dR =
∑ ∫ RJ (λ R) J (λ R)dR
Cn
1 n =1 0 0 m 0 n
∞
According to the orthogonal property of Bessel’s function, the integral on the right-hand side equals zero if m ≠ n but it is not zero if m = n. Therefore, we have
Cm =
∫ RJ (λ R)dR = 2 λ ∫ RJ (λ R)dR
0 1 0 m 0 2 0 m
1
m
J1 (λm ) J (λm ) + J12 (λm )
2 0
Changing notation from m to n, we get J1 (λn ) 2 Cn = 2 λn J 0 (λn ) + J12 (λn ) thus, the dimensionless temperature becomes ∞ 2 J1 (λn ) J 0 ( λn R ) − λn2 Fo θ= e 2 2 n =1 λn J 0 (λn ) + J1 (λn )
(3.251)
∑
(3.252)
One-dimensional heat conduction in a spherical coordinate system can be solved by introducing a new dependent variable. Let us consider a sphere with radius of ro and a uniform initial temperature of Ti. It is exposed to a fluid with a temperature of T∞ ( T∞ < Ti ) and the convective heat transfer coefficient between the fluid and finite slab is h. Assuming that there is no internal heat generation and constant thermophysical properties, the governing equation is 1 ∂ 2 (rT ) 1 ∂T = (3.253) , 0 < x < ro , t > 0 α ∂t r ∂r 2 subject to the following boundary and initial conditions ∂T = 0, r = 0 (axisymmetric) (3.254) ∂r ∂T −k = h(T − T∞ ), r = ro (3.255) ∂r T = Ti , 0 < r < ro , t = 0 (3.256) By using the same dimensionless variables defined in eq. (3.234), eqs. (3.253) – (3.256) can be nondimensionalized as 1 ∂ 2 ( Rθ ) ∂θ = (3.257) , 0 < R < 1, Fo > 0 ∂Fo R ∂R 2 ∂θ = 0, R = 0 (3.258) ∂R ∂θ − = Biθ , R = 1 (3.259) ∂R
252 Advanced Heat and Mass Transfer
θ = 1, 0 < R < 1, Fo = 0 (3.260) Defining a new dependent variable (3.261) U = Rθ eqs. (3.257) – (3.260) become ∂ 2U ∂U = , 0 < R < 1, Fo > 0 (3.262) ∂R 2 ∂Fo ∂θ = 0, R = 0 (3.263) ∂R ∂U − = (Bi-1)θ , R = 1 (3.264) ∂R U = R, 0 < R < 1, Fo = 0 (3.265) It can be seen that eqs. (3.262) – (3.265) are identical to the case of heat conduction in a finite slab with a non-uniform initial temperature. This problem can be readily solved by using the method of separation of variables (see Problem 3.28). After the solution is obtained, one can change the dependent variable back to θ and the result is ∞ 2 4[sin(λn ) − λn cos(λn )] 1 (3.266) θ= sin(λn R )e − λn Fo R n =1 λn [2λn − sin(2λn )] where the eigenvalue is the positive root of the following equation 1 − λn cot λn = Bi (3.267)
∑
Nonhomogeneous Problems
The transient one-dimensional conduction problems that we discussed so far are limited to the case that the problem is homogeneous and the method of separation of variables works. When the problem is not homogeneous due to a nonhomogeneous energy equation or boundary condition, the solution of a nonhomogeneous problem can be obtained by superposition of a particular solution of the nonhomogeneous problem and the general solution of the corresponding homogeneous problem.
Partial Solution Let us consider a finite slab with thickness of L and a uniform initial temperature of Ti as shown in Fig. 3.17. At time t = 0, the temperature on the left side of the slab is suddenly increased to T0 while the temperature on the right side of the slab is maintained at Ti. Assuming that there is no internal heat generation in the slab and the thermophysical properties of the slab are constants, the energy equation is ∂ 2T 1 ∂T = , 0 < x < L, t > 0 (3.268) ∂x 2 α ∂t
Chapter 3 Heat Conduction 253
T=T0
t→∞ t>0 t=0 0 T=Ti L x
Figure 3.17 Heat conduction under boundary condition of the first kind
subject to the following boundary and initial conditions T = T0 , x = 0 (3.269) T = Ti , x = L (3.270) T = Ti , 0 < x < L, t = 0 (3.271) By defining the following dimensionless variables T − Ti x αt θ= , X = , Fo = 2 (3.272) T0 − Ti L L eqs. (3.268) – (3.271) will be nondimensionalized as ∂ 2θ ∂θ , 0 < X < 1, Fo > 0 = (3.273) 2 ∂X ∂Fo (3.274) θ = 1, X = 0 θ = 0, X = 1 (3.275) θ = 0, 0 < X < 1, Fo = 0 (3.276) It can be seen that the problem is still nonhomogeneous after nondimensionalization because eq. (3.274) is not homogeneous. When t → ∞ (i.e., Fo → ∞ ), the temperature distribution can reach to steady state. If the steady state temperature is represented by θ s , it must satisfy the following equations: ∂ 2θ s = 0, 0 < X < 1 (3.277) ∂X 2 θ s = 1, X = 0 (3.278) θ s = 0, X = 1 (3.279) which have the following solution: θs = 1 − X (3.280) To obtain the generation of the problem described by eqs. (3.273) – (3.276), a method of partial solution (Myers, 1987) will be employed. In this methodology, it is assumed that the solution of a nonhomogeneous problem can be expressed as θ ( X , Fo) = θ s ( X ) + θ h ( X , Fo) (3.281)
254 Advanced Heat and Mass Transfer
where θ h represent the solution of a homogeneous problem. Substituting eqs. (3.273) – (3.276) and considering eqs. (3.277) – (3.279), we have ∂ 2θ h ∂θ h = (3.282) , 0 < X < 1, Fo > 0 ∂X 2 ∂Fo θ h = 0, X = 0 (3.283) θ h = 0, X = 1 (3.284) θ h = X − 1, 0 < X < 1, Fo = 0 (3.285) which represent a new homogeneous problem. This problem can be solved using the method of separation of variables (see Problem 3.31) and the result is 2 ∞ sin(nπ X ) − ( nπ )2 Fo θh = − ∑ e (3.286) π n =1 n The solution of the nonhomogeneous problem thus becomes 2 ∞ sin(nπ X ) − ( nπ )2 Fo θ =1− X − ∑ e (3.287) π n =1 n Variation of Parameter The partial solution only works if the steady-state solution exists. If the steady-state solution does not exist, we can use the method of variation of parameters to solve the problem. Let us consider a finite slab with thickness of L and a uniform initial temperature of Ti. At time t = 0, the left side is subject to a constant heat flux while the right side of the slab is adiabatic (see Fig. 3.18). Assuming that there is no internal heat generation in the slab and the thermophysical properties of the slab are constants, the energy equation is ∂ 2T 1 ∂T = , 0 < x < L, t > 0 (3.288) ∂x 2 α ∂t subject to the following boundary and initial conditions ∂T ′′ −k = q0 , x = 0 (3.289) ∂x
′′ q0
0
L
adiabatic
x
Figure 3.18 Heat conduction under boundary condition of the second kind
Chapter 3 Heat Conduction 255
∂T = 0, x = L (3.290) ∂x T = Ti , 0 < x < L, t = 0 (3.291) By defining the following dimensionless variables T − Ti x αt θ= , X = , Fo = 2 (3.292) ′′ q0 L / k L L eqs. (3.288) – (3.291) will be nondimensionalized as ∂ 2θ ∂θ = (3.293) , 0 < X < 1, Fo > 0 2 ∂X ∂Fo ∂θ = −1, X = 0 (3.294) ∂X ∂θ = 0, X = 1 (3.295) ∂X (3.296) θ = 0, 0 < X < 1, Fo = 0 This nonhomogeneous problem does not have a steady-state solution, and therefore the partial solution cannot be applied. We will use the method of variation of parameters (Myers, 1987) to solve this problem. This method requires the following steps: 1. Set up a homogeneous problem by dropping the nonhomogeneous terms. 2. Solve the homogeneous problem to get eigenvalue λn and eigenfunctions Θ n ( X ) . 3. Assuming the solution of the original nonhomogeneous problem has the form of θ ( X , Fo) =
∑ A (Fo)Θ ( X ) .
n n n =1
n
4. Solve for An(Fo) using orthogonal property of Θ n . 5. Obtain an ordinary differential equation (ODE) for An(Fo) and solve for An(Fo) from the ODE 6. Put together the final solution. We will solve this nonhomogeneous problem by following the above procedure. The corresponding homogeneous problem is: ∂ 2θ h ∂θ h = (3.297) , 0 < X < 1, Fo > 0 ∂X 2 ∂Fo ∂θ h = 0, X = 0 (3.298) ∂X ∂θ h = 0, X = 1 (3.299) ∂X θ h = 0, 0 < X < 1, Fo = 0 (3.300) Assuming the solution of the above homogeneous problem is θ h = Θ( X )Γ(Fo) (3.301)
256 Advanced Heat and Mass Transfer
eq. (3.297) becomes
Θ′′( X ) Γ′(Fo) = = −λ 2 (3.302) Θ( X ) Γ(Fo) Since the objective here is to get the eigenvalue and eigen functions, we do not need to solve for Γ and only need to solve for Θ . The eigenvalue problem is Θ′′ + λ 2 Θ = 0 (3.303) Θ′(0) = 0 (3.304) Θ′(1) = 0 (3.305) Solving eqs. (3.303) – (3.305) yields the following eigenvalues and eigen functions λn = nπ (3.306) Θn ( X ) = cos(nπ X ), n = 0,1, 2,... (3.307) Now, let us assume that the solution of the original nonhomogeneous problem is
θ ( X , Fo) =
∑ A (Fo) cos ( nπ X )
n n=0 ∞ 1
∞
(3.308)
Multiplying eq. (3.308) by cos ( mπ X ) and integrating the resulting equation in the interval of (0, 1), one obtains
∫
1
0
θ ( X , Fo) cos(mπ X )dX =
∑ A ∫ cos(mπ X ) cos(nπ X )dX
n n =1 0
(3.309)
The integral on the right-hand side of eq. (3.309) can be evaluated as m≠n ⎧0, 1 ⎪ cos( mπ X ) cos(nπ X )dX = ⎨1/ 2, m = n ≠ 0 (3.310) 0 ⎪1, m=n=0 ⎩
∫
thus, eq. (3.309) becomes A0 (Fo) =
∫ θ ( X , Fo)dX ,
0 1
1
m=0
(3.311) (3.312)
Am (Fo) = 2 θ ( X , Fo) cos(mπ X )dX , m ≠ 0
0
∫
Differentiating eq. (3.311) with respect to Fo, one obtains: 1 ∂θ dA0 dX = (3.313) 0 ∂Fo dFo Substituting eq. (3.293) into eq. (3.313) and integrating with respect to X yields 1 ∂ 2θ dA0 ∂θ ∂θ dX = = − =1 (3.314) 2 0 ∂X ∂X X =1 ∂X X = 0 dFo Integrating eq. (3.314) with respect to Fo, we have
∫
∫
A0 (Fo) = Fo + C1 = When Fo = 0, eq. (3.315) becomes
∫ θ ( X , Fo)dX
0
1
(3.315)
Chapter 3 Heat Conduction 257
A0 (0) = C1 = thus, we have
∫ θ ( X ,0)dX = 0
0
1
(3.316)
A0 (Fo) = Fo (3.317) Differentiating eq. (3.312) and considering eq. (3.293) yield 1 ∂θ 1 ∂ 2θ dAm =2 cos(mπ X )dX = 2 cos(mπ X )dX (3.318) 2 0 ∂Fo 0 ∂X dFo Using integration by parts twice, the following ODE is obtained: dAm = 2 − (mπ ) 2 Am (3.319) dFo 2 Multiplying eq. (3.319) by an integrating factor e( mπ ) Fo , we have 2 2 d⎡ A e( mπ ) Fo ⎤ = 2e( mπ ) Fo (3.320) ⎣m ⎦ dFo which can be integrated to get 2 2 Am = + C2 e − ( mπ ) Fo (3.321) 2 ( mπ ) where C2 is an integral constant that needs to be determined by an initial condition. For Fo = 0, eq. (3.312) becomes
∫
∫
Am (0) = 2 θ ( X ,0) cos( mπ X ) dX = 2 cos( mπ X ) dX = 0 (3.322)
0 0
∫
1
∫
1
Substituting eq. (3.321) into eq. (3.322), one obtains 2 C2 = − (mπ ) 2 therefore, we have 2 2 2 Am = − e − ( mπ ) Fo 2 2 (mπ ) (mπ ) Changing m back to n for notation, 2 2 2 An = − e − ( nπ ) Fo (3.323) 2 2 (nπ ) (nπ ) Substituting eqs. (3.317) and (3.323) into eq. (3.308), the solution becomes ∞ cos ( nπ X ) 2 ∞ cos ( nπ X ) − ( nπ )2 Fo 2 −2 (3.324) θ ( X , Fo) = Fo+ 2 e π n =1 n2 π n =1 n2 When the time (Fourier number) becomes large, the last term on the right-hand side will become zero and the solution is represented by the first two terms only. To simplify eq. (3.324), let us assume the solution at large Fo can be expressed as θ ( X , Fo) = Fo+f ( X ) (3.325) which is referred to as an asymptotic solution and it must satisfy eqs. (3.293) – (3.295). Substituting eq. (3.325) into eqs. (3.293) – (3.295), we have f ′′( X ) = 1 (3.326)
∑
∑
258 Advanced Heat and Mass Transfer
f ′(0) = −1, f ′(1) = 0 (3.327) Integrating eq. (3.326) and considering eq. (3.327), we obtain X2 f (X ) = − X +C (3.328) 2 where C cannot be determined from eq. (3.327) because both boundary conditions are for the first order derivative. To determine C, we can expand f(X) defined in eq. (3.328) into cosine Fourier series, i.e. ∞ X2 f (X ) = − X + C = a0 + an cos(nπ X ) (3.329) 2 n =1 After determining a0 and an, and considering that f(X) is identical to the second term on the right-hand side of eq. (3.324), we have ∞ cos ( nπ X ) X2 12 −X+ = 2 (3.330) 2 3 π n =1 n2 Substituting eq. (3.330) into eq. (3.324), the final solution becomes ∞ cos ( nπ X ) − ( nπ )2 Fo X2 12 −X+ − 2 (3.331) θ ( X , Fo) = Fo+ e 2 3 π n =1 n2
∑
∑
∑
Transient Heat Conduction in a Semi-Infinite Body
Let us consider heat conduction in a semi-infinite body (x > 0) with an initial temperature of 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 because 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 . The physical model of the problem is illustrated in Fig. 3.19, 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 (3.332) ∂x 2 α ∂t T=T0
t=0 0
T=Ti x
Figure 3.19 Heat conduction in a semi-infinite body
Chapter 3 Heat Conduction 259
T ( x, t ) = T0 x = 0, t > 0 (3.333) T ( x, t ) = Ti x → ∞, t > 0 (3.334) T ( x, t ) = Ti x > 0, t = 0 (3.335) which can be solved by using the method of separation of variables or integral approximate solution.
Separation of Variables Defining the following dimensionless variables T − T0 x αt θ= , X = , Fo = 2 (3.336) Ti − T0 L L where L is a characteristic length, eqs. (3.332) – (3.335) will be nondimensionalized as ∂ 2θ ∂θ = , X > 0, Fo > 0 (3.337) 2 ∂X ∂Fo (3.338) θ = 0, X = 0 θ = 1, X → ∞ (3.339) θ = 1, X > 0, Fo = 0 (3.340) Assuming that the temperature can be expressed as θ ( X , Fo) = Θ( X )Γ(Fo) (3.341) and substituting eq. (3.341) into eq. (3.337), one obtains Θ′′( X ) Γ′(Fo) = = −λ 2 (3.342) Θ( X ) Γ(Fo) whose general solutions are Θ = C1 cos λ X + C2 sin λ X (3.343)
Γ = C3e− λ Fo (3.344) where C1, C2, and C3 are integral constants. Substituting eq. (3.341) into eq. (3.338), the following boundary condition of eq. (3.343) is obtained Θ(0) = 0 (3.345) Substituting eq. (3.343) into eq. (3.345) yields C1 = 0 and eq. (3.343) becomes Θ = C2 sin λ X (3.346) The boundary condition at X → ∞ , eq. (3.339) cannot be used to get the eigen equation for λ and therefore, there will be no restriction on the value of λ for heat conduction in a semi-infinite body. Substituting eqs. (3.346) and (3.344) into eq. (3.341), the solution becomes 2 θ λ = C (λ )sin ( λ X ) e − λ Fo (3.347)
where C = C2 C3 . The general solution for the problem can be obtained by using linear combination of eq. (3.347) for all possible λ , i.e.,
2
260 Advanced Heat and Mass Transfer
θ=
∫
∞
0
C (λ )sin ( λ X ) e − λ Fo d λ
2
(3.348)
Substituting eq. (3.348) into eq. (3.340), one obtains 1=
∫
∞
λ =0
C (λ )sin ( λ X ) d λ
If we solve the problem by using Laplace transformation, we have ∞ ⎡2 ∞ ⎤ 1= sin ( λ X ) ⎢ sin ( λ X ′ ) dX ′⎥ d λ λ =0 ⎣ π X ′=0 ⎦ Comparing the above two equations, an expression of C is obtained: 2∞ C (λ ) = sin ( λ X ′ ) dX ′
∫
∫
π
∫
X ′= 0
The temperature distribution, eq. (3.348), becomes ∞ 2 2∞ θ= e − λ Fo sin ( λ X ′ ) dX ′ sin ( λ X ) d λdX ′
π
∫∫
X ′= 0
λ =0
which can be rewritten as ∞ 2 2∞ θ= e − λ Fo [ cos λ ( X − X ′) − cos λ ( X + X ′)]d λdX ′ Evaluating the integral with respect to λ yields
π
∫∫
X ′= 0
λ =0
∫
∫
thus
⎡ ( X − X ′) 2 ⎤ exp ⎢ − ⎥ λ =0 4Fo 4Fo ⎦ ⎣ ∞ 2 ⎡ ( X + X ′) 2 ⎤ π e − λ Fo cos λ ( X + X ′)d λ = exp ⎢ − ⎥ λ =0 4Fo 4Fo ⎦ ⎣
∞
e−λ
2
Fo
cos λ ( X − X ′)d λ =
π
∞ ⎧∞ ⎫ ⎡ ( X − X ′) 2 ⎤ ⎡ ( X + X ′) 2 ⎤ ⎪ ⎪ ⎨ ′ exp ⎢ − ⎥ dX ′ − X ′= 0 exp ⎢ − ⎥ dX ′⎬ (3.349) X =0 4Fo ⎦ 4Fo ⎦ 4π Fo ⎪ ⎪ ⎣ ⎣ ⎩ ⎭ Let us define a new variable X − X′ η= 4Fo The first integral in eq. (3.349), becomes ∞ ∞ ⎡ ( X − X ′) 2 ⎤ 2 exp ⎢ − ⎥ dX ′ = 4Fo − X / 4Fo exp(−η )dη X ′= 0 4Fo ⎦ ⎣
θ=
1
∫
∫
∫
∫
Similarly, the second integral can be evaluated by following a similar procedure: ∞ ∞ ⎡ ( X + X ′) 2 ⎤ 2 exp ⎢ − ⎥ dX ′ = 4Fo X / 4Fo exp(−η )dη X ′= 0 4Fo ⎦ ⎣
∫
∫
Substituting the above two equations into eq. (3.349), we have ∞ 1⎡∞ 2 X/ ⎤ θ= exp(−η 2 )dη − exp(−η 2 )dη ⎥ = X / 4Fo ⎣ ⎦ π ⎢ − X / 4Fo π0 which can be rewritten as
∫
∫
∫
4Fo
exp(−η 2 )dη
Chapter 3 Heat Conduction 261
X⎞ ⎟ ⎝ 4Fo ⎠ where erf in eq. (3.350) is the error function defined as: 2 z − z2 erf ( z ) = e dz
θ = erf ⎜
⎛
(3.350)
π
∫
0
(3.351)
Equation (3.350) can also be rewritten as dimensional form: T − T0 ⎛x⎞ = erf ⎜ (3.352) ⎟ Ti − T0 ⎝ 4α t ⎠ The surface heat flux can be obtained by applying Fourier’s law k (T0 − Ti ) ∂T ′′ q0 (t ) = −k = (3.353) ∂x x = 0 πα t The solution of heat conduction in a semi-infinite body under the boundary conditions of the second and third kinds can also be obtained by using the method of separation of variables (Ozisik, 1993).
Time-Dependent Boundary Condition Our discussions thus far have been limited to the case that the boundary condition is not a function of time. Periodic boundary conditions can be encountered in various applications ranging from heat conduction in a building during day and night to emerging technologies such as pulsed laser processing of materials. Let us reconsider the problem described by eqs. (3.332) – (3.335) but replace eq. (3.333) by T = Ti + f (t ) = Ti + A cos(ωt − β ), x = 0, t > 0 (3.354) where A is the amplitude of oscillation, ω is the angular frequency, and β is the phase delay. The method of separation of variables will not work since eq. (3.354) is not homogeneous. Introducing excess temperature ϑ = T − Ti , the governing equation and corresponding boundary and initial conditions become ∂ 2ϑ ( x, t ) 1 ∂ϑ ( x, t ) = x > 0, t > 0 (3.355) ∂x 2 α ∂t ϑ ( x, t ) = f (t ) = A cos(ωt − β ) x = 0, t > 0 (3.356) ϑ ( x, t ) = 0 x → ∞, t > 0 (3.357) ϑ ( x, t ) = 0 x > 0, t = 0 (3.358) Duhamel’s theorem can be used to handle the time-dependent boundary condition. Instead of solving eqs. (3.355) – (3.358) directly, we will start with a simpler auxiliary problem defined below: ∂ 2 Φ ( x, t ) 1 ∂Φ ( x, t ) = x > 0, t > 0 (3.359) ∂x 2 α ∂t Φ ( x, t ) = f (τ ) = A cos(ωτ − β ) x = 0, t > 0 (3.360) Φ ( x, t ) = 0 x → ∞, t > 0 (3.361) Φ ( x, t ) = 0 x > 0, t = 0 (3.362)
262 Advanced Heat and Mass Transfer
where τ in eq. (3.360) is treated as a parameter, rather than time. Duhamel’s theorem (Ozisik, 1993) stated that the solution to the original problem is related to the solution of the auxiliary problem by ∂t Φ ( x, t − τ ,τ ) dτ ϑ ( x, t ) = (3.363) ∂t τ = 0 which can rewritten using Leibniz’s rule t ∂ Φ ( x, t − τ ,τ )dτ + Φ ( x, t − τ ,τ ) τ =t ϑ ( x, t ) = (3.364) τ = 0 ∂t The second term on the right hand side is (3.365) Φ ( x, t − τ ,τ ) τ =t = Φ ( x,0,τ ) = 0
∫
∫
therefore, eq. (3.364) becomes ∂ Φ ( x, t − τ ,τ ) dτ (3.366) ∂t The solution of the auxiliary problem can be expressed as 2 ⎡ ⎛ x ⎞ ⎤ 2 f (τ ) ∞ e−η dη (3.367) Φ ( x, t ,τ ) = f (τ ) ⎢1- erf ⎜ ⎟⎥ = π x / 4α t ⎝ 4α t ⎠ ⎦ ⎣ The partial derivative appearing in eq. (3.366) can be evaluated as ⎡ ⎤ x x2 ∂ exp ⎢ − Φ ( x, t − τ ,τ ) = f (τ ) ⎥ (3.368) 3/ 2 ∂t 4πα (t − τ ) ⎣ 4α (t − τ ) ⎦
ϑ ( x, t ) =
∫
t
τ =0
∫
Substituting eq. (3.368) into eq. (3.366), the solution of the original problem becomes t ⎡ ⎤ x f (τ ) x2 exp ⎢ − ϑ ( x, t ) = (3.369) ⎥ dτ 3/ 2 4πα τ = 0 (t − τ ) ⎣ 4α (t − τ ) ⎦ Introducing a new independent variable x ξ= 4α (t − τ ) eq. (3.369) becomes ⎛ 2∞ x2 ⎞ 2 ϑ ( x, t ) = f ⎜t − (3.370) ⎟ exp(−ξ )d ξ x / 4α t 4αξ 2 ⎠ π ⎝
∫
∫
For the periodic boundary condition specified in eq. (3.356), we have ⎡⎛ ⎤ 2A ∞ x2 ⎞ cos ⎢ω ⎜ t − ϑ ( x, t ) = − β ⎥ exp(−ξ 2 )d ξ (3.371) 2⎟ π x / 4α t ⎣ ⎝ 4αξ ⎠ ⎦
∫
which can be rewritten as 2A ϑ ( x, t ) =
⎡⎛ ⎤ x2 ⎞ cos ⎢ω ⎜ t − − β ⎥ exp(−ξ 2 )d ξ 2⎟ π0 ⎣ ⎝ 4αξ ⎠ ⎦ (3.372) 2 x / 4α t ⎡⎛ ⎤ 2A x⎞ 2 cos ⎢ω ⎜ t − − ⎟ − β ⎥ exp(−ξ )d ξ 4αξ 2 ⎠ π0 ⎣⎝ ⎦
∫
∞
∫
Chapter 3 Heat Conduction 263
Evaluating the first integral on the right hand side of eq. (3.372) yields 1/ 2 ⎡ ⎛ ω ⎞1/ 2 ⎤ ⎡ ⎤ ⎛ω⎞ ϑ ( x, t ) = A exp ⎢ − x ⎜ ⎟ ⎥ cos ⎢ωt − x ⎜ ⎟ − β ⎥ ⎝ 2α ⎠ ⎢ ⎝ 2α ⎠ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ (3.373) 2 ⎡⎛ ⎤ 2 A x / 4α t x⎞ 2 cos ⎢ω ⎜ t − − ⎟ − β ⎥ exp(−ξ )d ξ 0 4αξ 2 ⎠ π ⎣⎝ ⎦ It can be seen that as t → ∞ , the second term will become zero and the first term represents the steady oscillation. 1/ 2 ⎡ ⎛ ω ⎞1/ 2 ⎤ ⎡ ⎤ ⎛ω⎞ ϑs ( x, t ) = A exp ⎢ − x ⎜ ⎟ ⎥ cos ⎢ωt − x ⎜ ⎟ − β ⎥ (3.374) ⎝ 2α ⎠ ⎢ ⎝ 2α ⎠ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦
∫
1/ 2 where A exp ⎡ − x (ω /(2α ) ) ⎤ represents the amplitude of oscillation at point x, ⎣ ⎦
and − x (ω /(2α ) ) in the cosine function represents the phase delay of oscillation at the point x relative to the oscillation of the surface temperature. It is apparent that the amplitude of the oscillation decreases as x increases, and the phase of oscillation delays with increasing x. The oscillation at x is not as strong as that at the surface and there is a delay from the time that the surface temperature changes to the time that the temperature at x responds to such change.
1/ 2
Integral Approximate Solution 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 effect of sudden surface temperature change gradually propagates into the semi-infinite body. It is useful here to introduce a concept similar to the thermal boundary layer for convective heat transfer – 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. 3.20). 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 ) x = δ (t ) =0 (3.375) ∂x T ( x, t ) = Ti x = δ (t ) (3.376) Integrating eq. (3.332) in the interval ( 0, δ ), one obtains
∂T ∂x
−
x =δ ( t )
∂T ∂x
=
x=0
α∫
1
δ (t )
0
∂T ( x, t ) dx ∂t
(3.377)
264 Advanced Heat and Mass Transfer
T T0
Ti
δ(t)
x
Figure 3.20 Heat conduction in a semi-infinite body with constant wall temperature.
The right-hand side of eq. (3.377) can be rewritten using Leibnitz’s rule, i.e., ∂T ∂T 1⎡d ⎛ δ dδ ⎤ ⎞ − = ⎢ ⎜ Tdx ⎟ − T x =δ (3.378) 0 ∂x x =δ ( t ) ∂x x = 0 α ⎣ dt ⎝ dt ⎥ ⎠ ⎦
∫
which represents the energy balance within the thermal penetration depth. Substituting eqs. (3.375) and (3.376) into eq. (3.378) yields ∂T d −α = (Θ − Tiδ ) (3.379) ∂x x = 0 dt where Θ(t ) =
∫
δ (t )
0
T ( x, t ) dx
(3.380)
Equation (3.379) is the integral energy equation of the conduction problem, and this equation pertains to the entire thermal penetration depth. It follows that a temperature distribution that satisfies eq. (3.379) does not necessarily satisfy differential eq. (3.332), which describes the energy balance at any and all points in the domain of the problem. In order to use the integral equation (3.379) 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 (3.381) 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. (3.333), (3.375) and (3.376) – one more condition must be identified so that
Chapter 3 Heat Conduction 265
all four constants in eq. (3.381) 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 (3.382) ∂t Substituting eq. (3.332) into eq. (3.382) yields ∂ 2T ( x, t ) =0 (3.383) x=0 ∂x 2 Substituting eq. (3.381) into eqs. (3.333), (3.375), (3.376) and (3.383) yields four equations for the constants in eq. (3.381). Solving for the four constants and substituting the results into eq. (3.381), the temperature distribution in the thermal penetration depth becomes 3 T ( x, t ) − Ti 3⎛ x ⎞ 1⎛ x ⎞ =1− ⎜ ⎟ + ⎜ ⎟ (3.384) T0 − Ti 2⎝δ ⎠ 2⎝δ ⎠ where the thermal penetration depth, δ , is still unknown. Substituting eq. (3.384) into eq. (3.379), an ordinary differential equation for δ is obtained: dδ 4α = δ t >0 (3.385) dt Since the thermal penetration depth equals zero at the beginning of the heat conduction, eq. (3.385) is subject to the following initial condition: δ =0 t =0 (3.386) The solution of eqs. (3.385) and (3.386) is δ = 8α t (3.387) which is consistent with the result of scale analysis, δ ∼ α t [see eq. (1.186)]. The temperature distribution in the thermal penetration depth can be obtained by eq. (3.384), 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.
266 Advanced Heat and Mass Transfer
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. 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.
3.3.3 Multidimensional Transient Heat Conduction Systems
The transient heat conduction problems discussed in the preceding subsection are for the case that the transient temperature varies in one dimension only. Analytical solution of multidimensional transient heat conduction using product solution will be discussed in this subsection. More complicated multidimensional transient heat conduction can be solved using the numerical method to be discussed in the next section. Let us consider transient heat conduction in a rectangular bar with dimensions of 2L1×2L2 and an initial temperature of Ti (see Fig. 3.21). At time t = 0, the rectangular bar is immersed in a fluid with temperature T∞. The convective heat transfer coefficient on the left and right sides of the rectangular bar is h1, while on the top and bottom it is h2. It is assumed that the thermal conductivity is independent from the temperature and there is no internal heat source. Since this is an axisymmetric problem, we can study a quarter of it. The energy equation for this problem is ∂ 2T ∂ 2T 1 ∂T + = (3.388) , 0 < x < L1 , 0 < y < L2 , t > 0 ∂x 2 ∂y 2 α ∂t with the following boundary and initial conditions: ∂T = 0, x = 0, 0 < y < L2 , t > 0 (3.389) ∂x y L2 h 1 , T∞ – L1 0 h 2 , T∞ h 1 , T∞ L1 x
-L2
h 2 , T∞
Figure 3.21 Two-dimensional transient heat conduction
Chapter 3 Heat Conduction 267
∂T = h(T − T∞ ), x = L1 , 0 < y < L2 , t > 0 (3.390) ∂x ∂T = 0, y = 0, 0 < x < L1 , t > 0 (3.391) ∂y ∂T −k = h(T − T∞ ), y = L2 , 0 < x < L1 , t > 0 (3.392) ∂y T = Ti , 0 < x < L1 , 0 < y < L2 , t > 0 (3.393) Defining the following dimensionless variables T − T∞ x y αt αt , X = , Y = , Fo1 = 2 , Fo 2 = 2 (3.394) θ= Ti − T∞ L1 L2 L1 L2 where both Fo1 and Fo2 are Fourier numbers based on different characteristic lengths. The dimensional time will be related to both Fourier numbers, i.e., t = t (Fo1 , Fo 2 ) (3.395) therefore, the right-hand side of eq. (3.388) becomes ⎡ ∂θ ∂Fo1 ∂T ∂θ ∂θ ∂Fo 2 ⎤ = (Ti − T∞ ) = (Ti − T∞ ) ⎢ + ⎥ ∂t ∂t ∂Fo 2 ∂t ⎦ ⎣ ∂Fo1 ∂t (3.396) ⎡ ∂θ α ∂θ α ⎤ = (Ti − T∞ ) ⎢ + 2 2⎥ ⎣ ∂Fo1 L1 ∂Fo 2 L2 ⎦ −k
Nondimensionalizing the left-hand side of eq. (3.388) and considering eq. (3.396), the dimensionless energy equation of the problem becomes ∂ 2θ ⎛ L1 ⎞ ∂ 2θ ∂θ ⎛ L1 ⎞ ∂θ +⎜ ⎟ = +⎜ ⎟ 2 2 ∂X ∂Fo1 ⎝ L2 ⎠ ∂Fo 2 ⎝ L2 ⎠ ∂Y The boundary and initial conditions, eqs. (3.389) – nondimensionalized as ∂θ = 0, X = 0, 0 < Y < 1, Fo1 > 0, Fo 2 > 0 ∂X ∂θ − = Bi1θ , X = 1, 0 < Y < 1, Fo1 > 0, Fo 2 > 0 ∂X ∂θ = 0, Y = 0, 0 < X < 1, Fo1 > 0, Fo 2 > 0 ∂Y ∂θ − = Bi 2θ , Y = 1, 0 < X < 1, Fo1 > 0, Fo 2 > 0 ∂Y θ = 1, 0 < X < 1, 0 < Y < 1, Fo1 = Fo 2 = 0 where hL hL Bi1 = 1 1 , Bi 2 = 2 2 k k are Biot numbers for different surfaces.
2 2
(3.397) (3.393) are
(3.398) (3.399) (3.400) (3.401) (3.402) (3.403)
268 Advanced Heat and Mass Transfer
The idea of product solution is that the solution of the two-dimensional problem can be expressed as the product of two one-dimensional problems, i.e., θ = ϕ ( X , Fo1 )ψ (Y , Fo 2 ) (3.404) where ϕ is the solution of the following problem ∂ 2ϕ ∂ϕ = 2 ∂X ∂Fo1 ∂ϕ = 0, X = 0, Fo1 > 0 ∂X ∂ϕ − = Bi1ϕ , X = 1, Fo1 > 0 ∂X ϕ = 1, 0 < X < 1, Fo1 = 0 ∂ 2ψ ∂ψ = 2 ∂Y ∂Fo 2 (3.405) (3.406) (3.407) (3.408) (3.409)
and ψ satisfies
∂ψ = 0, Y = 0, Fo 2 > 0 (3.410) ∂Y ∂ψ (3.411) − = Bi 2ψ , Y = 1, Fo 2 > 0 ∂Y ψ = 1, 0 < Y < 1, Fo 2 = 0 (3.412) It can be demonstrated that eqs. (3.397) – (3.402) can be satisfied by eq. (3.404) if ϕ and ψ are solutions of eqs. (3.405) – (3.408) and (3.409) – (3.412), respectively. The solution of eqs. (3.405) – (3.408) is: ∞ 2 4sin λn (3.413) ϕ= cos ( λn X ) e − λn Fo1 n =1 2λn + sin 2λn where
∑
= cot λn Bi1 and the solution of eqs. (3.409) – (3.412) is ∞ 2 4sinν m ψ= cos (ν mY ) e −ν m Fo2 ν ν m =1 2 m + sin 2 m where
λn
(3.414)
∑
(3.415)
= cotν m (3.416) Bi 2 Substituting eqs. (3.413) and (3.415) into eq. (3.404), the solution of the twodimensional problem is obtained. ∞ ∞ 16sin λn sinν m cos ( λn X ) cos (ν mY ) − ( λn2 Fo1 +ν m Fo2 ) 2 (3.417) θ= e (2λn + sin 2λn )(2ν m + sin 2ν m ) n =1 m =1
νm
∑∑
Chapter 3 Heat Conduction 269
The product solution constructs the solution of a two-dimensional transient conduction as the product of the solutions of two one-dimensional problems. It is different from the method of separation of variables that changes the solution of a partial differential equation into the product of the solutions of ordinary differential equations. The product solution only works for transient problems and it will not work with steady-state problems. The only region that the product solution works is the intersection space of the two one-dimensional problems. The product solution will also work for a three-dimensional brick where the solution can be expressed as the solution of three one-dimensional problems (see Problem 3.42). Another example that product solution can be applied to is transient conduction in a short cylinder where the solution can be expressed as the product of solutions of one-dimensional transient conduction in a plane wall and a cylinder (Problem 3.43).
3.4 Numerical Simulation of Heat Conduction Problems
The steady and transient heat conductions discussed in the preceding two sections are for linear problems with simple geometric configurations and boundary conditions. Even for the cases where analytical solution is available, the interpretation of the solution is often difficult because the solution is expressed as a very complicated series. Numerical solution therefore becomes a desirable approach for heat conduction under non-linear, complex geometric configurations, or complex boundary conditions. Instead of obtaining the analytical expression of the temperature distribution, the results of numerical solutions are given at discrete points. The numerical solution involves three steps: (1) discretization of the computational domain, (2) discretization of the governing equations, and (3) solution of the algebraic equations (Murthy et al., 2006). The discretization of the computational domain is a process that divides the computational domain into many control volumes as shown in Fig. 3.22 (Patankar, 1980). Each subdomain is bounded by faces represented by the dashed lines. The physical properties of the cell are defined at the grid points. In Practice A, the faces are always located midway between the grid points. On the contrary, the grid points are always located in the center of the control volume in Practice B. Obviously, if a uniform grid is used, the two practices result in identical grid size and therefore, any discussion on the differences between two practices are meaningful only if the grid size is not uniform. Since the faces in Practice A are located in the midway between grid points, the heat flux across the faces can be more accurately calculated. The disadvantage of Practice A is that the properties of the entire control volume are represented by a point not located in the center of the control volume, which will result in inaccuracy. On the other hand, the properties of the entire control volume in Practice B are represented by the grid point at the center which is a better representation. In addition, Practice
270 Advanced Heat and Mass Transfer
(a) Practice A
(b) Practice B
Figure 3.22 Discritization of the computational domain.(Patankar, 1980)
B can easily handle discontinuity of thermophysical properties or boundary conditions. The discretization of governing equations can be done by local or point-wise representation of the partial differential equations (finite difference method; FDM), or weighted integral of the partial differential equations (finite element method; FEM). Patankar (1980) proposed a finite volume method (FVM) involving obtaining the discretized equation by performing integration on the governing equation over the small region. While the resultant algebraic equations for FVM and FDM are often similar, the FVM can guarantee conservation of the mass, momentum and energy on each cell, regardless of the size of the cell. The FVM can also be very easily extended to convective heat transfer because mature numerical methods have been developed in the last three decades. Other notable methods for numerical solution of heat conduction problems include boundary element method (BEM) that volume integrals are converted to surface integrals, as well as spectral methods that approximate solutions by a series of functions. The drawback of these methods is that they are very difficult to extend to convective heat transfer. To maintain uniformity of the numerical method for all heat transfer modes and build a solid foundation for students to learn other methods, we will focus on FVM. This section will cover numerical solutions of one-dimensional steady and transient conduction problems, as well as multidimensional transient heat conduction problems. Its extension to convection and radiation problems will be discussed in detail in later chapters of this text.
Chapter 3 Heat Conduction 271
3.4.1 One-Dimensional Steady-State Conduction
Discretization of Governing Equations
To keep our discussion as general as possible, let us consider a one-dimensional steady-state heat conduction problem with temperature-dependent thermal conductivity and internal heat generation (Tao, 2001): 1 d⎛ dT ⎞ (3.418) ⎜ kA( x) ⎟+S =0 A( x) dx ⎝ dx ⎠ where A(x) is the area of heat conduction. Equation (3.418) is valid for heat conduction in Cartesian ( A( x) = 1 ), cylindrical ( x = r , A(r ) = r ), and spherical ( x = r , A(r ) = r 2 ) coordinate systems, as well as heat transfer from extended surfaces (A(x) is the cross-sectional area of the fin). The source term, S, can be caused by internal heat source/sink, convection from the extended surface, or metabolic heat generation and blood perfusion in bioheat transfer problems. The source term in eq. (3.418) is often a function of temperature as in the cases of heat transfer from extended surfaces and bioheat transfer problems. It is desirable to reflect this temperature-dependence in our discretized governing equation. While the dependence of the source term on the temperature can be of a very complicated form (e.g., radiation heat loss from the extended surface), it is ideal to express the source term as a linear function of temperature because the discretized equation will be a linear algebraic equation. We can thus express the source term as S = SC + S P T (3.419) where SC is a constant and SP is a coefficient. Substituting eq. (3.419) into eq. (3.418), and multiplying the resulting equation by A(x) yields d⎛ dT ⎞ ⎜ kA( x) ⎟ + ( SC + S PT ) A( x) = 0 dx ⎝ dx ⎠ Integrating the above equation over the control volume P (shaded area in Fig. 3.23), one obtains e ⎛ dT ⎞ ⎛ dT ⎞ kA − ⎜ kA + ( SC + S PT ) Adx = 0 (3.420) ⎜ ⎟ ⎟ w ⎝ dx ⎠e ⎝ dx ⎠ w Assuming the temperature T and conduction area A for each control volume are represented by their value at a grid point, the integral in eq. (3.420) becomes
∫
∫ (S
w
e
C
+ S PT ) Adx = SC AP (Δx) P + S PTP AP (Δx) P
Assuming the temperature distribution between any two neighboring grid points is piecewise linear, the first two terms on the left-hand side of eq. (3.420) become TE − TP ⎛ dT ⎞ ⎜ kA ⎟ = ke Ae (δ x)e ⎝ dx ⎠e
272 Advanced Heat and Mass Transfer
(δx)w (δx)ww (δx)w+
(δx)e (δx)ee P (Δx)P E (δx)e+
Control Volume W
Figure 3.23 Control volume for one-dimensional heat conduction
TP − TW ⎛ dT ⎞ ⎜ kA ⎟ = k w Aw (δ x) w ⎝ dx ⎠ w where ke and kw are thermal conductivities at the faces of the control volume. To ensure that the heat flux across the faces of the control volume is continuous, it is important that the following harmonic mean conductivity at the faces are used (see Problem 3.43), i.e., (δ x)e (δ x)e − (δ x)e + (δ x) w (δ x) w− (δ x) w+ = + = + (3.421) , ke kP kE kw kW kP For a uniform grid system, eq. (3.421) becomes 2k K 2k P K E ke = , kw = W P (3.422) kP + K E kW + K P Substituting the above equations into eq. (3.420), the following discretized equation is obtained aPTP = aE TE + aW TW + b (3.423) where kA kA (3.424) aE = e e , aW = w w (δ x)e (δ x) w aP = aW + aE − S P AP (Δx) P (3.425) b = SC AP (Δx) P (3.426)
Boundary Conditions
If the temperature at the boundary is known (boundary condition of the first kind), no specific treatment on the boundary condition is necessary and the temperature at internal grid points can be obtained by solving the above discretized equations. For boundary conditions of the second and third kind, additional equations are necessary for the boundary temperatures. For Practice A, the width of the control volume at the boundary is equal to half of that of the
Chapter 3 Heat Conduction 273
(δx)1
(δx)M
1 (Δx)1
2
M-1 (Δx)M
M
(a) Left boundary (b) Right boundary Figure 3.24 Control volumes at boundaries for Practice A
internal control volumes (see Fig. 3.24). If the heat flux at the left boundary is known, the energy balance for the control volume at the left boundary is T −T ′′ qB A1 + k1 A1 2 1 + ( Sc + S PT1 ) A1 (Δx)1 = 0 (δ x)1 which can be expressed as
a1T1 = aET2 + b
where k1 A1 (δ x)1 aP = aE − S P A1 (Δx)1 ′′ b = qB A1 + SC A1 (Δx)1 If the left boundary is subject to a convective condition, we have ′′ qB = h(T f − T1 ) aE =
(3.427) (3.428) (3.429) (3.430) (3.431)
It can be shown that eqs. (3.427) and (3.428) are still valid but eqs. (3.429) and (3.430) should be modified to aP = aE − S P A1 (Δx)1 + hA1 (3.432) b = hT f A1 + SC A1 (Δx)1 (3.433) The boundary condition on the right hand side can be obtained by following a similar procedure (see Problem 3.46). If the computational domain is discretized by using Practice B, the size of the control volume for the grid at the boundary is zero (see Fig. 3.25), thus the discretized equation for the boundary condition at the left side can be obtained by setting the control volume size (Δx)1 to zero in eqs. (3.429) – (3.430) or eqs. (3.432) – (3.433). Similarly, the discretized equation for the boundary condition on the right side for Practice B can be obtained by setting (Δx) M in the discretized equation for Practice A.
274 Advanced Heat and Mass Transfer
(δx)1
(δx)M
1
2 (Δx)2
M-1 (Δx)M-1
M
(a) Left boundary (b) Right boundary Figure 3.25 Control volumes at boundaries for Practice B
Solution of Discretized Equations
When the matrix for the coefficient for discretized algebraic equations of onedimensional conduction is written, all of the nonzero components are aligned along three diagonals of the matrix. While the discretized algebraic equations can be solved by using Gaussian elimination method, a simple version – referred to as TDMA (TriDiagonal Matrix Algorithm) – can be used (Patankar, 1980). The discretized equations for every point, eq. (3.423), can be rewritten as aiTi = biTi +1 + ciTi −1 + d i , i = 1, 2,...M (3.434) which indicates that the temperature at any point is related to temperatures of its immediate neighbor points only. Obviously, for the left and boundaries, we have c1 = 0, bM = 0 To get the temperature, we hope that our temperature at the ith grid is related to the neighbor on the right, i.e., Ti = PTi +1 + Qi (3.435) i th For the (i − 1) grid point, eq. (3.435) becomes Ti −1 = Pi −1Ti + Qi −1 (3.436) Substituting eq. (3.436) into eq. (3.434), we get aiTi = biTi +1 + ci ( Pi −1Ti + Qi −1 ) + bi (3.437) which can be rearranged to bi d +cQ Ti = Ti +1 + i i i −1 (3.438) ai − ci Pi −1 ai − ci Pi −1 Comparing eq. (3.435) and (3.438) yields bi d +cQ Pi = , Qi = i i i −1 (3.439) ai − ci Pi −1 ai − ci Pi −1 For the first grid point at the left (i = 1), c1 = 0 and eq. (3.439) becomes b d P = 1 , Q1 = 1 (3.440) 1 a1 a1
Chapter 3 Heat Conduction 275
For the last point at the right (i = M), bM = 0 and it follows that PM = 0 . Equation (3.435) for i = M becomes TN = QN (3.441) The solution of the algebraic equations can thus be completed in the following steps: (1) obtain P1 and Q1 from eq. (3.440), (2) calculate Pi and Qi from eq. (3.439), (3) obtain temperature at the right side by eq. (3.441), and (4) calculate temperatures for all grid points by eq. (3.435). It should be noted that if the thermal conductivity is independent from temperature, the temperatures at different grid points can be obtained by solving the discretized equations using TDMA once. If the thermal conductivity depends on temperature, the problem needs to be solved by an iteration procedure: (1) assume tentative thermal conductivity at all grid points, (2) obtain coefficients for discretized equation based on the tentative conductivity, (3) solve the discretized equation using TDMA, and (4) update the thermal conductivities at all points and solve for the temperatures. This procedure needs to be repeated until there is no notable change of temperature distributions between two consecutive iterations.
3.4.2 One-Dimensional Transient Conduction
Discretization of Governing Equations
For transient conduction problems, eq. (3.418) should be modified to: 1 d⎛ dT ⎞ ∂T (3.442) ⎜ kA( x) ⎟ + S = ρc A( x) dx ⎝ dx ⎠ ∂t Similar to the case of steady-state conduction, it is assumed that the source term is a linear function of temperature. Substituting eq. (3.419) into eq. (3.442) and multiplying the resulting equation by A(x) yields d⎛ dT ⎞ ∂T ⎜ kA( x) ⎟ + ( SC + S PT ) A( x) = ρ cA( x) ∂t dx ⎝ dx ⎠ Integrating the above equation with respect to t in the interval of (t, t+Δt) and with respect to x over the control volume P (shaded area in Fig. 3.19), we have t +Δt w d t +Δt w dT ⎞ ⎛ ( SC + S PT ) A( x)dxdt ⎜ kA( x) ⎟dxdt + t t e dx ⎝ e dx ⎠ (3.443) w t +Δt ∂T ρ cA( x) dtdx = e t ∂t The first integral on the left-hand side becomes t +Δt w d t +Δt ⎡ dT ⎞ ⎛ ⎛ dT ⎞ ⎛ dT ⎞ ⎤ kA( x) dxdt = ⎢⎜ kA ⎜ ⎟ ⎟ − ⎜ kA ⎟ ⎥dt t e dx ⎝ t dx ⎠ ⎣⎝ dx ⎠e ⎝ dx ⎠ w ⎦
∫∫
∫∫
∫∫
∫∫ ∫
t
∫
=
t +Δt
⎡ T −T ⎤ TE − TP − k w Aw P W ⎥dt ⎢ ke Ae (δ x)e (δ x) w ⎦ ⎣
276 Advanced Heat and Mass Transfer
TP TPt + Δt TPt f=1 f = 0.5 f=0
t
t+Δt
Figure 3.26 Variation of temperature within one time step
which requires knowledge of how temperature varies in the time interval of (t, t+Δt). There are commonly three kinds of choices (see Fig. 3.26): (1) TP (t ) ≡ TPt where TPt is the temperature of grid point P at time t, (2) TP (t ) ≡ TPt +Δt where
TPt +Δt is the temperature of grid point P at time t+Δt, or (3) TP (t ) linearly varies
from TPt at time t to TPt + Δt at time t+Δt. The following result is valid for all three variations: where the above three cases correspond to f = 0 , f = 1 , and f = 1/ 2 . Similar expressions for TE and TW can be obtained and we have t +Δt w d ⎧⎡ T −T ⎤ TE − TP dT ⎞ ⎪ ⎛ − kw Aw P W ⎥ ⎜ kA( x) ⎟dxdt = ⎨ f ⎢ ke Ae t e dx ⎝ dx ⎠ (δ x)e (δ x) w ⎦ ⎪⎣ ⎩
∫
t +Δt
t
TP (t )dt = ⎡ fTPt +Δt + (1 − f )TPt ⎤ Δt ⎣ ⎦
∫∫
0 ⎫ ⎡ T 0 − TW ⎤ ⎪ T 0 − TP0 + (1 − f ) ⎢ ke Ae E − k w Aw P ⎥ ⎬ Δt (δ x)e (δ x) w ⎦ ⎪ ⎣ ⎭ where the superscript t is replaced by 0 and the superscript t+Δt has been dropped, i.e., temperatures for grid point P at time t and t+Δt are represented by TP0 and at time t+Δt is represented by TP, respectively. The second integral on the left-hand side of eq. (3.443) can be expressed as
∫∫
t
t +Δt
w
e
( SC + S PT ) A( x)dxdt = {SC + S P [ fTP0 + (1 − f )TP ]} AP (Δx) P Δt
The integral on the right hand of eq. (3.443) can be evaluated as w t +Δt ∂T ρ cA( x) dtdx = ( ρ c) P AP (TP − TP0 )(Δx) P e t ∂t Substituting the above equations into eq. (3.443) the discretized equation for onedimensional transient conduction becomes
∫∫
Chapter 3 Heat Conduction 277
+ {SC + S P [ fTP0 + (1 − f )TP ]} AP (Δx) P Δt = ( ρ c) P AP (TP − TP0 )(Δx) P (3.444) which can be rearranged to obtain (Tao, 2001) 0 aPTP = aE [ fTE + (1 − f )TE0 ] + aW [ fTW + (1 − f )TW ] + b where kA kA aE = e e , aW = w w (δ x)e (δ x) w
0 aP = faE + faW + aP − fS P AP (Δx) P
⎧ ⎪ ⎨f ⎪ ⎩
⎡ T −T TE − TP − k w Aw P W ⎢ ke Ae (δ x)e (δ x) w ⎣
0 ⎤ ⎡ T 0 − TW ⎤ ⎫ T 0 − TP0 ⎪ + (1 − f ) ⎢ ke Ae E − k w Aw P ⎥ ⎥ ⎬ Δt (δ x)e (δ x) w ⎦ ⎪ ⎦ ⎣ ⎭
(3.445) (3.446) (3.447) (3.448) (3.449)
b = ⎡ a − (1 − f ) aE − (1 − f )aW + (1 − f ) S P AP (Δx) P ⎤ T + SC AP (Δx) P ⎣ ⎦
0 P 0 P 0 aP =
ρ c ( Δx ) P
Δt
Explicit and Implicit Schemes
Different schemes can be obtained when f in eqs. (3.445) – (3.449) takes different values. When f = 0 , i.e., the temperature at grid point P in the entire interval of (t, t+Δt) is equal to the temperature at time t, eq. (3.445) becomes 0 aPTP = aE TE0 + aW TW + b (3.450) which indicates that the temperature for grid point P at time t + Δt can be explicitly obtained from temperatures at the previous time step, because it is not related to temperatures of grid points E and W at time t+Δt. The scheme for the case of f = 0 is referred to as the explicit scheme. While an explicit scheme is computationally very convenient, it is not unconditionally stable. For the case of constant thermophysical properties, without an internal heat source, and with uniform grid size, eq. (3.450) becomes ⎡ 2αΔt ⎤ 0 α Δt 0 0 TP = (TE + TW ) + ⎢1 − T (3.451) 2 2⎥ P ( Δx) ⎣ (Δx) ⎦ Although the criterion of the stability for eq. (3.451) can be rigorously derived from Von Neumann analysis, we will use a simple reasoning to obtain this criterion. If the temperatures at E and W remain unchanged, we would expect that TP is higher if TP0 is higher and TP is lower if TP0 is lower. This can be assured only if the coefficient of TP0 in eq. (3.451) is positive, i.e., α Δt 1 < Fo Δ = (3.452) (Δx) 2 2 where FoΔ is grid Fourier number. Equation (3.452) indicates that for a system with uniform grid size of Δx, the time step Δt must be small enough to ensure the grid Fourier number is less than 0.5. Equation (3.451) can be rewritten in terms of the grid Fourier’s number as
278 Advanced Heat and Mass Transfer
0 TP = Fo Δ (TE0 + TW ) + (1 − 2Fo Δ )TP0 (3.453) If f = 0.5 , i.e., the temperature at grid point P varies linearly in the interval of (t, t+Δt). Equation (3.445) becomes a a 0 (3.454) aPTP = E [TE + TE0 ] + W [TW + TW ] + b 2 2 which is referred to as the Crank-Nicolson (C-N) scheme. Although the C-N scheme has often been described as unconditionally stable, its result is physically meaningful only if the grid Fourier number is less than 1. For the case that Fo Δ > 1 , C-N scheme is mathematically stable but the solution physically does not make sense. When f = 1 , eqs. (3.445), (3.447) and (3.448) become aPTP = aE TE + aW TW + b (3.455)
0 aP = aE + aW + aP − S P AP (Δx) P
(3.456) (3.457)
b = ⎡ a + S P AP (Δx) P ⎤ + SC AP (Δx) P ⎣ ⎦
0 P
which is referred to as a fully-implicit scheme. It is truly unconditionally stable and physically meaningful under any grid Fourier number. This advantage makes the fully-implicit scheme top choice for numerical solution of conduction problems. Since fully-implicit schemes are stable under any time step, one can set Δt → ∞ and the discretized equations for transient conduction become identical to those for the steady-state conduction problem. Thus, we can run the program for transient heat conduction with implicit schemes for very large time steps, say 1020, to get a solution for a steady state conduction problem. Since the temperature of grid point P at time t+Δt is related to the temperatures of its neighbor points E and W, the TDMA discussed in the previous subsection can be employed.
3.4.3 Multidimensional Transient Conduction
Two-Dimensional Transient Conduction
The energy equation for a two-dimensional transient conduction problem with temperature-dependent thermal conductivity and internal heat source can be written as ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂T (3.458) ⎟ + S = ρc ⎜k ⎟ + ⎜k ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂t To discretize eq. (3.458), integrating eq. (3.458) with respect to t in the interval of (t, t+Δt) and over the control volume P (shaded area in Fig. 3.27) gives (Tao, 2001)
Chapter 3 Heat Conduction 279
N (δx)w n Δy W w Δx Control Volume P s S
Figure 3.27 Control volume for two-dimensional heat conduction
(δx)e
(δy)n
e
E (δy)s
∫ ∫∫
t s
t +Δt
n
∂ ⎛ ∂T ⎜k w ∂x ⎝ ∂x
e e w
⎞ ⎟ dxdydt + ⎠
∫ ∫∫
t e
t +Δt
w
n
s
∂ ⎛ ∂T ⎜k ∂y ⎝ ∂y
⎞ ⎟ dydxdt ⎠
∂T dtdxdy ∂t Assuming the heat fluxes are uniform on all faces of the control volume and employing a fully-implicit scheme, the above equation becomes ⎡ TE − TP ⎡ T − TP T −T ⎤ T − TS ⎤ − k w P W ⎥ ΔyΔt + ⎢ kn N − ks P ⎢ ke ⎥ ΔxΔt (δ x) w ⎦ (δ y ) s ⎦ ⎣ (δ x)e ⎣ (δ y ) n +
∫ ∫∫
t s
t +Δt
n
( SC + S PT )dxdydt =
∫∫∫
s e
n
w
t +Δt
t
ρc
+ ( SC + S PT )ΔxΔyΔt = ( ρ c) P (TP − TP0 )ΔxΔy which can be rearranged to obtain (Tao, 2001) aPTP = aE TE + aW TW + aN TN + aS TS + b where k Δy k Δy k Δx k Δx aE = e , aW = w , aN = n , aS = s (δ x)e (δ x) w (δ y ) n (δ y ) s
0 aP = aE + aW + aN + aS + aP − S P ΔxΔy
(3.459) (3.460) (3.461) (3.462) (3.463)
b = a T + SC ΔxΔy ( ρ c ) P Δx Δ y 0 aP = Δt
00 PP
Three-Dimensional Transient Conduction
The discretized equation for a three-dimensional transient conduction problem with temperature-dependent thermal conductivity and an internal heat
280 Advanced Heat and Mass Transfer
source can be obtained by integrating the energy equation with respect to t in the interval of (t, t+Δt) and over the three-dimensional control volume P (formed by considering two additional neighbors at top, T, and bottom, B). The final form of the discretized equation is (Patankar, 1980) aPTP = aETE + aW TW + aN TN + aS TS + aT TT + aBTB + b (3.464) where k ΔyΔz k ΔyΔz k ΔxΔz aE = e , aW = w , aN = n (δ x)e (δ x) w (δ y ) n (3.465) k s ΔxΔz k t Δx Δ y kb ΔxΔy , aT = , aB = aS = (δ y ) s (δ z )t (δ z )b
0 aP = aE + aW + aN + aS + aT + aB + aP − S P ΔxΔyΔz
(3.466) (3.467) (3.468)
b = a T + S C Δ x Δ y Δz ( ρ c ) P Δx Δ y Δ z 0 aP = Δt
00 PP
Solution of Algebraic Equations
To obtain the temperature distribution for one time step, it is necessary to solve the discretized algebraic equations for all grid points simultaneously. Direct solution of these equations (e.g., using Cramer’s rule) requires a large amount of computational time and storage. If the problem is linear, we only need to solve the set of linear algebraic equations once for every step. If, however, the problem is nonlinear, solving the algebraic equations once only gives us the temperature distribution based on the temporary coefficients of the algebraic equations and several iterations must be performed to get the temperature distribution for one time step. For this case, it is imperative to develop an efficient algorithm to solve the algebraic equations. A simple but efficient algorithm – line-by-line method – that combines direct method (TDMA) and Gauss-Seidel iteration method will be discussed here (Patankar, 1980). We will demonstrate this method using a two-dimensional conduction problem. The discretized equation (3.459) can be rewritten as Ai , jTi , j = BijTi +1, j + Ci , jTi −1, j + Di , jTi , j +1 + Ei , jTi , j −1 + Ei , j (3.469) where i = 1, 2,...M , j = 1, 2,...N . The line-by-line method can be performed either in the x- or y-direction. To apply this method in the x-direction, we can choose the jth row to apply the TDMA [see Fig. 3.28(a)]. To do that, it is assumed that the temperature at the ( j − 1)th and ( j + 1)th row are known from the previous iteration. For the first iteration at the current time step, they can be taken as values at the previous time step. Equation (3.469) can then be rewritten as Ai , jTi , j = BijTi +1, j + Ci , jTi −1, j + ( Di , jTi , j +1 + Ei , jTi , j −1 + Ei , j ) (3.470) where ~ represents the values from the previous iteration. If we treat all terms in the parentheses on the right hand of eq. (3.470) as the constant term, di, in eq.
Chapter 3 Heat Conduction 281
N j+1 j j–1 3 2 i =1 2 3 M 1 i–1 i i+1
Figure 3.28 Line-by-line method
(3.434), eq. (3.470) can be solved using TDMA to get Ti,j in the jth row. We can perform this procedure from j = 1 to N to get the temperature for the whole domain updated. This can also be done in the reverse direction from j = N to 1 . The information about the boundary conditions in the x-direction can be propagated into the entire domain, but the information on the boundary condition in the y-direction can only propagate one grid into the computational domain. To speed up the process, we can also choose the ith column to update its temperature using TDMA [see Fig. 3.28(b)]. In this case, eq. (3.469) can then be rewritten as Ai , jTi , j = Di , jTi , j +1 + Ei , jTi , j −1 + ( BijTi +1, j + Ci , jTi −1, j + Ei , j ) (3.471) which can be applied either from i = 1 to M or i = M to 1 . The line-by-line method is a combination of direct method (TDMA) and iteration method. Once we complete the procedure along the x- direction twice (from j = 1 to N and then from j = N to 1 ) and along the y-direction twice (from i = 1 to M and then from i = M to 1 ) – referred to as one iteration – the information from boundary conditions on all four sides can be propagated throughout the computational domain. The numerical practices by the authors and other researchers indicate that a converged solution can be obtained after two to four iterations. The above line-by-line method can be easily extended to the case of threedimensional heat conduction problem. While our discussion here is for solution of algebraic equations resulting from the discretization of the heat conduction problem, it can also be used to solve the discretized equations for convection problems because – as will become evident in the next two chapters – the discretized equation for convection problems will be identical to eqs. (3.459) or (3.464).
282 Advanced Heat and Mass Transfer
3.5 Melting and Solidification
3.5.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. Figure 3.29 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. T T0 T
Ti , Tm
Ti , Tm
T0 s(t)
(a) Melting ( T0 > Tm )
Figure 3.29 One-region melting and solidification.
x
s(t)
(b) Solidification ( T0 < Tm )
x
Chapter 3 Heat Conduction 283
T T0 Tm Ti
T Ti Tm
T0 s(t) (a) Melting ( T0 > Tm > Ti )
Figure 3.30 Two-region melting and solidification.
x
s(t) (b) Solidification ( T0 < Tm < Ti )
x
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. Figure 3.30 shows the temperature distribution of onedimensional two-region melting and solidification problems. 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. 3.31, 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 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. 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. As shown for the one-dimensional melting problem illustrated in Fig. 3.30, the solid-liquid interface separates the liquid and solid phases. The temperatures of the liquid and solid phases near the interface must equal the
284 Advanced Heat and Mass Transfer
T Ti Tm2 Tm1
x
Figure 3.31 Solidification of a multicomponent PCM.
temperature of the interface, which is at the melting point, T m . Therefore, the boundary conditions at the interface can be expressed as T ( x, t ) = Ts ( x, t ) = Tm , x = s (t ) (3.472) where T ( 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 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, ρ , 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. 3.30(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 > ρ ) 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 ρ ) = ρ u p
(3.473)
Chapter 3 Heat Conduction 285
The extra velocity induced by the density change can be obtained by rearranging eq. (3.473) as ρ −ρ uρ = s up (3.474)
ρs
which is also valid for case ρ s > ρ , except the extra velocity becomes negative. Another necessary boundary condition is the energy balance at the solidliquid interface. If the enthalpy of the liquid and solid phases at the melting point are h and hs , the energy balance at the solid-liquid interface can be expressed as: ′′ q′′ − qs = ρ u p h − ρ s (u p − uρ )hs x = s (t ) (3.475)
′′ where q′′ and qs are the heat fluxes in the x-direction in the liquid and solid phases, respectively. Substituting eq. (3.474) into eq. (3.475), one obtains ds ′′ (3.476) q′′ − qs = ρ hs x = s (t ) dt where hs = h − 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 ) q′′ = − k (3.477) ∂x x = s ( t )
′′ qs = − k s ∂Ts ( x, t ) ∂x (3.478)
x = s (t )
The energy balance at the solid-liquid interface for a melting problem can be obtained by substituting eqs. (3.477) and (3.478) into eq. (3.476), i.e., ∂T ( x, t ) ∂T ( x, t ) ds (t ) ks s x = s (t ) (3.479) −k = ρ hs dt ∂x ∂x 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 ) −k = ρ s hs ks s x = s (t ) (3.480) ∂x ∂x dt The only difference between eqs. (3.479) or (3.480) is the density on the right-hand 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. (3.479) or (3.480). It should be noted that density change causes advection in the liquid phase, which further complicates the problem. The above boundary conditions are valid for one-dimensional problems only. For multi-dimensional phase change problems, the boundary conditions at the interface can be expressed as T ( x, y, t ) = Ts ( x, y, t ) = Tm x = s ( y, t ) (3.481)
286 Advanced Heat and Mass Transfer
∂Ts ( x, y, t ) ∂T ( x, y, t ) −k = ρ hs vn x = s ( y, t ) (3.482) ∂n ∂n where n is a unit vector along the normal direction of the solid-liquid interface, and vn is the solid-liquid interface velocity along the n-direction. It is apparent that eq. (3.482) 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 ) Equation (3.482) can then become the following form (see Problem 3.55; Ozisik, 1993): ⎡ ⎛ ∂s ⎞ 2 ⎤ ⎡ ∂T ∂T ⎤ ∂s x = s ( y, t ) (3.483) ⎢1 + ⎜ ⎟ ⎥ ⎢ k s s − k ⎥ = ρ hs ∂t ∂x ⎦ ⎢ ⎝ ∂y ⎠ ⎥ ⎣ ∂x ⎣ ⎦ Similarly, for a three-dimensional melting problem with an interface described by z = s ( x, y , t ) the energy balance at the interface is ⎡ ⎛ ∂s ⎞ 2 ⎛ ∂s ⎞ 2 ⎤ ⎡ ∂T ∂T ⎤ ∂s = ρ hs z = s ( x, y, t ) (3.484) ⎢1 + ⎜ ⎟ + ⎜ ⎟ ⎥ ⎢ k s s − k ∂x ⎠ ⎝ ∂y ⎠ ⎥ ⎣ ∂x ∂x ⎥ ∂t ⎦ ⎢⎝ ⎣ ⎦ For solidification problems, it is necessary to replace the liquid density ρ in eqs. (3.483) and (3.484) 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. This section 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. ks
3.5.2 Exact Solution
Governing Equations of the Solidification Problem
The physical model of the solidification problem to be investigated in this subsection is shown in Fig. 3.30(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
Chapter 3 Heat Conduction 287
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 = (3.485) ∂x 2 α1 ∂t T1 ( x, t ) = T0 x = 0, t > 0 (3.486) For the liquid phase, the governing equations are ∂ 2T2 1 ∂T2 = s (t ) < x < ∞, t > 0 (3.487) 2 ∂x α 2 ∂t T2 ( x, t ) → Ti x → ∞, t > 0 (3.488) T2 ( x, t ) = Ti x > 0, t = 0 (3.489) The boundary conditions at the interface are T1 ( x, t ) = T2 ( x, t ) = Tm x = s (t ), t > 0 (3.490) ∂T ∂T ds (3.491) k1 1 − k2 2 = ρ hs x = s (t ), t > 0 ∂x ∂x dt Before obtaining the solution of the above problem, a scale analysis of the energy balance eq. (3.491) is performed. At the solid-liquid interface, the scales of derivatives of solid and liquid temperature are ∂T1 Tm − T0 ∼ ∂x s ∂T2 Ti − Tm ∼ ∂x s The scale of the interface velocity is ds s ∼ dt t Substituting the above equations into eq. (3.491), one obtains T −T T −T s k1 m 0 − k2 i m ∼ ρ hs s s t The scale of the location of the solid-liquid interface is obtained by rearranging the above equation, i.e., ⎛ k T −T ⎞ s2 ∼ Ste ⎜1 − 2 i m ⎟ (3.492) α1t ⎝ k1 Tm − T0 ⎠ where Ste = (3.493) hs 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, hs . For c p1 (Tm − T0 )
288 Advanced Heat and Mass Transfer
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 hs 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. (3.492) 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. (3.492) can be simplified as s2 ∼ Ste (3.494) α1t which means that the interfacial velocity increases with increasing ΔT = Tm − T0 or decreasing hs .
Dimensionless Form of the Governing Equations
The governing eqs. (3.485) – (3.491) 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⎪ ⎪ ⎬ (3.495) ⎪ c p1 (Tm − T0 ) k2 α2 ⎪ Nα = Nk = Ste = ⎪ ⎪ k1 hs α1 ⎭ where L 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 (3.496) ∂X 2 ∂τ θ1 ( X ,τ ) = 1 X = 0, τ > 0 (3.497)
∂ 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τ (3.498) (3.499) (3.500) (3.501) (3.502)
Chapter 3 Heat Conduction 289
θ
1.0
θ1(x,τ)
s(τ)
θi
θ2(x,τ)
x
Figure 3.32 Dimensionless temperature distribution in the PCM.
Dimensionless temperature distribution in a PCM can be qualitively illustrated by Fig. 3.32. It can be seen that the dimensionless temperature distribution is similar to that of a melting problem. Equations (3.496) – (3.502) 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.
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. (3.496) – (3.502), i.e., ∂ 2θ ∂θ = 0 < X < S (τ ), τ > 0 (3.503) ∂X 2 ∂τ θ ( X ,τ ) = 1 X = 0, τ > 0 (3.504) θ ( X ,τ ) = 0 X = S (τ ), τ > 0 (3.505) ∂θ 1 dS − = (3.506) X = S (τ ), τ > 0 ∂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 ) (3.507)
290 Advanced Heat and Mass Transfer
where B is an unspecified constant. It can be verified that eq. (3.507) satisfies eqs. (3.503) and (3.504). The constant B in eq. (3.507) can be determined by substituting eq. (3.507) into eq. (3.505), i.e., 0 = 1 + Berf ( S / 2τ 1/ 2 ) Since B is a constant, S / 2τ 1/ 2 must also be a constant in order for the above equation to be satisfied. This constant can be represented by λ , so λ = S / 2τ 1/ 2 (3.508) Thus, 1 B=− erf (λ ) By substituting the above expression of B into eq. (3.507), the dimensionless temperature distribution in phase 1 is obtained: erf ( X / 2τ 1/ 2 ) θ ( X ,τ ) = 1 − (3.509) erf (λ ) Substituting eq. (3.509) and eq. (3.508) into eq. (3.506), the following equation is obtained for the constant λ : 2 Ste (3.510) λ eλ erf (λ ) = 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. (3.509) and (3.508) respectively. The exact solution of the one-region phase change problem can also be obtained by using a similarity solution (see Problem 3.55).
π
Example 3.4 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. 3.33. 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 melting problem. The governing equation, and the corresponding initial and boundary conditions of the problem are as follows
Chapter 3 Heat Conduction 291
T
′ q0′
T(x,t)
Ti=Tm
s(t)
Figure 3.33 Melting in a semi-infinite region under constant heat flux
x
∂ 2T 1 ∂T = 0 < x < s (t ), t > 0 (3.511) ∂x 2 α ∂t ∂T ( x, t ) ′′ x = 0, t > 0 −k = q0 (3.512) ∂x T ( x, t ) = Tm x = s (t ), t > 0 (3.513) ∂T ds x = s (t ), t > 0 −k = ρ hs (3.514) dt ∂x 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 , X = , S = , τ = 2 , Ste = θ= (3.515) ′′ q0 L/k L L L hs where L is a characteristic length of the problem, eqs. (3.511) – (3.514) are nondimensionalized as ∂ 2θ ∂θ = (3.516) 0 < X < S (τ ), τ > 0 ∂X 2 ∂τ ∂θ ( X ,τ ) = −1 (3.517) X = 0, τ > 0 ∂X X = S (τ ), τ > 0 θ ( X ,τ ) = 0 (3.518) ∂θ 1 dS X = S (τ ), τ > 0 − = (3.519) ∂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:
292 Advanced Heat and Mass Transfer
⎛X⎞ ⎟ ⎝2 τ ⎠ where ierfc is a differential complementary error function: 1 − z2 ierfc( Z ) = e − Zerfc( Z )
θ ( X ,τ ) = A + 2 τ ierfc ⎜
(3.520)
π
(3.521)
The unknown variable A in eq. (3.520) can be obtained by substituting eq. (3.520) into eq. (3.518), i.e., ⎛S⎞ A = −2 τ ierfc ⎜ (3.522) ⎟ ⎝2 τ ⎠ Thus, the temperature distribution in the liquid phase is ⎡ ⎛X⎞ ⎛ S ⎞⎤ θ ( X ,τ ) = 2 τ ⎢ierfc ⎜ (3.523) ⎟ − ierfc ⎜ ⎟⎥ ⎝2 τ ⎠ ⎝ 2 τ ⎠⎦ ⎣ Substituting eq. (3.523) into eq. (3.519) yields an ordinary differential equation about the location of the solid-liquid interface, i.e., dS ⎛S⎞ = Ste ⋅ erfc ⎜ (3.524) ⎟ dτ ⎝2 τ ⎠ which is subject to the following initial condition: S (τ ) = 0 τ =0 (3.525) A closed-form expression of S (τ ) cannot be obtained from eq. (3.524) because the separation of variables in eq. (3.524) is impossible. Equation (3.524) 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. (3.520) does not satisfy the partial differential equation (3.516) unless A in eq. (3.520) is a constant. However, A in eq. (3.520) is not a constant because it is a function of dimensionless time τ [see eq. (3.522)]. 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 semi-infinite 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 3.4 (see Problem 3.56). A very simplified solution, which will be given in Example 3.5, can provide the first order estimation to the phase change problem.
Chapter 3 Heat Conduction 293
Example 3.5 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 a singlecomponent PCM with a 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 to freeze the entire slab. 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′′ = ρ hs dt Assuming steady-state conduction in the frozen layer, the latent heat released by freezing must be conducted through the frozen layer, i.e., −1 ⎛ s 1⎞ q′′ = (Tm − T∞ ) ⎜ + ⎟ ⎝k h⎠ Combining the above two equations, we have Tm − T∞ ⎛ s 1⎞ dt = ⎜ + ⎟ ds (3.526) ρ hs ⎝k h⎠ The time it takes to freeze the entire slab, tf, can be obtained by integrating eq. (3.526)
Symmetric line Solid Liquid
Tm T0 h , T∞
0
s
L
x
Figure 3.34 Freezing of a slab under convective cooling.
294 Advanced Heat and Mass Transfer
∫
i.e.,
tf
0
Tm − T∞ dt = ρ hs tf =
∫
L
0
⎛ s 1⎞ ⎜ + ⎟ ds ⎝k h⎠
⎛ L2 L ⎞ +⎟ (3.527) ⎜ Tm − T∞ ⎝ 2k h ⎠ which can be nondimensionalized as 1 ⎛1 1 ⎞ (3.528) τf = ⎜+ ⎟ Ste ⎝ 2 Bi ⎠ where τ f = α t f / L2 , Ste = c p (Tm − T∞ ) / hs , Bi = hL / k are dimensionless time, Stefan number and Biot number, respectively. When Bi → ∞ , the problem becomes the boundary condition of the first kind [ T (0, t ) = T∞ ] and eq. (3.528) becomes 1 τf = (3.529) 2Ste
Exact Solution of the Two-Region Problem
ρ hs
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 (3.496) – (3.502) 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 ) (3.530) (3.531) where A and B in eqs. (3.530) and (3.531) are unspecified constants, and erfc is the complementary error function, defined as erfc( z ) = 1 − erf ( z ) (3.532) It should be noted that eqs. (3.530) and (3.531) satisfy eqs. (3.496) – (3.500). The constants A and B can be determined by using boundary condition (3.501), i.e., 1 + Aerf ( S / 2τ 1/ 2 ) = 0
θ 2 ( X ,τ ) = θ i + Berfc[ X / 2( Nατ 1/ 2 )]
θi + Berfc[ S / 2( Nατ 1/ 2 )] = 0 Since A and B are constants, S / 2τ 1/ 2 must also be a constant in order for the above two equations to be satisfied. This constant can be represented by λ , so λ = S / 2τ 1/ 2 (3.533)
Thus, the constants A and B can be determined as
Chapter 3 Heat Conduction 295
A=−
B=−
1 erf (λ )
θi 1 erfc(λ/Nα/ 2 )
The temperature distributions in both phases can be obtained by substituting the above equations for A and B into eqs. (3.530) and (3.531), i.e., erf[ X / 2τ 1/ 2 ] θ1 ( X ,τ ) = 1 − (3.534) erf (λ ) (3.535) ⎣ Substituting eqs. (3.534), (3.535), and (3.533) into eq. (3.502), the following equation is obtained for the constant λ : Nθ λπ e−λ e − λ/Nα + k1/ 2i = (3.536) 1/ 2 erf (λ ) Nα erfc(λ/Nα ) Ste Equation (3.536) 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. (3.534), (3.535), and (3.533), respectively.
2
θ 2 ( X ,τ ) = θ i ⎢1 −
⎡
erf[ X / 2( Nατ )1/ 2 ] ⎤ ⎥ 1 erf (λ/Nα/ 2 ) ⎦
3.5.3 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. In this subsection, the integral approximate method (see Section 3.3.2) will be employed to various melting/solidification problems.
One-Region Problem
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. (3.503) – (3.506). 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
296 Advanced Heat and Mass Transfer
the PCM. Therefore, the thickness of the liquid phase is identical to the thermal penetration depth. Integrating eq. (3.503) with respect to X in the interval of (0, S ) , one obtains ∂θ ( X ,τ ) ∂X where Θ(τ ) = −
X = S (τ )
∂θ ( X ,τ ) ∂X
S (τ )
=
X =0
d Θ(τ ) dτ
(3.537)
∫
0
θ ( X ,τ )dX
(3.538)
Substituting eq. (3.506) into eq. (3.537) yields the integral equation of the one-region problem: 1 dS ∂θ ( X ,τ ) d Θ(τ ) − − = (3.539) ∂X Ste dτ dτ X =0 Assume now that the temperature distribution in the liquid phase is the following second-degree polynomial: ⎛ X −S ⎞ ⎛ X −S⎞ (3.540) ⎟ + A2 ⎜ ⎟ ⎝S⎠ ⎝S⎠ where A0, A1, and A2 are three unknown constants. Equations (3.504) and (3.505) can be used to determine the constants in eq. (3.540). However, eq. (3.506) is not suitable for determining the constant in eq. (3.540) because dS /dτ in eq. (3.506) is unknown. Another appropriate boundary condition at the solid-liquid interface is needed. Differentiating eq. (3.505), one obtains ∂θ ∂θ dθ = dX + dτ = 0 X = S (τ ) ∂X ∂τ i.e., ∂θ dS ∂θ + =0 X = S (τ ) ∂X dτ ∂τ Substituting eqs. (3.503) and (3.506) into the above equation yields
2
θ ( X ,τ ) = A0 + A1 ⎜
∂ 2θ ⎛ ∂θ ⎞ (3.541) Ste⎜ ⎟= 2 ⎝ ∂X ⎠ ∂X Equation (3.541) is an additional boundary condition at the solid-liquid interface; it can be used to determine the coefficients in eq. (3.540). Thus, the constants in eq. (3.540) can be determined by eqs. (3.504), (3.505) and (3.541). After the constants in eq. (3.540) 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 ⎟ (3.542) Ste Ste ⎝S⎠⎜ ⎠ ⎝ ⎠ Substituting eq. (3.542) into the integral equation (3.539) leads to an ordinary differential equation for S (τ ) :
2
Chapter 3 Heat Conduction 297
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 3.35 Comparison between integral and exact solutions of the one-region conduction controlled melting problem.
dS 1 − 1 + 2 Ste + 2 Ste =6 (3.543) dτ 5 + 1 + 2Ste + 2Ste The initial condition of eq. (3.543) is S (τ ) = 0 τ =0 (3.544) The right-hand side of eq. (3.543) is a constant and its solution is 1 S = 2λτ 2 (3.545) where S ⎡ 1 − 1 + 2 Ste + 2Ste ⎤ 2 λ = ⎢3 (3.546) ⎥ ⎣ 5 + 1 + 2Ste + 2 Ste ⎦ Figure 3.35 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.
1
Example 3.6 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
298 Advanced Heat and Mass Transfer
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 3.5 except the heat flux at the surface is a function of time, i.e., ′′ (3.547) q′′(t ) = f (t ) q0 ′ 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θ ∂θ = (3.548) 0 < X < S (τ ), τ > 0 ∂X 2 ∂τ ∂θ = − f (τ ) (3.549) X = 0, τ > 0 ∂X X = S (τ ), τ > 0 θ ( X ,τ ) = 0 (3.550) ∂θ 1 dS − = X = S (τ ), τ > 0 (3.551) ∂X Ste dτ where the nondimensional variables in the above eqs. (3.548) – (3.551) are defined by eq. (3.515). Since eq. (3.549) is a nonhomogeneous boundary condition, the exact solution in Example 3.4 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. (3.548) 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 (τ ) = (3.552) Ste dτ dτ where
Θ(τ ) = Integrating eq. (3.552) with respect to τ in the interval (0, τ ) yields τ S (τ ) + Θ(τ ) = f (τ )dτ (3.554) 0 Ste Assuming that the temperature distribution in the liquid phase is
∫
S (τ )
0
θ ( X ,τ )dX
(3.553)
∫
⎛ X −S ⎞ ⎛ X −S⎞ θ ( X ,τ ) = A0 + A1 ⎜ ⎟ + A2 ⎜ ⎟ ⎝S⎠ ⎝S⎠ where the constants, A0, A1, and A2 are three unspecified constants that can be determined from eqs. (3.549), (3.550) and (3.541). After all unknown constants are determined, the temperature distribution in the liquid phase becomes
2
Chapter 3 Heat Conduction 299
1 ⎢1 − (1 + 4μ ) θ ( X ,τ ) = ⎣ 2 Ste
⎡
1/ 2 ⎤ ⎥ ⎦
⎛ X − S ⎞ 1 ⎢1 − (1 + 4μ ) ⎣ ⎜ ⎟+ Ste ⎝S⎠8
⎡
1/ 2 ⎤ 2 ⎥ ⎦
⎛ X −S⎞ ⎜ ⎟ ⎝S⎠ (3.555) (3.556)
2
where
μ = f (τ ) S (τ ) Ste
Substituting eq. (3.555) into eq. (3.554), one obtains (3.557) 6 For the case of constant heat flux, i.e., f (τ ) = 1, eq. (3.557) can be simplified to S ⎡ SteS + 5 + (1 + 4SteS )1/ 2 ⎤ = 6 Steτ (3.558) ⎢ ⎥ ⎣ ⎦
0
⎢ ⎣
Ste 2 f (τ )
∫
τ
f (τ )dτ =
μ⎡
μ + 5 + (1 + 4μ )1/ 2 ⎤ ⎥ ⎦
Two-Region Problem
The physical model of a melting problem is shown in Fig. 3.36: 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 = (3.559) 2 ∂x α 2 ∂t
∂T2 1 ′′ = − q0 x = 0, 0 < t < tm ∂x k2 T2 ( x, t ) → Ti x → ∞, 0 < t < tm
(3.560) (3.561) (3.562)
T2 ( x, t ) = Ti , 0 < x < ∞, t = 0
300 Advanced Heat and Mass Transfer
T T1(x,t) ′ q0′ T2(x,t) T2(x,tm)
Tm
s(t)
δm
δ(t)
x
Figure 3.36 Melting in a subcooled semi-infinite body under constant heat flux
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⎞ (3.563) T2 ( x, t ) = Ti + ⎜1 − ⎟ 2k2 ⎝ δ ⎠ where δ is the thermal penetration depth. It can be obtained by substituting eq. (3.563) into the integral equation, which can in turn be obtained by integrating eq. (3.559) in the interval of (0, δ ). The result is (3.564) The highest temperature of the semi-infinite body occurs at its surface (x=0) and can be expressed as q′′δ (3.565) Ts (t ) = Ti + 0 2k 2 Melting occurs when the surface temperature reaches the melting point, Tm , and the corresponding thermal penetration depth is 2k (T − T ) δm = 2 m i (3.566) ′′ q0 Then the duration of preheating can be obtained by combining eqs. (3.564) and (3.566), 2k 2 (T − T ) 2 tm = 2 m 2 i (3.567) ′′ 3α 2 q0 The temperature distribution at time tm is
⎛ x⎞ T2 ( x, tm ) = Ti + (Tm − Ti )⎜ 1 − ⎟ ⎝ δm ⎠
2
δ = 6α 2t
(3.568)
Chapter 3 Heat Conduction 301
Governing Equations for the Melting Stage After melting starts, the governing equations in the different phases must be specified separately. The temperature in the liquid phase satisfies ∂ 2T1 1 ∂T1 0 < x < s (t ), t > tm = (3.569) ∂x 2 α1 ∂t ∂T1 ( x, t ) 1 ′′ x = 0, t > tm = − q0 (3.570) k1 ∂x The governing equation and corresponding boundary and initial conditions for the solid phase are ∂ 2T2 1 ∂T2 = (3.571) s (t ) < x < ∞, t > tm 2 ∂x α 2 ∂t T2 ( x, t ) → Ti x → ∞, t > tm (3.572)
⎛ x⎞ (3.573) 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 (3.574) ∂T ∂T ds (3.575) k2 2 − k1 1 = ρ hs 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 ⎪ hs hs hs ⎪
2
X=
x ′′ α1 ρ1hs /q0
S=
s ′′ α1 ρ1hs /q0
Δ=
δ
′′ α1 ρ1hs /q0 α Nα = 2 α1
Δm =
δm ′′ α1 ρ1hs /q0
τ=
′′ (α1 ρ1hs /q0 ) 2
α1t
⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
(3.576)
where Sc is a subcooling parameter that signifies the ratio of the sensible preheat needs to raise the temperature of the body before melting takes place, eqs. (3.569) – (3.575) become ∂ 2θ1 ∂θ1 = 0 < X < S (τ ), τ > τ m (3.577) ∂X 2 ∂τ ∂θ1 ( X ,τ ) X = 0, τ > τ m = −1 (3.578) ∂X ∂ 2θ 2 1 ∂θ 2 = (3.579) S (τ ) < X < ∞, τ > τ m 2 ∂X Nα ∂τ θ 2 ( X ,τ ) → − Sc X → ∞, τ > τ m (3.580)
302 Advanced Heat and Mass Transfer
2 ⎡⎛ ⎤ X⎞ X > 0, τ = τ m θ 2 ( X ,τ ) = Sc ⎢⎜1 − (3.581) ⎟ − 1⎥ ⎢⎝ Δ m ⎠ ⎥ ⎣ ⎦ θ1 ( X ,τ ) = θ 2 ( X ,τ ) = 0 X = S (τ ), τ > τ m (3.582) ∂θ 1 ∂θ 2 dS − 1+ = (3.583) X = S (τ ), τ > τ m ∂X Nα ∂X dτ where Δ m and τ m can be obtained by substituting eq. (3.576) into eqs. (3.566) and (3.567), i.e., Δ m = 2 Nα Sc (3.584) 2 τ m = Nα Sc 2 (3.585) 3 Integral Approximate Solution Figure 3.37 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. (3.579) over the interval ( S , Δ ), then applying the definition of the thermal penetration depth and eq. (3.582), yields the integral equation of the solid phase: ∂θ 1d −2 = (3.586) ( Θ2 + ScΔ ) ∂X X = S Nα dτ
Figure 3.37 Dimensionless temperature distribution for melting in a semi-infinite body under constant heat flux.
Chapter 3 Heat Conduction 303
where Θ2 =
∫
Δ
S
θ 2 dX
(3.587)
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 where the constants B0, B1, and B2 can be obtained from the boundary condition, eq. (3.582) and the definition of the thermal penetration depth ( θ 2 X =Δ = − Sc and ∂θ 2 / ∂X phase is ⎡⎛ Δ − X ⎞ 2 ⎤ (3.588) ⎟ − 1⎥ ⎢⎝ Δ − S ⎠ ⎥ ⎣ ⎦ Substituting eq. (3.588) into eq. (3.586) yields a relationship between the location of the solid-liquid interface and the thermal penetration depth: 6 Nα dS d Δ =2 + (3.589) Δ−S dτ dτ Integrating the differential equation of liquid phase eq. (3.577) over the interval of (0, S) and applying the boundary conditions eqs. (3.578) and (3.582) yields the integral equation of the liquid phase: ∂θ1 d Θ1 +1 = (3.590) ∂X X = S dτ where
X =Δ
= 0). The final form of the temperature distribution in the solid
θ 2 ( X ,τ ) = Sc ⎢⎜
Θ1 =
∫
S
0
θ1dX
(3.591)
Assuming that the temperature profile in the liquid phase is a second-degree polynomial function, ⎛ X −S ⎞ ⎛ X −S⎞ (3.592) ⎟ + A2 ⎜ ⎟ ⎝S⎠ ⎝S⎠ where A0, A1, and A2 need to be determined using appropriate boundary conditions. Equations (3.578) and (3.582) provided two conditions; and the third condition needs to be obtained by differentiating eq. (3.582), i.e., ∂θ ∂θ dθ1 = 1 dX + 1 dτ = 0 X = S (τ ) ∂X ∂τ i.e., ∂θ1 dS ∂θ1 + =0 (3.593) X = S (τ ) ∂X dτ ∂τ Substituting eqs. (3.577) and (3.583) into eq. (3.593) yields 2 1 ∂θ1 ∂θ 2 ∂ 2θ1 ⎛ ∂θ1 ⎞ − = (3.594) ⎜ ⎟ 2 ⎝ ∂X ⎠ Nα ∂X ∂X ∂X
2
θ1 = A0 + A1 ⎜
304 Advanced Heat and Mass Transfer
which is an additional boundary condition at the solid-liquid interface, and can be used to determine the coefficients in eq. (3.592). After the constants in eq. (3.592) are determined, the temperature distribution in the liquid phase becomes
θ1 ( X ,τ ) =
where
S ⎛ X − S ⎞ 1 X 2 − S2 ⎜ ⎟− p 2⎝ S ⎠ 2 S2
2
2
(3.595)
⎛ N ScS 1 ⎞ ⎛ N ScS 1 ⎞ p=⎜ α − ⎟+ ⎜ α − ⎟ +S (3.596) ⎝ Δ−S 2⎠ ⎝ Δ−S 2⎠ Substituting the integral eqs. (3.590) and (3.586) into eq. (3.583) yields dΘ dS d + (3.597) ( Θ2 + ScΔ ) + 1 = 1 dτ dτ dτ Integrating both sides of eq. (3.597) with respect to τ within ( τ m ,τ ), one obtains S + [Θ2 (τ ) − Θ 2 (τ m )] + Sc(Δ − Δ m ) + [Θ1 (τ ) − Θ1 (τ m )] = (τ − τ m ) (3.598) Substituting eqs. (3.587), (3.588), (3.591), and (3.595) into eq. (3.598), the resulting interface equation is 12 S + ( p + 3 + 2Sc) S + Sc(Δ − Δ m ) = 3(τ − τ m ) (3.599) 2 Equation (3.589) is a first-order ordinary differential equation, while eq. (3.599) 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τ (3.600) ⎜ ⎟
⎝ ⎠
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
Chapter 3 Heat Conduction 305
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 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. 3.38). 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 (3.601) + ρ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. (3.601) are T = Ti , t = 0 (3.602)
Figure 3.38 Physical model for ablation.
306 Advanced Heat and Mass Transfer
T = Tm , x = 0 (3.603) ∂T T = Ti , = 0, x→∞ (3.604) ∂x The energy balance at the surface is ∂T , x=0 (3.605) q′′ − ρU a hs = − k ∂x x = 0 Shortly after ablation begins, the process enters steady-state and the ablation velocity becomes a constant. Equation (3.601) can be simplified as U dT d 2T =− a (3.606) 2 dx α dx where the thermal properties are assumed to be independent of temperature. Equation (3.606) can be treated as a first-order differential equation of dT/dx, and its solution is U − ax dT = C1e α (3.607) dx where C1 is an integrating constant. Equation (3.607) can be further integrated to yield U Cα − ax T = − 1 e α + C2 (3.608) Ua The constants C1 and C2 can be determined from eqs. (3.603) – (3.604), and the temperature distribution in the ablating material becomes
T = Ti + (Tm − Ti )e α Substituting eq. (3.609) into eq. (3.605) yields q′′ − ρU a hs = ρ c pU a (Tm − Ti )
− Ua x
(3.609) (3.610)
which can be rearranged to obtain the ablation velocity: q′′ Ua = (3.611) ρ hs (1 + Sc) where Sc is the subcooling parameter defined as c p (Tm − Ti ) Sc = (3.612) hs The fraction of heat removed by ablation over the total heat flux is ′ qa′ ρU a hs = (3.613) q′′ q′′ Substituting the ablation velocity from eq. (3.611) into eq. (3.613), one obtains ′ qa′ 1 (3.614) = q′′ 1 + Sc If the subcooling parameter is Sc = 0.2, 83% of the total heat flux will be removed by ablation.
Chapter 3 Heat Conduction 307
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 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 cylindrical coordinate system 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. 3.39, 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
Figure 3.39 Solidification of an infinite liquid PCM around an ID-cooled cylinder.
308 Advanced Heat and Mass Transfer
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 (3.615) ri < r < s (t ) t > 0 ⎜r ⎟= r ∂r ⎝ ∂r ⎠ α ∂t T (r , t ) = T0 r = ri t > 0 (3.616) T (r , t ) = Tm r = s t > 0 (3.617) ∂T ds k r =s t >0 = ρ hs (3.618) dt ∂r where the subscript 1 for solid phase has been dropped for ease of notation. Defining dimensionless variables, c (T − T ) T −T αt r s R= S= (3.619) θ= m τ = 2 Ste = p m 0 Tm − T0 ri ri ri hs eqs. (3.615) – (3.618) become 1 ∂ ⎛ ∂θ ⎞ ∂θ (3.620) 1 < R < S (τ ) τ > 0 ⎜R ⎟= R ∂R ⎝ ∂R ⎠ ∂τ θ ( R, τ ) = 1 R = 1 τ > 0 (3.621) θ ( R,τ ) = 0 R = S (τ ) τ > 0 (3.622) ∂θ 1 dS − = R=S τ >0 (3.623) ∂R Ste dτ The above eqs. (3.620) – (3.623) are also valid for melting around a hollow cylinder when used with the appropriate dimensionless variables. From elementary heat transfer, it is known 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 ⎞ (3.624) θ =1+ϕ ⎜ ⎟ − (1 + ϕ )⎜ ⎟ ⎝ ln S ⎠ ⎝ ln S ⎠ where ϕ is an unknown variable. Equation (3.624) automatically satisfies eqs. (3.621) and (3.622), and ϕ can be obtained by differentiating eq. (3.622) ∂θ ∂θ dθ = dR + dτ = 0, R = S ∂R ∂τ which can be rearranged to get ∂θ dS ∂θ + =0 (3.625) ∂R dτ ∂τ Substituting eqs. (3.620) and (3.622) into eq. (3.625) yields the following expression: ⎛ ∂θ ⎞ 1 ∂ ⎛ ∂θ ⎞ − Ste⎜ ⎟+ ⎜R ⎟=0 R=S τ >0 ⎝ ∂R ⎠ R ∂R ⎝ ∂R ⎠
2 2
(3.626)
Chapter 3 Heat Conduction 309
ϕ:
Substituting eq. (3.624) into eq. (3.626) gives the following expression for 2 +ϕ =
1 + 2Ste − 1 (3.627) Ste Substituting eqs. (3.624) and (3.627) into eq. (3.623) leads to an ordinary differential equation for the locations of the solid-liquid interface dS 1 + 2Ste − 1 = (3.628) dτ S ln S which is subjected to the following initial condition: S (τ ) = 1, τ =0 (3.629) Integrating equation (3.628) over the time interval (0,τ ) results in the following equation for the location of the solid-liquid interface: 12 1 S ln S − ( S 2 − 1) = 1 + 2 Ste − 1 τ (3.630) 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 the integral approximate method to the solution of solidification around a cylindrical tube convectively cooled from inside.
(
)
Example 3.7 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. (3.615) – (3.618) except that eq. (3.616) needs to be replaced by the following expression: ∂T = hi (T − Ti ) r = ri t > 0 k (3.631) ∂r The dimensionless governing equations of the problem are the same as eqs. (3.620) – (3.623), but eq. (3.621) needs to be replaced by the dimensionless form of eq. (3.631), i.e., ∂θ (3.632) = Bi (θ − 1), R = 1 τ = 0 ∂R where Bi in eq. (3.632) is the Biot number defined as hr Bi = i (3.633) k
310 Advanced Heat and Mass Transfer
and θ is defined as Tm − T (3.634) Tm − Ti It is also assumed that the temperature distribution has a secondorder logarithmic function of the form
θ=
⎛ ln R ⎞ ⎛ ln R ⎞ (3.635) ⎟ − ( A + B)⎜ ⎟ ⎝ ln S ⎠ ⎝ ln S ⎠ The two unknown variables A and B in eq. (3.635) can be determined by substituting eq. (3.635) into eqs. (3.632) and (3.626), with the result B = − Bi (1 − A) ln S (3.636)
2
θ = A + B⎜
1 + 2 ASte − 1 (3.637) Ste Substituting eq. (3.636) into eq. (3.637), an equation for A is obtained as follows: 1 + 2 ASte − 1 2 A − Bi (1 − A) ln S = (3.638) Ste Substituting eq. (3.635) into eq. (3.623), the ordinary differential equation for the location of the solid-liquid interface is obtained: dS 1 + 2 ASte − 1 = (3.639) dτ S ln S which is subject to the initial condition specified by eq. (3.629). The temperature of the inner surface of the tube is θ (1,τ ) = A (3.640) The solid-liquid interface location can be obtained by numerical solution of eq. (3.639) and (3.640). It should be noted that A becomes 1 if the Biot number becomes infinite [see eq. (3.638)]. In that case the temperature of the inner surface becomes 1 and eq. (3.639) reduces to eq. (3.628). 2A + B =
3.5.4 Numerical Simulation
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
Chapter 3 Heat Conduction 311
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. 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, 1990). 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.
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. 3.40. It is assumed that the
312 Advanced Heat and Mass Transfer
densities of the liquid and solid phase are identical ( ρ s = ρ ). The energy equation can be written as ∂h ∂ ⎛ ∂T ⎞ ρ = ⎜k (3.641) ⎟ ∂t ∂x ⎝ ∂x ⎠ The enthalpy, h, is a function of temperature, T: ⎧ T < Tm ⎪ c ps (T − Tm ) ⎨ (3.642) h(T ) = ⎪ ⎪ c (T − T ) + h m s T > Tm ⎪p ⎩ T T0 T(x,t) Tm Ti
x=0
s
x=L
x
Figure 3.40 Conduction-controlled melting in a finite slab.
h h (T )
hs
0 Tm
T
Figure 3.41 Enthalpy of the PCM as a function of temperature.
Chapter 3 Heat Conduction 313
The variation of enthalpy h with temperature T can be plotted as shown in Fig. 3.41. At the melting point of the PCM, the enthalpies of solid and liquid phases at the melting point are 0 and hs , 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 (3.643) k (T ) = ⎨ ⎪k T > Tm ⎪ ⎩ It can be verified that the energy equation for liquid and solid phases can be obtained simply by substituting eqs. (3.642) and (3.643) into eq. (3.641). 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. (3.642) yields ⎧ h≤0 ⎪Tm + h / c ps
⎨ T = ⎪Tm ⎪ ⎪ ⎪m ⎩ ⎪
0 < h < hs h ≥ hs
(3.644)
T + ( h − hs ) / c p
The initial and boundary conditions of the melting problem are T ( x, t ) = Ti < Tm t =0 (3.645) T ( x, t ) = T0 > Tm x=0 (3.646) ∂T ( x, t ) =0 x=L (3.647) ∂x The discretization of space and time is shown in Fig. 3.42, 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 (3.641) is discretized by using an explicit scheme which uses forward difference in time and central difference in space: 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 (3.648) 2 Δt ( Δx )
Figure 3.42 Discretization of space and time for the solution of the two-region melting problem.
314 Advanced Heat and Mass Transfer
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 k j−1 =
2
(3.649) (3.650)
2k j −1k j k j −1 + k j
Equation (3.648) 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⎟ ⎟ ⎟ ⎠
(3.651) 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. (3.644), i.e., ⎧ Tm + h n +1 / c ps h n +1 ≤ 0 (solid) j j ⎪ n +1 n +1 Tj = ⎨ Tm 0 < h j < hs (interface) j = 2,..., N − 1 (3.652) ⎪T + (h n +1 − h ) / c (liquid) h n +1 ≥ hs j s p j ⎩m The initial and boundary conditions are discretized as T j0 = Ti j = 1, 2, ..., N (3.653) T1n +1 = T0
n +1 N n +1 N −1
(3.654)
T =T (3.655) th The location of the solid-liquid interface at the (n+1) 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 < hs , the solid-liquid interface will be n Δx hM+1 + Δx s n +1 = ( M − 1)Δx − (3.656) 2 hs The solution procedure can be summarized as follows:
1. Determine the initial enthalpy at every node h0 by using eq. (3.642). j 2. Calculate the enthalpy after the first time step at nodes j = 2, ..., N − 1 by using eq. (3.651). 3. Determine the temperature after the first time step at node j = 1, ..., N by using eqs. (3.654), (3.652), and (3.655).
Chapter 3 Heat Conduction 315
4. Find a control volume in which the enthalpy falls between 0 and hs and determine the location of the solid-liquid interface by using eq. (3.656). 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 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 ,α )Δt 1 ≤ (3.657) 2 2 ( Δx ) where α s and α are thermal diffusivity of solid and liquid, respectively. Equation (3.657) is valid for one-dimensional problems only. In addition 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.97)] ∂ ( ρ h) ∂ ( ρ uh) ∂ ( ρ vh) ∂ ( ρ wh) ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ + + + = ⎜k ⎟+ ⎜k ⎟+ ⎜k ⎟ ∂t ∂x ∂y ∂z ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ (3.658)
316 Advanced Heat and Mass Transfer
with the state equation dh = c(T ) (3.659) ∂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 by eq. (3.644). Substituting eq. (3.644) into eq. (3.658), the energy equation can be transformed into the following form (Cao and Faghri, 1989): ∂ ( ρ h) ∂ ( ρ uh) ∂ ( ρ vh) ∂ ( ρ wh) + + + ∂t ∂x ∂y ∂z (3.660) 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 ⎪ Γ ( h ) = ⎨0 (3.661) 0 < h < hs ⎪k / c h > hs p ⎩
⎧0 h≤0 ⎪ 0 < h < hs (3.662) S ( h ) = ⎨0 ⎪ hk / c h > h p s ⎩ 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.
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
Chapter 3 Heat Conduction 317
Figure 3.43 Dependence of specific heat to temperature for equivalent heat capacity model.
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 T < Tm − ΔT ⎧c ps ⎪ c ps + c p ⎪h c p (T ) = ⎨ s + Tm − ΔT < T < Tm + ΔT (3.663) 2 ⎪ 2 ΔT T > Tm + ΔT ⎪c p ⎩ which assumes that the temperature of the PCM is changed from Tm − ΔT to Tm + ΔT when latent heat, hs , 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. 3.43. 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 (3.664) = ⎜k ⎟+ ⎜k ⎟ + ⎜k ⎟ ∂t ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ where the thermal conductivity in eq. (3.664) 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
318 Advanced Heat and Mass Transfer
T < Tm − ΔT ⎧ks ⎪ k − ks ⎪ k (T ) = ⎨ k s + (T − Tm + ΔT ) Tm − ΔT < T < Tm + ΔT 2ΔT ⎪ T > Tm + ΔT ⎪k ⎩
(3.665)
The advantage of the equivalent heat capacity model is its simplicity. Equation (3.664) 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.
Temperature-Transforming Model
The temperature-transforming model proposed by Cao and Faghri (1990) 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 ⎞ (3.666) ρ = ⎜k ⎟+ ⎜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. 3.44. This relationship can be analytically expressed as (see Problem 3.73) ⎧ c ps (T − Tm ) + c ps ΔT T < Tm − ΔT ⎪ c ps + c p h⎞ h ⎪⎛ c ps + c p h(T ) = ⎨⎜ + s ⎟ (T − Tm ) + ΔT + s Tm − ΔT < T < Tm + ΔT 2 2ΔT ⎠ 2 2 ⎪⎝ ⎪ c p (T − Tm ) + c ps ΔT + hs T > Tm + ΔT ⎩ (3.667) By defining specific heat in the mushy zone as c ps + c p (3.668) cm = 2 eq. (3.667) becomes
Chapter 3 Heat Conduction 319
Figure 3.44 Dependence of enthalpy on temperature for phase change occurring over a range of temperature.
⎧c ps (T − Tm ) + c ps ΔT ⎪ h⎞ h ⎪⎛ h(T ) = ⎨⎜ cm + s ⎟ (T − Tm ) + cm ΔT + s 2ΔT ⎠ 2 ⎪⎝ ⎪c p (T − Tm ) + c ps ΔT + hs ⎩ Equation (3.669) can be rewritten as h(T ) = c p (T )(T − Tm ) + b(T )
T < Tm − ΔT Tm − ΔT < T < Tm + ΔT (3.669) T > Tm + ΔT (3.670)
where c p (T ) and b(T ) can be determined from eq. (3.669): T < Tm − ΔT ⎧c ps ⎪ h ⎪ (3.671) c p (T ) = ⎨cm + s Tm − ΔT < T < Tm + ΔT 2ΔT ⎪ T > Tm + ΔT ⎪c p ⎩ T < Tm − ΔT ⎧c ps ΔT ⎪ h ⎪ (3.672) b(T ) = ⎨cm ΔT + s Tm − ΔT < T < Tm + ΔT 2 ⎪ T > Tm + ΔT ⎪c ps ΔT + hs ⎩ Substituting eq. (3.670) into eq. (3.666) yields ∂ (c pT ) ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂b (3.673) ρ = ⎜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. (3.665). The energy equation has been transformed into a
320 Advanced Heat and Mass Transfer
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. (3.673) reduces to a linear equation. The temperature-transforming model and the equivalent heat capacity model differ in significant ways. Comparing eq. (3.673) and eq. (3.664), one can conclude that the equivalent heat capacity model is a special case of eq. (3.673), with ∂b/∂t = 0 and cp independent of the spatial variables x, y, z and time. This is 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 (3.673) 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 3.45 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. (3.495).
Ste = 0.2 Ks = 1 Cs = 1
S
τ1/2
Figure 3.45 Comparison of interface location for the Stefan problem obtained by exact and numerical solutions using temperature transforming model.
Chapter 3 Heat Conduction 321
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).
3.6 Microscale Heat Conduction
Microscale heat transfer has drawn the attention of many researchers due to its importance in nanotechnology, information technology, and biotechnology. 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.
3.6.1 Extensions of Classic Model
Hyperbolic Model
The classical heat conduction theory based on Fourier’s law assumes that thermal disturbance propagates with an infinite speed. As heat conduction is accomplished by successive collision of the energy carriers (phonons or electrons), the propagation of thermal disturbance is always at a finite speed. This is particularly important for those processes involving extremely short times, cryogenic temperatures, or high heat fluxes. To account for the finite propagation speed of thermal wave, the Cattaneo-Vernotte thermal wave model can be used ∂q′′ = − k ∇T (3.674) q′′ + τ ∂t where τ is the relaxation time that can be interpreted as the time scale at which intrinsic length scale of thermal diffusion ( α t ) is equal to the intrinsic length
322 Advanced Heat and Mass Transfer
scale of thermal wave ( ct ) (Tzou and Li, 1993; Tzou, 1997), where c is the thermal propagation speed. Thus the relaxation time is τ = α / c2 (3.675) Substituting eq. (3.674) into eq. (3.1), the following energy equation is obtained 1 ∂T τ ∂ 2T 1⎛ ∂q′′′ ⎞ (3.676) + = ∇ 2T + ⎜ q′′′ + τ ⎟ 2 k⎝ ∂t ⎠ α ∂t α ∂t Since the second order derivative of temperature with respect to time appears in eq. (3.676), it is referred to as the hyperbolic heat conduction model. Mathematically, eq. (3.676) is a thermal wave equation and the thermal diffusivity appears as a dumping effect of thermal propagation. With appropriate relaxation times, eq. (3.676) can be used to describe temperatures of different energy carriers – such as phonon and electron temperatures – as will be discussed later.
Dual-Phase Lag Model
Equation (3.674) can be viewed as the first order approximation of the following equation, q′′(r, t + τ ) = − k ∇T (r, t ) (3.677) which indicates that there is a delay between the heat flux vector and the temperature gradient. For the same point in the conduction medium, the temperature gradient is established at time t, but the heat flux vector will be established at a later time t + τ , i.e., the relaxation, τ , can be interpreted as the time delay from the onset of the temperature gradient to the heat flux vector. While the thermal wave model assumes that the temperature gradient always precedes the heat flux, Tzou (1997) proposed a dual-phase lag model that allows either the temperature gradient (cause) to precede heat flux vector (effect) or the heat flux vector (cause) to precede the temperature gradient (effect), i.e., q′′(r, t + τ q ) = − k ∇T (r, t + τ t ) (3.678) where τ q is the phase lag for the heat flux vector, while τ T is the phase lag for the temperature gradient. If τ q > τ T , the local heat flux vector is the result of the temperature gradient at the same location but an early time. On the other hand, if τ q < τ T , the temperature gradient is the result of the heat flux at an early time. The first order approximation of eq. (3.678) is: ∂q′′ ∂ ⎡ ⎤ (3.679) = − k ⎢∇T + τ T (∇T ) ⎥ q′′ + τ q ∂t ∂t ⎣ ⎦ Substituting eq. (3.679) into eq. (3.1), the energy equation based on the dualphase lag model is 1 ∂T τ q ∂ 2T 1⎛ ∂ ∂q′′′ ⎞ + = ∇ 2T + τ T (∇ 2T ) + ⎜ q′′′ + τ q ⎟ (3.680) 2 ∂t ∂t ⎠ k⎝ α ∂t α ∂t
Chapter 3 Heat Conduction 323
which reduces to the parabolic conduction equation (3.10) if both τ q and τ T are zero. In the absence of phase lag for temperature gradient ( τ T = 0 ), eq. (3.680) is reduced to the hyperbolic conduction model, eq. (3.676).
3.6.2 Two-Step Model
During laser-metal interaction, the laser energy is first deposited into electrons on the metal surface, where two competing processes occur (Hohlfeld, 2000). One is ballistic motion of the excited electrons into deeper parts of the metal with velocity close to the Fermi velocity (~106 m/s). Another is collision between the excited electrons and electrons around the Fermi level – an electron temperature is defined upon establishment of equilibrium among hot electrons. These hot electrons are then diffused into deeper parts of the electron gas at a speed (<104 m/s) much lower than that of the ballistic motion. Meanwhile, the hot electrons are cooled by transferring their energy to the lattice through electron-phonon coupling. The nonequilibrium between electrons and lattice has been observed experimentally (Eesley, 1986; Elsayed-Ali et al., 1987) and can be described by the two-temperature model, which was originally proposed by Anisimov et al. (Anisimov et al., 1974) and rigorously derived by Qiu and Tien (1993) from the Boltzmann transport equation.
Parabolic Two-Step Model
Assuming heat conduction in the electron can be described by Fourier’s law and neglecting heat conduction in the lattice, the energy equations of the free electrons and lattices (phonons) are ∂T Ce e = ∇ ⋅ ( ke ∇Te ) − G (Te − Tl ) + q′′′ (3.681) ∂t ∂T Cl l = G (Te − Tl ) (3.682) ∂t where the volumetric heat capacity of lattice is Cl = ρcp, and the volumetric heat capacity of electrons is π 2 ne k B (3.683) Ce = Te = BeTe 2μ F where ne is the number density of electrons, kB is the Boltzmann constant, and μF is Fermi energy. Equation (3.683) indicates that the volumetric heat capacity of the electron is proportional to the electron temperature. It should be noted that the volumetric heat capacity of electrons is much less than that of the lattice even at very high electron temperature. At nonequilibrium condition, thermal conductivity of the electrons depends on the temperatures of both electrons and lattice, i.e.,
324 Advanced Heat and Mass Transfer
⎛T ⎞ ke = keq ⎜ e ⎟ ⎝ Tl ⎠
(3.684)
where keq(T) is the thermal conductivity of the electron when the electrons and lattice are in thermal equilibrium. The electron-lattice coupling factor, G, is to account for the rate of energy exchange between electrons and phonons and it can be estimated by 22 9 ne k BTD vF (3.685) G= 16 Λ (Tl )Tl μ F where TD is Debye temperature, vF is Fermi velocity, and Λ is the electron mean free path. Neglecting conduction in the lattice is justified by the fact that the thermal conductivity of the lattice is two orders of magnitude smaller than that of the free electrons (Klemens and Williams, 1986). The heat conduction model represented by eqs. (3.681) and (3.682) is referred to as a parabolic two-step model because Fourier’s law was used to describe heat transfer in the electron gas.
Example 3.8 Assuming all properties of electrons and lattice are independent from temperatures, obtain a single energy equation for lattice temperature by combining eqs. (3.681) and (3.682). Also compare the result with the dual-phase lag model. Solution: Solving for Te from eq. (3.682) yields C ∂T Te = Tl + l l (3.686) G ∂t Substituting eq. (3.686) into eq. (3.681), we have C ∂T ⎞ C ∂T ⎞ ∂⎛ ⎛ ⎛ C ∂T ⎞ Ce ⎜ Tl + l l ⎟ = ke ∇ 2 ⎜ Tl + l l ⎟ − G ⎜ l l ⎟ + q′′′ (3.687) ∂t ⎝ G ∂t ⎠ G ∂t ⎠ ⎝ ⎝ G ∂t ⎠ which can be rearranged as Ce + Cl ∂T Ce Cl ∂ 2T C∂ q′′′ + = ∇ 2Tl + l (∇ 2T ) + (3.688) 2 ke G ∂t ke ∂t Gke ∂t where the subscript “l” for lattice has been dropped. Comparing eq. (3.688) with the energy equation for the dual-phase lag model, eq. (3.680), it is apparent that they have almost identical form except the partial derivative of heat source with respective to time is not present in eq. (3.688). The thermophysical properties in the dual-phase lag model are related to the properties appearing in the two-temperature model by ke C Ce Cl k = ke , α = , τT = l , τ q = (3.689) Ce + Cl G G (Ce + Cl ) The ratio of two phase-lag times is
Chapter 3 Heat Conduction 325
C τ T Ce + Cl = =1+ l τq Ce Ce which indicates that τ T is always greater than τ q .
(3.690)
Hyperbolic Two-Step Model
If we consider the hyperbolic effect on the conduction in the electron gas, the energy equation for the electron gas is ∂T Ce e = −∇ ⋅ q′′ − G (Te − Tl ) + q′′′ (3.691) ∂t where ∂q′′ = − ke ∇Te (3.692) q′′ + τ e ∂t while the energy equation for the lattice is still eq. (3.682). Equations (3.691) and (3.692) can be combined to yield ∂T ∂ 2T Ce e + Ceτ e 2e = ∇ ⋅ (ke∇Te ) ∂t ∂t (3.693) ∂ ∂q′′′ −G (Te − Tl ) − τ e ⎡G (Te − Tl ) ⎤ + q′′′ + τ e ⎦ ∂t ⎣ ∂t The conduction model represented by eqs. (3.693) and (3.682) is referred to as a hyperbolic two-step model. Qiu and Tien (1993) simulated picosecond lasermetal interaction and concluded that the parabolic two-step model predicted the general temperature response but it failed to predict the finite speed of energy propagation. Therefore, the hyperbolic two-step model can provide better accuracy for ultrafast laser-metal interaction.
Dual-Parabolic Two-Step Model
The contribution of heat conduction in phonons was neglected in the above two models. If it is assumed that the heat conduction in the phonons can be modeled using the classical Fourier’s law, the energy equations of the lattices (phonons) are ∂T Cl l = ∇ ⋅ ( kl ∇Tl ) + G (Te − Tl ) (3.694) ∂t The bulk thermal conductivity of metal measured at equilibrium, keq, is the sum of electron thermal conductivity, ke, and the lattice thermal conductivity, kl. Since the mechanism for heat conduction in metal is diffusion of free electron, ke is usually dominant. For gold, ke is 99% of keq, while kl only contributes to 1% of keq (Klemens and Williams, 1986). Since both eqs. (3.681) and (3.694) are parabolic, this model is referred to as dual-parabolic two-step model. For the
326 Advanced Heat and Mass Transfer
case that phonon temperature gradient is significant, inclusion of conduction in the phonons is essential.
Dual-Hyperbolic Two-Step Model
For the case that heat conduction in both electrons and phonons needs to be considered using hyperbolic model, the energy equation for the lattice becomes ∂T Cl l = −∇ ⋅ q′′ + G (Te − Tl ) (3.695) l ∂t where ∂q′′ q′′ + τ l l = − kl ∇Tl (3.696) l ∂t Combining (3.691) and (3.692) to eliminate q′′ yield l ∂Tl ∂ 2Tl ∂ + Clτ l 2 = ∇ ⋅ (kl ∇Tl ) + G (Te − Tl ) + τ l ⎡G (Te − Tl ) ⎤ (3.697) Cl ⎦ ∂t ∂t ∂t ⎣ Equation (3.697) together with eq. (3.693) become governing equations for the dual-hyperbolic two-step model. Chen and Beraun (2001) applied the dualhyperbolic two-step model to simulate ultrashort laser pulse interactions with metal film. They found that the electron temperatures obtained from the dualhyperbolic model and the hyperbolic model are very close. However, the lattice temperatures obtained from the two models differ significantly.
3.6.3 Microscale Phase Change
Short-pulsed 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 electronlattice 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 controlled heat transfer at the solidliquid interface as suggested by eq. (3.479). 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
Chapter 3 Heat Conduction 327
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. Short-pulsed laser melting of a 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 (see Fig. 3.46) is solved by Kuo and Qiu (1996). The problem can be approximated to be onedimensional because the radius of the laser beam is significantly larger than the film thickness of the gold. The dual-parabolic two-step model can be used to describe this process and the energy equations of the free electrons and the lattice are ∂T ∂T ⎤ ∂⎡ ′′′ (3.698) Ce (Te ) e = ⎢ ke (Te , Tl ) e ⎥ − G (Te − Tl ) + qlaser ∂t ∂x ⎣ ∂x ⎦ ∂T ∂ ⎛ ∂T ⎞ (3.699) Cl l = ⎜ kl l ⎟ + G (Te − Tl ) ∂t ∂x ⎝ ∂x ⎠ The source terms in eq. (3.698) can be described by the following equation: 2 ⎡x ⎛t⎞⎤ 1− R ′′′ (3.700) qlaser = 0.94 J exp ⎢ − − 2.77 ⎜ ⎟ ⎥ ⎜⎟ t pδ ⎢δ ⎝ tp ⎠ ⎥ ⎣ ⎦ 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).
Figure 3.46 Laser melting of thin film.
328 Advanced Heat and Mass Transfer
For the conventional melting process, the velocity of the solid-liquid interface is obtained by energy balance as specified by eq. (3.482). 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 T − T ⎞⎤ ⎛ h ⎞⎡ ds = Vs exp ⎜ − s ,a ⎟ ⎢1 − exp ⎜ − s ,a m l , I ⎟ ⎥ (3.701) ⎜ ⎟ dt ⎝ kbTm ⎠ ⎢ ⎝ kbTm Tl , I ⎠ ⎥ ⎣ ⎦ where Vs is the maximum interface velocity that can be approximated as the speed of sound in the liquid phase. hs , a is the latent heat of fusion per atom, kb is the Boltzmann constant, and Tl , I 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 Te ( x, −2t p ) = Tl ( x, −2t p ) = Ti (3.702) The boundary conditions of the problem can be specified by assuming that the heat loss from the film surface can be neglected, i.e., ∂Te ∂T ∂T ∂T =e =l =l =0 (3.703) ∂x x = 0 ∂x x = L ∂x x = 0 ∂x x = L Equations (3.698) and (3.699) 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, Te = Tl = T , one can add eqs. (3.698) and (3.699) to obtain the following one-step model: ∂T ∂ ⎛ ∂T ⎞ ′′ (3.704) = ⎜k C ⎟ + ql′aser ∂t ∂x ⎝ ∂x ⎠ where C = Ce + Cl and k = ke + kl are the heat capacity and thermal conductivity of the thin film, respectively. In arriving at eq. (3.704), 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. (3.704) are T ( x, −2t p ) = Ti (3.705) ∂T ∂T = =0 (3.706) ∂x x = 0 ∂x x = L The one-step model represented by eqs. (3.704) – (3.706) is presented here in comparison with the two-step model.
Chapter 3 Heat Conduction 329
Figure 3.47 Comparison between two-step and one-step models ( t p and
= 20 ps , L = 100 nm ,
R = 0.6 ; Kuo and Qiu, 1996).
Numerical solutions for both two-step and one-step models were obtained by Kuo and Qiu (1996). Figure 3.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
330 Advanced Heat and Mass Transfer
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. Anisimov, S. I., Kapeliovich, B. L., Perel’man, T. L., 1974, “Electron Emission from Metal Surface Exposed to Ultrashort Laser Pulses,” Sov. Phys. JETP, Vol. 39, pp. 375-377. Arpaci, V.S., 1966, Conduction Heat Transfer, Addison-Wesley Pub. Co., Reading, MA. 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. 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, 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., 1990, “A Numerical Analysis of Phase Change Problem including Natural Convection,” ASME Journal of Heat Transfer, Vol. 112, pp. 812-815. Chen, J. K., and Beraun, J. E., 2001, “Numerical Study of Ultrashort Laser Pulse Interactions with Metal Films,” Numerical Heat Transfer A: Applications, Vol. 40, pp. 1-20. 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.
Chapter 3 Heat Conduction 331
Crank, J., 1984, Free and Moving Boundary Problems, Clarendon Press, Oxford. Eckert, E.R.G., and Drake, R.M., 1987, Analysis of Heat and Mass Transfer, Hemisphere, Washington, DC. Eesley, G. L., 1986, “Generation of Non-equilibrium Electron and Lattice Temperatures in Copper by Picosecond Laser Pulses,” Physical Review B, Vol. 33, pp. 2144-2155. 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. Elsayed-Ali, H. E., Norris, T. B., Pessot, M. A., and Mourou, G. A., 1987, “Time-Resolved Observation of Electron-Phonon Relaxation in Copper,” Physical Review Letters, Vol. 58, pp. 1212-1215. Fletcher, L.S., 1988, “Recent Developments in Contact Heat Transfer,” ASME J. Heat Transfer, Vol. 110, pp. 1059-1070. Holman, J.P., 2002, Heat Transfer, 9th ed., McGraw-Hill, New York, NY. Hohlfeld, J., Wellershoff, S. S., Gudde, J., Conrad, U., Jahnke, V., and Matthias, E., 2000, “Electron and Lattice Dynamics Following Optical Excitation of Metals,” Chemical Physics, Vol. 251, pp. 237-258. 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. Klemens, P. G., and Williams, R. K., 1986, “Thermal Conductivity of Metals and Alloys,” Int. Metals Review, 31, pp. 197-215. 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. Murthy, J.Y., Minkowycz, W.J., Sparrow, E.M., and Mathur, S.R., “Survey of Numerical Methods,” Handbook of Numerical Heat Transfer, 2nd ed., eds., Minkowycz, W.J., Sprrow, E.M., and Murthy, J.Y., pp. 3-51, Wiley, New York. Ozisik, M.N., 1993, Heat Conduction, 2nd ed., Wiley-Interscience, New York. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemishpere, Washington, DC. 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.
332 Advanced Heat and Mass Transfer
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. Rohsenow, W.W, Hartnett, J.P., and Cho, Y.I., 1998, Handbook of Heat Transfer, 3rd ed., McGraw-Hill, New York. 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. Tao, W.Q., 2001, Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese). Tzou, D.Y., and Li, J., 1993, “Thermal Wave Emanating from a Fast-Moving Heat Source with a Finite Dimension,” ASME Journal of Heat Transfer, Vol. 115, pp. 526-532. 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. 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., 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.
Problems
3.1. The vector form of the energy equation for heat conduction in an anisotropic medium is eq. (3.4). Show that the energy equation becomes (3.5) in a Cartesian coordinate system. Simplify the energy equation (3.4) for heat conduction in an isotropic material for a cylindrical coordinate system. What is the energy equation for heat conduction of an isotropic material in a spherical coordinate system. A double-layer glass window is consistent of two layers of glass with thickness of 5 mm each separated by an air gap. The inner surface temperature of the inner glass is 15 ºC while the outer surface of the outer glass is at -15 ºC. Assuming there is no natural convection in the air gap,
3.2. 3.3. 3.4.
Chapter 3 Heat Conduction 333
what is the heat loss per unit area by this window if the thickness of air gap is 5 mm? If there is only one layer of glass, what is the heat loss? 3.5. Derive the energy equation, (3.47), for generalized 1-D heat conduction with internal heat generation by performing an energy balance for the differential control volume A(s)ds as shown in Fig. 3.2. If the internal heat source in the finite slab is q′′′ = Ax(1 − x) , how would you modify eqs. (3.54) and (3.55) Consider a hollow cylinder with inner and outer diameters of ri and ro, respectively. If the inner and outer surface temperatures are T1 and T2, respectively, and the intensity of internal heat source is uniformly q′′′, what is the temperature distribution in the cylinder? Repeat the above problem for a hollow sphere. For a straight fin with an adiabatic tip, show that the temperature distribution, total heat transfer, and fin efficiency can be expressed as eqs. (3.76) – (3.78) by solving eq. (3.70) with appropriate boundary conditions.
3.6. 3.7.
3.8. 3.9.
3.10. For a straight fin with a convection tip, show that the temperature distribution, total heat transfer, and fin efficiency can be expressed as eqs. (3.79) – (3.81) by solving eq. (3.70) with appropriate boundary conditions. 3.11. The diameter of the straight fin shown in Fig. 3.8 is equal to D at the base and gradually decreases following a parabolic curve to zero at the tip. If the local diameter at x is expressed by D(1 − x / L) 2 , obtain the temperature distribution, heat transfer rate, and fin efficiency. 3.12. A cylindrical biological tissue with a radius of 10 mm and an initial temperature of 37 ºC is exposed to air with a temperature of 25ºC. The convective heat transfer coefficient between the tissue and fluid is 2 W/m2K. The properties of the biological tissue are ρ = 1000 kg/m3, k = 0.628 W/m-K, c = 4187 J/kg-K. The thermal physical properties and temperature of blood are: ρb = 1.06×103 kg/m3, cpb = 3860 J/(kg•K), ωb = 1.87×10-3 ′′′ m3/(m3tissue•s), Tb = 37 oC. The metabolic heat generation is qm = 3 3 1.19×10 W/m . Find the heat loss per unit length from at the surface of the tissue and the tissue surface temperature. 3.13. If the biological tissue in the above problem is exposed to water at 25ºC and the convective heat transfer coefficient between the tissue and the water is 200 W/m2-K, what is the heat loss per unit length of the tissue and the tissue surface temperature? 3.14. If the left and bottom of the two-dimensional domain shown in Fig. 3.10 are adiabatic, the right side is kept at T1, and the top boundary has variable temperature f(x), what is the temperature distribution in the twodimensional domain?
334 Advanced Heat and Mass Transfer
3.15. The left, right and bottom of the two-dimensional domain shown in Fig. 3.10 are in contact with fluid at temperature T∞ and the convective heat transfer coefficient at these three boundaries are h. The top boundary has variable temperature f(x). Find the temperature distribution in the twodimensional domain. 3.16. Show that if ϕ ( x) satisfies eqs. (3.148) – (3.150) and ψ ( x, y ) satisfies eqs. (3.151) – (3.155), ϑ ( x, y ) = ϕ ( x) + ψ ( x, y ) will satisfy eqs. (3.142) – (3.146), i.e., solution of an nonhomogeneous problem can be obtained by superposing solutions of two simpler problems. 3.17. Show that The solution of heat conduction with internal heat generation can also be obtained by solving a 1-D heat conduction problem in the ydirection with internal heat source and a 2-D heat conduction problem without internal heat generation, i.e., ϑ ( x, y ) = ϕ ( y ) + ψ ( x, y ) . Formulate the equations for ϕ ( y ) and ψ ( x, y ) and obtain the solution of eqs. (3.142) – (3.146) using the equations that you formulated. 3.18. A rectangular domain with a dimension of 2 L × 2 H and uniform heat generation is in contact with fluid at T∞ and the convective heat transfer coefficient is h (see Fig. P3.1). Obtain the temperature distribution in the rectangular domain. y H h , T∞ –L 0 h, T∞ q′′′ h, T∞ L x
–H
h , T∞
Figure P3.1
3.19. A horizontal steam pipe is buried underground and the distance from the ground surface to the center of the steam pipe is 0.5 m (see Fig. 3.12). The length and outer diameter of the steam pipe are 300 m and 15 cm, respectively. If the surface temperature of the steam pipe is 95 ºC and the ground surface temperature is 10 ºC, what is the heat loss from the steam pipe per unit time? The thermal conductivity of the soil en as 0.5 W/m-K. 3.20. The steam pipe in the previous problem is made of carbon steel with inner diameter of 0.4 m. If the steam temperature is 100 ºC and the convective
Chapter 3 Heat Conduction 335
heat transfer coefficient between the steam and the inner surface of the pipe is 20 W/m2-K, what is the heat loss from the steam pipe? 3.21. If the ground surface in Fig. 3.12 is insulated, i.e., (∂T / ∂y ) y = 0 = 0 , the heat loss from the buried cylinder can be derived by placing an imagery heat source at y = − a . If the heat loss will be expressed as q = kS (T − T∞ ) , where T∞ is the ground temperature far away from the buried cylinder, find the shape factor. 3.22. An arbitrarily shaped object with volume V, surface area As, and a uniform initial temperature of Ti is shown in Fig. 3.15. At time t = 0, the arbitrarily shaped object is cooled by radiative heat exchange with surrounding at a temperature of T∞. Assuming lumped capacitance method is valid and the surface emissivity of the object is ε, obtain the transient temperature of the object. 3.23. Steel spheres 14 mm in diameter are annealed by heating to 1200 K and then slowly cooling to 400 K in an air environment with ambient temperature of T∞ = 300 K and h = 10 W/m2-K. Assuming the properties of steel to be k = 50 W/m-K, ρ = 8000 kg/m3, and cp = 650 J/kg-K, estimate the time required for the cooling process. 3.24. Show that the separation constant in eq. (3.209) cannot be a complex number. 3.25. A finite slab with thickness of 2L is initially at a uniform temperature Ti. For times t > 0 the boundary at x = − L and x = L are kept at T∞. Give PDE and the corresponding boundary and initial conditions of the problem. Nondimensionlize the equations by introducing appropriate dimensionless variables. Solve the dimensionless equations to get the solution of dimensionless temperature. 3.26. A long cylinder with radius of ro has a uniform initial temperature of Ti. From time t = 0, the surface temperature of the cylinder is reduced to and kept at T∞ ( T∞ < Ti ). It is assumed that there is no internal heat generation and thermophysical properties are constants. What is the transient temperature distribution in the cylinder? 3.27. A long hollow cylinder with inner radius of ri and outer radius of ro has a uniform initial temperature of Ti. From time t = 0, both inner and outer surfaces are in contact with fluid at a temperature of T∞ ( T∞ < Ti ). The convective heat transfer coefficient at both inner and outer surface is h. Find the transient temperature distribution in the cylinder. 3.28. Solve eqs. (3.262) – (3.265) using the method of separation of variables and show that the solution of eqs. (3.257) – (3.260) is eq. (3.266).
336 Advanced Heat and Mass Transfer
3.29. A sphere with radius of ro has a uniform initial temperature of Ti. From time t = 0, the surface temperature of the sphere is reduced to and kept at T∞ ( T∞ < Ti ). Assuming that there is no internal heat generation and thermophysical properties are constants, find the transient temperature distribution in the sphere. 3.30. A spherical shell with inner radius of ri and outer radius of ro has a uniform initial temperature of Ti. From time t = 0, both inner and outer surfaces temperature are reduced to and kept at T∞ ( T∞ < Ti ). There is no internal heat generation and thermophysical properties are constants. Find the transient temperature distribution in the sphere. 3.31. Solve eqs. (3.282) – (3.285) using the method of separation of variables and show that its solution is eq. (3.286). 3.32. A finite slab ( 0 ≤ x ≤ L ) has an initial temperature of Ti. At time t=0, the temperatures at x = 0 and x = L are suddenly changed to T1 and T2 respectively ( T1 > T2 > Ti ). Obtain the temperature distribution using the method of partial solution. 3.33. A finite slab ( 0 ≤ x ≤ L ) has an initial temperature of Ti and a uniform internal heat source of q′′′ . After time t = 0, the left side ( x = 0 ) is wellinsulated while the temperatures on the right side (x = L) is kept at Ti. Obtain the temperature distribution using the method of partial solution. 3.34. Solve the previous problem using the method of variation of parameter. 3.35. A finite slab ( 0 ≤ x ≤ L ) has an initial temperature of T ( x,0) = f ( x) and a non-uniform, time dependent, internal heat source of q′′′( x, t ) . After time t = 0, the temperatures on left side and right sides are T (0, t ) = g 0 (t ) and T (0, t ) = g L (t ) , respectively. Obtain the temperature distribution using the method of partial solution. 3.36. The initial temperature of the previous problem is uniformly equal to Ti, ′′′ and the internal heat source is a function of time only, q′′′ = q0 sin(ωt ) . The temperature of both left and right sides are kept at initial temperature. What is the temperature distribution? 3.37. Two semi-infinite bodies A and B are initially at temperature TA,i and TB,i ( TA,i > TB ,i ) are brought into contact at time t = 0. The thermophysical properties of these two bodies are different and the contact thermal resistance can be neglected. The condition of thermal equilibrium requires that the surface temperature of these two bodies change to T0 (between TA,i and TB,i) from t = 0 and remain unchanged throughout the process. Find this surface temperature and transient temperatures in both A and B.
Chapter 3 Heat Conduction 337
3.38. For the heat conduction problem shown in Fig. 3.13, if the boundary condition on the left is constant heat flux, eq. (3.333), will be replaced by ′′ q′′ = q0 x = 0, t > 0 (3.707) Differentiating eqs. (3.332), (3.334), and (3.335) with respect to x, and multiplying the resulting equations by − k , the governing equations and boundary conditions will be in term of heat flux. These equations together with eq. (3.707) – referred to as flux formulation – will have the same form like eqs. (3.332) – (3.335) whose solution can be used to obtain the transient heat flux distribution. The transient temperature distribution for heat conduction in semi-infinite medium with constant heat flux can be obtained by integrating the transient heat flux with respect to x. Obtain the solution of heat conduction in semi-infinite medium with constant heat flux using this approach. 3.39. A semi-infinite body (x > 0) has an initial temperature of Ti . From time t = 0, the surface temperature of the semi-infinite body is increased linearly, i.e., T (0, t ) = Ti + ct , where c is a constant. Obtain the transient temperature distribution using Duhamel’s theorem. 3.40. 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. 3.41. 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. 3.42. Consider transient heat conduction in a 3-D block with dimension of 2L1×2L2×2L3 and an initial temperature of Ti. At time t = 0, the block is immersed into a fluid with temperature T∞. The convective heat transfer coefficients on the three directions are h1, h2 and h3, respectively. It is assumed that the thermal conductivity is independent from the temperature and there is no internal heat source. Obtain the solution of transient temperature distribution using the method of product solution. 3.43. A cylinder with length of 2L and radius of r0 has an initial temperature of Ti. At time t = 0, the cylinder is immersed into a fluid with temperature T∞. The convective heat transfer coefficients on both ends of the cylinder are h1 while on the cylindrical surface is h2. It is assumed that the thermal conductivity is independent from the temperature and there is no internal heat source. Obtain the transient temperature distribution using the method of product solution.
338 Advanced Heat and Mass Transfer
3.44. When an iterative procedure is employed to solve the discretized equation, (3.423), (3.455), (3.459), or (3.464), the iteration will converge if aP ≥ ∑ anb and > is satisfied at least at one point (the subscript nb represents all neighbor points). Discuss the effect on selection of Sp in eq. (3.419) on the numerical solution. 3.45. Show that the heat flux across the faces of the control volume is continuous if harmonic mean conductivity defined in eq. (3.421) at the faces is used. 3.46. Derive the discretized algebraic equation for the grid point at the right boundary that has boundary condition of the second or third kind as shown in Fig. 3.24(b). 3.47. Write a computer program to solve one-dimensional steady-state heat conduction problem discussed in Section 3.4.1. Practice B should be used to discretized the computational domain. Use your computer program to solve heat conduction in an annular fin as shown in Fig. 3.9. 3.48. Use the control volume shown in Fig. 3.24 to derive the discretized equations boundary conditions (for both 2nd and 3rd kinds) at both left and left ends for one-dimensional transient conduction problem. The discretized equations should be able to accommodate f = 0, 0.5 and 1 . If the Practice B is used to discretize the computational domain, how would you modify your discretized equations? 3.49. The stability criterion requires that the grid Fourier number for the explicit scheme less than 0.5 for the internal grid point [see eq. (3.452)]. Analyze the discretized equations for the boundary point that you obtained from the previous problem (for f = 0 ) and obtain the stability criterion for the boundary point. 3.50. In a system with three grid points, E, P, and W, the temperature at E and W are kept at 100 °C and 50 °C, respectively. The initial temperature of the grid point P is 0 °C. Use grid Fourier number of 0.1, 0.25, 0.4, 0.6, 1.0 and eq. (3.453) to get temperature at point P at different times and show that reasonable results cannot e obtained if the grid Fourier number is greater than 0.5. 3.51. Extend the computer program that you developed fro Problem 3.47 to the one-dimensional transient conduction and apply your program to solve transient heat conduction in the annular fin as shown in Fig. 3.9. The initial temperature of the annular fin is at the ambient fluid temperature. 3.52. If the temperatures at all boundaries of a multidimensional heat conduction problem are known, we only need to solve temperatures for the internal points. For the boundary condition of the second and third kind, on the other hand, the discretized equations for the bouindary points must be solved. Thererfore, number of discretized algebraic equations that need
Chapter 3 Heat Conduction 339
solved for different boundary conditions will be different. To overcome this problem, one can use Practice B and treat the energy added to or extracted from the boundary as internal heat source of the first internal cell near the boundary so that the boundary temperature will not appear in the discretized equation for the internal points. This approach – referred to as additional source term method – will enable us to focus on internal points only so that the number of discretized algebraic equations to be solved is the same for all boundary conditions. Derive the discretized equation for the internal grid point P shown in Fig. P3.2. The temperature at boundary point W must not appear in the discretized equation.
N ′′ qB
W
P
E
S
Figure P3.2
3.53. For a two-dimensional melting/solidification problem, the energy balance equation at the interface is
ks ∂Ts ( x, y, t ) ∂T ( x, y, t ) −k = ρ hs 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 + ⎜ ⎟ ⎥ ⎢ k s s − k ⎥ = ρ hs ∂t ∂x ⎦ ⎢ ⎝ ∂y ⎠ ⎥ ⎣ ∂x ⎣ ⎦ x = s ( y, t )
3.54. 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 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. (3.496) – (3.502), which are nondimensional governing equations for solidification
340 Advanced Heat and Mass Transfer
in a semi-infinite region, are also valid for this problem provided that subscripts 1 and 2 represent liquid and solid, respectively. 3.55. The governing equation of the one-region phase change problem, eq. (3.503), 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 one-region problem using the similarity solution. 3.56. 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. 3.57. Solve for melting in a semi-infinite body with convective heating using integral approximate method. The temperature distribution in the liquid phase can be assumed as a second order polynomial function. 3.58. 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. 3.59. The dimensionless governing equations of melting in a subcooled semiinfinite solid with constant heat flux heating were given in eqs. (3.577) – (3.583). Solve the problem using a semi-exact solution: exact solution in the liquid and integral approximate solution in the solid. 3.60. 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, hs = 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? 3.61. 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. (3.605) – and obtain the
Chapter 3 Heat Conduction 341
ablation velocity with the convective and radiation cooling effects considered. 3.62. The surface heat flux in the ablation problem discussed in Section 3.5.3 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. Ablative Material Substrate L
Adiabatic
′′ qw
T1 k1 cp,1 ρ1 y y=0
Tδ qδ
T2 k2 cp,2 ρ2
x=0 y=s
Figure P3.3
x
3.63. 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. P3.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
342 Advanced Heat and Mass Transfer
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. 3.64. 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. 3.65. A line heat source with intensity of q′ (W/m) is located at r = 0 in an infinite solid as shown in Fig. P3.4. The initial temperature of the solid, Ti, 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 ∞ e− z dz . form of Ei[ − r 2 /(4α t )] , where Ei ( z ) = z z
∫
Line heat source
Tm Ti s(t) r
Figure P3.4
0
3.66. One-region melting in a semi-infinite body with constant surface temperature can be described by eqs. (3.503) – (3.506). In a quasistationary 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 (3.503) 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.
Chapter 3 Heat Conduction 343
3.67. 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. 3.68. 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. 3.69. Show that the energy equation for the enthalpy model, eq. (3.641),satisfies the energy equations for both the liquid and solid phases, eqs. (3.485) and (3.487), as well as, energy balance, eq. (3.491). 3.70. In the enthalpy method, we obtained an equation for the numerical solution: h n +1 = j Δt
ρ ( Δ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 ,α )Δt 1 ≤ 2 2 ( Δx ) 3.71. 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. 3.72. 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?
344 Advanced Heat and Mass Transfer
3.73. Show that the enthalpy-temperarture relation in Fig. 3.44 can be represented analytically by eq. (3.667). 3.74. The energy balance at the solid-liquid interface for a melting problem given by eq. (3.479) 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. The heat flux and temperature gradient are related by ∂q′′( x, t ) ∂T ( x, t ) = −k q′′( x, t ) + τ ∂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. 3.75. 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. P3.5). Modify the energy balance that you obtained in the previous problem to account for the contribution of sensible heat to the interfacial energy balance.
Figure P3.5
3.76. Derive the energy balance at the solid-liquid interface for one-dimensional melting problem if the heat conductions in both phases need to be described by the dual-phase lag model. The phase lag times, τq and τT, for both phases can be assumed to be the same. 3.77. For problems with boundary conditions of the second kind, it is desirable to express the energy equation for hyperbolic conduction in term of heat
Chapter 3 Heat Conduction 345
flux. Eliminate temperature from eqs. (3.1) and (3.674) to obtain the energy equation with heat flux as unknown. 3.78. Obtain energy equation for dual-phase lag model in term of heat flux by eliminating temperature from eqs. (3.1) and (3.679). 3.79. Assuming all properties of electrons and lattice are independent from temperatures in the hyperbolic two-step model, combine eqs. (3.693) and (3.682) to obtain a single energy equation for lattice temperature and compare the result with the dual-phase lag model. 3.80. A thin free-standing gold film is subject to a ultrafast laser pulase heating. The laser pulse can be described using eq. (3.700). Write a computer program using the parabolic and hyperbolic two-step model and using the thermophysical properties in Qiu and Tien (1993) to compare the results for these two models. 3.81. Consider the control volume shown in Fig. 3.23 and use approprite numerical scheme to dicrertize the energy equation for the dual-phase lag model, eq. (3.680). 3.82. Nanofluid is obtained by adding a small amount of nanoparticles into base fluid. Pulsed laser is used to heat a layer of nanofluid with thickness of L. The laser can only interact with the nanoparticles and the base fluid is considered to be transparent. If the pulse duration is shorter than the time that it takes for the nanoparticles and base fluid to reach thermal equilibrium, the two-step model can also be used to investigate heat conduction in the nanofluid. What are the governing equations and corresponding boundary conditions for this problem? 3.83. Combine the energy equations for nanoparticles and base fluid that were obtained for the previous problem and obtain a single equation for the base fluid. Compare your equation with the dual-phase lag model and discuss the relationship between the phase lag time and the properties of the nanofluid.
346 Advanced Heat and Mass Transfer