# Numerical Solution of 1-D Steady Conduction

(Difference between revisions)
 Revision as of 00:46, 4 December 2009 (view source)← Older edit Current revision as of 08:42, 3 July 2010 (view source) (8 intermediate revisions not shown) Line 1: Line 1: - ==3.5 Melting and Solidification== + ==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 [[#References|(Tao, 2001)]]: - ===3.5.1 Introduction=== + - +
$\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = 0 \qquad \qquad(1)$
- 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 (${C_{18}}{H_{38}}$) 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, ${T_i}$, equals the melting point, ${T_m}$, 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,  ${T_0}$ , is greater than ${T_m}$ for melting and is below ${T_m}$ for solidification. + where $A(x)$ is the area of heat conduction. Equation (1) is valid for heat conduction in Cartesian ($A(x) = 1$), cylindrical ($x = r,{\rm{ }}A(r) = r$), and spherical ($x = r, [itex]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. (1) 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 - [[Image:Chapter3_(3).png|thumb|400 px|alt=One-region melting and solidification|
(a) Melting () (b) Solidification ()

Figure 1: One-region melting and solidification]] +
$S = {S_C} + {S_P}T \qquad \qquad(2)$
- [[Image:Chapter3_(3).png|thumb|400 px|alt=Two-region melting and solidification.|
(a) Melting () (b) Solidification ()

Figure 2: Two-region melting and solidification.]] + where ${S_C}$ is a constant and ${S_P}$ is a coefficient. - For the melting process, if the initial temperature of the PCM, ${T_i}$, is below the melting point of the PCM, ${T_m}$, (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 one-dimensional two-region melting and solidification problems. + Substituting eq. (2)  into eq. (1), and multiplying the resulting equation by $A(x)$ yields - For a multi-component PCM, the solid-liquid phase change process occurs over a range of temperatures (${T_m}$1, ${T_{m2}}$), instead of a single temperature. The PCM is liquid if its temperature is above ${T_{m2}}$ and solid when its temperature is below ${T_{m1}}$. Between the solid and liquid phases there is a ''mushy zone'' where the temperature falls within the phase change temperature range (${T_{m1}}$, ${T_{m2}}). 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. + [itex]\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + ({S_C} + {S_P}T)A(x) = 0$
+ [[Image:Chapter3_(21).gif|thumb|400 px|alt=Control volume for one-dimensional heat conduction|'''Figure 1: Control volume for one-dimensional heat conduction''']] + Integrating the above equation over the control volume P (shaded area in Fig. 1), one obtains - 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. +
${\left( {kA\frac{{dT}}{{dx}}} \right)_e} - {\left( {kA\frac{{dT}}{{dx}}} \right)_w} + \int_w^e {({S_C} + {S_P}T)Adx} = 0 \qquad \qquad(3)$
- 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 + 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) becomes - [[Image:Chapter3_(4).png|thumb|400 px|alt=Solidification of a multicomponent PCM|Figure 1: Solidification of a multicomponent PCM]] +
$\int_w^e {({S_C} + {S_P}T)Adx} = {S_C}{A_P}{(\Delta x)_P} + {S_P}{T_P}{A_P}{(\Delta 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) become - temperature of the interface, which is at the melting point, ${T_m}.$ Therefore, the boundary conditions at the interface can be expressed as +
${\left( {kA\frac{{dT}}{{dx}}} \right)_e} = {k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}}$
+ +
${\left( {kA\frac{{dT}}{{dx}}} \right)_w} = {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}$
-
${T_\ell }(x,t) = {T_s}(x,t) = {T_m},\quad \quad x = s(t) \qquad \qquad( )$
(3.472) + where ${k_e}$ and ${k_w} 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., - where [itex]{T_\ell }(x,t)$ and ${T_s}(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, ${\rho _\ell },$ is less than the density of the solid paraffin, ${\rho _s}.$ When water freezes, however, its volume increases, so the density of ice is less than that of liquid water. +
$\frac{{{{(\delta x)}_e}}}{{{k_e}}} = \frac{{{{(\delta x)}_{e - }}}}{{{k_P}}} + \frac{{{{(\delta x)}_{e + }}}}{{{k_E}}},{\rm{ }}\frac{{{{(\delta x)}_w}}}{{{k_w}}} = \frac{{{{(\delta x)}_{w - }}}}{{{k_W}}} + \frac{{{{(\delta x)}_{w + }}}}{{{k_P}}} \qquad \qquad(4)$
- 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., ${\rho _s} > {\rho _\ell }$) 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_\rho }$. The density change must satisfy the conservation of mass at the interface, i.e., + - + -
${\rho _s}({u_p} - {u_\rho }) = {\rho _\ell }{u_p} \qquad \qquad( )$
+ - (3.473) + - The extra velocity induced by the density change can be obtained by rearranging eq. (3.473) as + For a uniform grid system, eq. (4) becomes -
${u_\rho } = \frac{{{\rho _s} - {\rho _\ell }}}{{{\rho _s}}}{u_p} \qquad \qquad( )$
+
${k_e} = \frac{{2{k_P}{K_E}}}{{{k_P} + {K_E}}},{\rm{ }}{k_w} = \frac{{2{k_W}{K_P}}}{{{k_W} + {K_P}}} \qquad \qquad(5)$
- (3.474) + - which is also valid for case ${\rho _s} > {\rho _\ell },$ except the extra velocity becomes negative. Another necessary boundary condition is the energy balance at the solid-liquid interface. If the enthalpy of the liquid and solid phases at the melting point are ${h_\ell }$ and ${h_s}$, the energy balance at the solid-liquid interface can be expressed as: + Substituting the above equations into eq. (3), the following discretized equation is obtained - + -
${q''_\ell } - {q''_s} = {\rho _\ell }{u_p}{h_\ell } - {\rho _s}({u_p} - {u_\rho }){h_s}\quad \quad x = s(t) \qquad \qquad( )$
+ - (3.475) + - where ${q''_\ell }$ and ${q''_s}$ 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 +
${a_P}{T_P} = {a_E}{T_E} + {a_W}{T_W} + b \qquad \qquad(6)$
- + -
${q''_\ell } - {q''_s} = {\rho _\ell }{h_{s\ell }}\frac{{ds}}{{dt}}\quad \quad x = s(t) \qquad \qquad( )$
+ - (3.476) + - where ${h_{s\ell }} = {h_\ell } - {h_s}$ 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: + where -
${q''_\ell } = \mathop {\left. { - {k_\ell }\frac{{\partial {T_\ell }(x,t)}}{{\partial x}}} \right|}{x = s(t)} \qquad \qquad( )$
+
${a_E} = \frac{{{k_e}{A_e}}}{{{{(\delta x)}_e}}},{\rm{ }}{a_W} = \frac{{{k_w}{A_w}}}{{{{(\delta x)}_w}}} \qquad \qquad(7)$
- (3.477) + -
${q''_s} = \mathop {\left. { - {k_s}\frac{{\partial {T_s}(x,t)}}{{\partial x}}} \right|}{x = s(t)} \qquad \qquad( )$
+
${a_P} = {a_W} + {a_E} - {S_P}{A_P}{(\Delta x)_P} \qquad \qquad(8)$
- (3.478) + - 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., +
$b = {S_C}{A_P}{(\Delta x)_P} \qquad \qquad(9)$
-
${k_s}\frac{{\partial {T_s}(x,t)}}{{\partial x}} - {k_\ell }\frac{{\partial {T_\ell }(x,t)}}{{\partial x}} = {\rho _\ell }{h_{s\ell }}\frac{{ds(t)}}{{dt}}\quad \quad x = s(t) \qquad \qquad( ) + ==Boundary Conditions== - (3.479) + [[Image:Chapter3_(22).gif|thumb|400 px|alt=Control volumes at boundaries for Practice A| (a) Left boundary (b) Right boundary '''Figure 2: Control volumes at boundaries for Practice A''']] + [[Image:Chapter3_(23).gif|thumb|400 px|alt=Control volumes at boundaries for Practice A| (a) Left boundary (b) Right boundary '''Figure 3: Control volumes at boundaries for Practice B''']] - For a solidification process, the energy balance equation at the interface can be obtained by a similar procedure: + 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 internal control volumes (see Fig. 2). If the heat flux at the left boundary is known, the energy balance for the control volume at the left boundary is - + - [itex]{k_s}\frac{{\partial {T_s}(x,t)}}{{\partial x}} - {k_\ell }\frac{{\partial {T_\ell }(x,t)}}{{\partial x}} = {\rho _s}{h_{s\ell }}\frac{{ds(t)}}{{dt}}\quad \quad x = s(t) \qquad \qquad( )$
+ - (3.480) + - 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. +
${q''_B}{A_1} + {k_1}{A_1}\frac{{{T_2} - {T_1}}}{{{{(\delta x)}_1}}} + ({S_c} + {S_P}{T_1}){A_1}{(\Delta x)_1} = 0$
- 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 + - + which can be expressed as -
${T_\ell }(x,y,t) = {T_s}(x,y,t) = {T_m}\quad \quad x = s(y,t) \qquad \qquad( )$
+ - (3.481) + -
${k_s}\frac{{\partial {T_s}(x,y,t)}}{{\partial n}} - {k_\ell }\frac{{\partial {T_\ell }(x,y,t)}}{{\partial n}} = {\rho _\ell }{h_{s\ell }}{v_n}\quad \quad x = s(y,t) \qquad \qquad( )$
+
${a_1}{T_1} = {a_E}{T_2} + b \qquad \qquad(10)$
- (3.482) + - where ${\mathbf{n}}$ is a unit vector along the normal direction of the solid-liquid interface, and ${v_n}$ is the solid-liquid interface velocity along the n-direction. + where - 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 [[#References|(see Problem 3.55; Ozisik, 1993)]]: + - + -
$\left[ {1 + \mathop {\left( {\frac{{\partial s}}{{\partial y}}} \right)}^2 } \right]\left[ {{k_s}\frac{{\partial {T_s}}}{{\partial x}} - {k_\ell }\frac{{\partial {T_\ell }}}{{\partial x}}} \right] = {\rho _\ell }{h_{s\ell }}\frac{{\partial s}}{{\partial t}}\quad \quad x = s(y,t) \qquad \qquad( )$
+ - (3.483) + - Similarly, for a three-dimensional melting problem with an interface described by +
${a_E} = \frac{{{k_1}{A_1}}}{{{{(\delta x)}_1}}} \qquad \qquad(11)$
- + -
$z = s(x,y,t)$
+ - + - the energy balance at the interface is + -
$\left[ {1 + \mathop {\left( {\frac{{\partial s}}{{\partial x}}} \right)}^2 + \mathop {\left( {\frac{{\partial s}}{{\partial y}}} \right)}^2 } \right]\left[ {{k_s}\frac{{\partial {T_s}}}{{\partial x}} - {k_\ell }\frac{{\partial {T_\ell }}}{{\partial x}}} \right] = {\rho _\ell }{h_{s\ell }}\frac{{\partial s}}{{\partial t}}\quad \quad z = s(x,y,t) \qquad \qquad( )$
+
${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} \qquad \qquad(12)$
- (3.484) + - For solidification problems, it is necessary to replace the liquid density {\rho _\ell } in eqs. (3.483) and (3.484) with the solid-phase density ${\rho _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. + [itex]b = {q''_B}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(13) - ==3.5.2 Exact Solution== + If the left boundary is subject to a convective condition, we have - + - ===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 {T_i},$ which exceeds the melting point ${T_m},$ is in a half-space $x > 0.$ At time $t = 0,$ the temperature at the boundary $x = 0$ is suddenly decreased to a temperature ${T_0},$ which is below the melting point of the PCM. Solidification occurs from the time $t = 0. This is a two-region solidification problem as the temperatures of both the liquid and solid phases are unknown and must be determined. It is assumed that the densities of the PCM for both phases are the same. Natural convection in the liquid phase is neglected, and therefore the heat transfer mechanism in both phases is pure conduction. + [itex]{q''_B} = h({T_f} - {T_1}) \qquad \qquad(14)$
- The temperature in the solid phase must satisfy + -
$\frac{{{\partial ^2}{T_1}}}{{\partial {x^2}}} = \frac{1}{{{\alpha _1}}}\frac{{\partial {T_1}}}{{\partial t}}\quad \quad 0 < x < s(t),\quad t > 0 \qquad \qquad( )$
+ It can be shown that eqs. (10) and (11) are still valid but eqs. (12)  and (13) should be modified to - (3.485) + -
${T_1}(x,t) = {T_0}\quad \quad x = 0,\quad t > 0 \qquad \qquad( )$
+
${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} + h{A_1} \qquad \qquad(15)$
- (3.486) + - For the liquid phase, the governing equations are +
$b = h{T_f}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(16)$
-
$\frac{{{\partial ^2}{T_2}}}{{\partial {x^2}}} = \frac{1}{{{\alpha _2}}}\frac{{\partial {T_2}}}{{\partial t}}\quad \quad s(t) < x < \infty ,\quad t > 0 \qquad \qquad( )$
+ The boundary condition on the right hand side can be obtained by following a similar procedure (see Problem 3.46). - (3.487) + 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), thus the discretized equation for the boundary condition at the left side can be obtained by setting the control volume size ${(\Delta x)_1}$ to zero in eqs. (12) – (13) or eqs. (15) – (16). Similarly, the discretized equation for the boundary condition on the right side for Practice B can be obtained by setting {(\Delta x)_M}[/itex] in the discretized equation for Practice A. + + ==Solution of Discretized Equations== -
${T_2}(x,t) \to {T_i}\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad( )$
+ When the matrix for the coefficient for discretized algebraic equations of one-dimensional 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 [[#References|(Patankar, 1980)]]. The discretized equations for every point, eq. (6), can be rewritten as - (3.488) + -
${T_2}(x,t) = {T_i}\quad \quad x > 0,\quad t = 0 \qquad \qquad( )$
+
${a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}{T_{i - 1}} + {d_i},{\rm{ }}i = 1,2,...M \qquad \qquad(17)$
- (3.489) + - The boundary conditions at the interface are + 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 -
${T_1}(x,t) = {T_2}(x,t) = {T_m}\quad \quad x = s(t),\quad t > 0 \qquad \qquad( )$
+
${c_1} = 0,{\rm{ }}{b_M} = 0$
- (3.490) + + To get the temperature, we hope that our temperature at the {i_{th}}[/itex] grid is related to the neighbor on the right, i.e., + +
${T_i} = {P_i}{T_{i + 1}} + {Q_i} \qquad \qquad(18)$
-
${k_1}\frac{{\partial {T_1}}}{{\partial x}} - {k_2}\frac{{\partial {T_2}}}{{\partial x}} = \rho {h_{s\ell }}\frac{{ds}}{{dt}}\quad \quad x = s(t),\quad t > 0 \qquad \qquad( )$
+ For the ${(i - 1)^{{\rm{th}}}}$ grid point, eq. (18) becomes - (3.491) + +
${T_{i - 1}} = {P_{i - 1}}{T_i} + {Q_{i - 1}} \qquad \qquad(19)$
- 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 + Substituting eq. (19) into eq. (17), we get -
$\frac{{\partial {T_1}}}{{\partial x}} \sim \frac{{{T_m} - {T_0}}}{s}$
+
${a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}({P_{i - 1}}{T_i} + {Q_{i - 1}}) + {b_i} \qquad \qquad(20)$
- + -
$\frac{{\partial {T_2}}}{{\partial x}} \sim \frac{{{T_i} - {T_m}}}{s}$
+ - + - The scale of the interface velocity is + -
$\frac{{ds}}{{dt}} \sim \frac{s}{t}$
+ which can be rearranged to - + - Substituting the above equations into eq. (3.491), one obtains + -
${k_1}\frac{{{T_m} - {T_0}}}{s} - {k_2}\frac{{{T_i} - {T_m}}}{s} \sim \rho {h_{s\ell }}\frac{s}{t}$
+
${T_i} = \frac{{{b_i}}}{{{a_i} - {c_i}{P_{i - 1}}}}{T_{i + 1}} + \frac{{{d_i} + {c_i}{Q_{i - 1}}}}{{{a_i} - {c_i}{P_{i - 1}}}} \qquad \qquad(21)$
- + - The scale of the location of the solid-liquid interface is obtained by rearranging the above equation, i.e., + - + -
$\frac{{{s^2}}}{{{\alpha _1}t}} \sim {\rm{Ste}}\left( {1 - \frac{{{k_2}}}{{{k_1}}}\frac{{{T_i} - {T_m}}}{{{T_m} - {T_0}}}} \right) \qquad \qquad( )$
+ - (3.492) + - where + Comparing eq. (18) and (21) yields -
${\rm{Ste}} = \frac{{{c_{p1}}({T_m} - {T_0})}}{{{h_{s\ell }}}} \qquad \qquad( )$
+
${P_i} = \frac{{{b_i}}}{{{a_i} - {c_i}{P_{i - 1}}}},{\rm{ }}{Q_i} = \frac{{{d_i} + {c_i}{Q_{i - 1}}}}{{{a_i} - {c_i}{P_{i - 1}}}} \qquad \qquad(22)$
- (3.493) + - is the Stefan number. Named after J. Stefan, a pioneer in discovery of the solid-liquid 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_{{p_1}}}({T_m} - {T_0}),$ to latent heat, ${h_{s\ell }}$. For a latent heat thermal energy storage system, the Stefan number is usually very small because the temperature difference in such a system is small, while the latent heat ${h_{s\ell }}$ 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. + For the first grid point at the left $(i = 1)$, ${c_1} = 0$ and eq. (22) becomes - It can be seen from eq. (3.492) that the effect of heat conduction in the liquid phase can be neglected if ${T_i} - {T_m} \ll {T_m} - {T_0}$ or ${k_2} \ll {k_1}$. In that case, eq. (3.492) can be simplified as + -
$\frac{{{s^2}}}{{{\alpha _1}t}} \sim {\rm{Ste}} \qquad \qquad( )$
+
${P_1} = \frac{{{b_1}}}{{{a_1}}},{\rm{ }}{Q_1} = \frac{{{d_1}}}{{{a_1}}} \qquad \qquad(23)$
- (3.494) + - + - which means that the interfacial velocity increases with increasing $\Delta T = \left| {{T_m} - {T_0}} \right|$ or decreasing ${h_{s\ell }}$. + - + - ===Dimensionless Form of the Governing Equations=== + - + - The governing eqs. (3.485) – (3.491) can be nondimensionalized by introducing the following dimensionless variables: + - + -
$\left. \begin{array}{l} + - \theta = \frac{{{T_m} - T}}{{{T_m} - {T_0}}}\,\,\,\,{\theta _i} = \frac{{{T_m} - {T_i}}}{{{T_m} - {T_0}}}\,\,\,X = \frac{x}{L}\quad S = \frac{s}{L}\quad \tau = \frac{{{\alpha _1}t}}{{{L^2}}} \\ + - {N_\alpha } = \frac{{{\alpha _2}}}{{{\alpha _1}}}\,\,\,\,\,\,{N_k} = \frac{{{k_2}}}{{{k_1}}}\,\,\,\,{\rm{Ste}} = \frac{{{c_{{p_1}}}({T_m} - {T_0})}}{{{h_{s\ell }}}} \\ + - \end{array} \right\} \qquad \qquad( )$
+ - (3.495) + - + - 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: + - + -
$\frac{{{\partial ^2}{\theta _1}}}{{\partial {X^2}}} = \frac{{\partial {\theta _1}}}{{\partial \tau }}\quad \quad 0 < X < S(\tau ),\quad \tau > 0 \qquad \qquad( )$
+ - (3.496) + -
${\theta _1}(X,\tau ) = 1\quad \quad X = 0,\quad \tau > 0 \qquad \qquad( ) + For the last point at the right (i = M)$, ${b_M} = 0 and it follows that [itex]{P_M} = 0$. Equation (18) for $i = M becomes - (3.497) + - [itex]\frac{{{\partial ^2}{\theta _2}}}{{\partial {X^2}}} = \frac{1}{{{N_\alpha }}}\frac{{\partial {\theta _2}}}{{\partial \tau }}\quad \quad S(\tau ) < X < \infty ,\quad \tau > 0 \qquad \qquad( )$
+
${T_N} = {Q_N} \qquad \qquad(24)$
- (3.498) + -
${\theta _2}(X,\tau ) \to {\theta _i}\quad \quad X \to \infty ,\quad \tau > 0 \qquad \qquad( )$
+ The solution of the algebraic equations can thus be completed in the following steps: (1) obtain {P_1}[/itex] and ${Q_1}$ from eq. (23), (2) calculate ${P_i}$ and ${Q_i} from eq. (22), (3) obtain temperature at the right side by eq. (24), and (4) calculate temperatures for all grid points by eq. (18). - (3.499) + - [itex]{\theta _2}(X,\tau ) = {\theta _i}\quad \quad X > 0,\quad \tau = 0 \qquad \qquad( )$
+ 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.500) + -
${\theta _1}(X,\tau ) = {\theta _2}(X,\tau ) = 0\quad \quad X = S(\tau ),\quad \tau > 0 \qquad \qquad( )$
+ ==References== - (3.501) + -
$- \frac{{\partial {\theta _1}}}{{\partial X}} + {N_k}\frac{{\partial {\theta _2}}}{{\partial X}} = \frac{1}{{Ste}}\frac{{dS}}{{d\tau }}\quad \quad X = S(\tau ),\quad \tau > 0 \qquad \qquad( )$
+ Patankar, S.V., 1980, ''Numerical Heat Transfer and Fluid Flow'', Hemisphere, Washington, DC. - (3.502) + - [[Image:Chapter3_(4).jpg|thumb|400 px|alt=Dimensionless temperature distribution in the PCM|Figure : Dimensionless temperature distribution in the PCM]] + Tao, W.Q., 2001, ''Numerical Heat Transfer'', 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese). - 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. + ==Further Reading== - ===Exact Solution of the One-Region Problem=== + ==External Links==

## 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):

$\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = 0 \qquad \qquad(1)$

where A(x) is the area of heat conduction. Equation (1) is valid for heat conduction in Cartesian (A(x) = 1), cylindrical (x = r,A(r) = r), and spherical (x = r, A(r) = r2) 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. (1) 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 = {S_C} + {S_P}T \qquad \qquad(2)$

where SC is a constant and SP is a coefficient.

Substituting eq. (2) into eq. (1), and multiplying the resulting equation by A(x) yields

$\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + ({S_C} + {S_P}T)A(x) = 0$
Figure 1: Control volume for one-dimensional heat conduction

Integrating the above equation over the control volume P (shaded area in Fig. 1), one obtains

${\left( {kA\frac{{dT}}{{dx}}} \right)_e} - {\left( {kA\frac{{dT}}{{dx}}} \right)_w} + \int_w^e {({S_C} + {S_P}T)Adx} = 0 \qquad \qquad(3)$

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) becomes

$\int_w^e {({S_C} + {S_P}T)Adx} = {S_C}{A_P}{(\Delta x)_P} + {S_P}{T_P}{A_P}{(\Delta 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) become

${\left( {kA\frac{{dT}}{{dx}}} \right)_e} = {k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}}$
${\left( {kA\frac{{dT}}{{dx}}} \right)_w} = {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_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.,

$\frac{{{{(\delta x)}_e}}}{{{k_e}}} = \frac{{{{(\delta x)}_{e - }}}}{{{k_P}}} + \frac{{{{(\delta x)}_{e + }}}}{{{k_E}}},{\rm{ }}\frac{{{{(\delta x)}_w}}}{{{k_w}}} = \frac{{{{(\delta x)}_{w - }}}}{{{k_W}}} + \frac{{{{(\delta x)}_{w + }}}}{{{k_P}}} \qquad \qquad(4)$

For a uniform grid system, eq. (4) becomes

${k_e} = \frac{{2{k_P}{K_E}}}{{{k_P} + {K_E}}},{\rm{ }}{k_w} = \frac{{2{k_W}{K_P}}}{{{k_W} + {K_P}}} \qquad \qquad(5)$

Substituting the above equations into eq. (3), the following discretized equation is obtained

${a_P}{T_P} = {a_E}{T_E} + {a_W}{T_W} + b \qquad \qquad(6)$

where

${a_E} = \frac{{{k_e}{A_e}}}{{{{(\delta x)}_e}}},{\rm{ }}{a_W} = \frac{{{k_w}{A_w}}}{{{{(\delta x)}_w}}} \qquad \qquad(7)$
${a_P} = {a_W} + {a_E} - {S_P}{A_P}{(\Delta x)_P} \qquad \qquad(8)$
$b = {S_C}{A_P}{(\Delta x)_P} \qquad \qquad(9)$

## Boundary Conditions

(a) Left boundary (b) Right boundary

Figure 2: Control volumes at boundaries for Practice A
(a) Left boundary (b) Right boundary

Figure 3: Control volumes at boundaries for Practice B

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 internal control volumes (see Fig. 2). If the heat flux at the left boundary is known, the energy balance for the control volume at the left boundary is

${q''_B}{A_1} + {k_1}{A_1}\frac{{{T_2} - {T_1}}}{{{{(\delta x)}_1}}} + ({S_c} + {S_P}{T_1}){A_1}{(\Delta x)_1} = 0$

which can be expressed as

${a_1}{T_1} = {a_E}{T_2} + b \qquad \qquad(10)$

where

${a_E} = \frac{{{k_1}{A_1}}}{{{{(\delta x)}_1}}} \qquad \qquad(11)$
${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} \qquad \qquad(12)$
$b = {q''_B}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(13)$

If the left boundary is subject to a convective condition, we have

${q''_B} = h({T_f} - {T_1}) \qquad \qquad(14)$

It can be shown that eqs. (10) and (11) are still valid but eqs. (12) and (13) should be modified to

${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} + h{A_1} \qquad \qquad(15)$
$b = h{T_f}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(16)$

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), 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. (12) – (13) or eqs. (15) – (16). 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.

## Solution of Discretized Equations

When the matrix for the coefficient for discretized algebraic equations of one-dimensional 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. (6), can be rewritten as

${a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}{T_{i - 1}} + {d_i},{\rm{ }}i = 1,2,...M \qquad \qquad(17)$

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

${T_i} = {P_i}{T_{i + 1}} + {Q_i} \qquad \qquad(18)$

For the (i − 1)th grid point, eq. (18) becomes

${T_{i - 1}} = {P_{i - 1}}{T_i} + {Q_{i - 1}} \qquad \qquad(19)$

Substituting eq. (19) into eq. (17), we get

${a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}({P_{i - 1}}{T_i} + {Q_{i - 1}}) + {b_i} \qquad \qquad(20)$

which can be rearranged to

${T_i} = \frac{{{b_i}}}{{{a_i} - {c_i}{P_{i - 1}}}}{T_{i + 1}} + \frac{{{d_i} + {c_i}{Q_{i - 1}}}}{{{a_i} - {c_i}{P_{i - 1}}}} \qquad \qquad(21)$

Comparing eq. (18) and (21) yields

${P_i} = \frac{{{b_i}}}{{{a_i} - {c_i}{P_{i - 1}}}},{\rm{ }}{Q_i} = \frac{{{d_i} + {c_i}{Q_{i - 1}}}}{{{a_i} - {c_i}{P_{i - 1}}}} \qquad \qquad(22)$

For the first grid point at the left (i = 1), c1 = 0 and eq. (22) becomes

${P_1} = \frac{{{b_1}}}{{{a_1}}},{\rm{ }}{Q_1} = \frac{{{d_1}}}{{{a_1}}} \qquad \qquad(23)$

For the last point at the right (i = M), bM = 0 and it follows that PM = 0. Equation (18) for i = M becomes

${T_N} = {Q_N} \qquad \qquad(24)$

The solution of the algebraic equations can thus be completed in the following steps: (1) obtain P1 and Q1 from eq. (23), (2) calculate Pi and Qi from eq. (22), (3) obtain temperature at the right side by eq. (24), and (4) calculate temperatures for all grid points by eq. (18).

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.

## References

Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.

Tao, W.Q., 2001, Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese).