# Numerical Solution of 1-D Steady Conduction

(Difference between revisions)
 Revision as of 21:10, 7 December 2009 (view source)← Older edit Current revision as of 08:42, 3 July 2010 (view source) (6 intermediate revisions not shown) Line 1: Line 1: + ==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)]]: 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)]]: -
$\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = 0 \qquad \qquad( )$
+
$\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = 0 \qquad \qquad(1)$
- (3.418) + - 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,{\rm{ }}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 + 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,$ $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 -
$S = {S_C} + {S_P}T \qquad \qquad( )$
+
$S = {S_C} + {S_P}T \qquad \qquad(2)$
- (3.419) + - where ${S_C}$ is a constant and SP is a coefficient. + where ${S_C}$ is a constant and ${S_P}$ is a coefficient. - Substituting eq. (3.419)  into eq. (3.418), and multiplying the resulting equation by $A(x)$ yields + 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$
$\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. 3.23), one obtains + 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( )$
+
${\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)$
- (3.420) + - 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 + 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}$
$\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.420) become + 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)_e} = {k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}}$
- [[Image:Chapter3_(21).gif|thumb|400 px|alt=Control volume for one-dimensional heat conduction|Figure 1: Control volume for one-dimensional heat conduction]] -
${\left( {kA\frac{{dT}}{{dx}}} \right)_w} = {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}$
${\left( {kA\frac{{dT}}{{dx}}} \right)_w} = {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}$
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 ${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., -
$\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( )$
+
$\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)$
- (3.421) + - For a uniform grid system, eq. (3.421) becomes + 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( )$
+
${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.422) + - Substituting the above equations into eq. (3.420), the following discretized equation is obtained + 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( )$
+
${a_P}{T_P} = {a_E}{T_E} + {a_W}{T_W} + b \qquad \qquad(6)$
- (3.423) + where where -
${a_E} = \frac{{{k_e}{A_e}}}{{{{(\delta x)}_e}}},{\rm{ }}{a_W} = \frac{{{k_w}{A_w}}}{{{{(\delta x)}_w}}} \qquad \qquad( )$
- (3.424) -
${a_P} = {a_W} + {a_E} - {S_P}{A_P}{(\Delta x)_P} \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.425) + -
$b = {S_C}{A_P}{(\Delta x)_P} \qquad \qquad( )$
+
${a_P} = {a_W} + {a_E} - {S_P}{A_P}{(\Delta x)_P} \qquad \qquad(8)$
- (3.426) + - ===Boundary Conditions=== +
$b = {S_C}{A_P}{(\Delta x)_P} \qquad \qquad(9)$
- 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 + ==Boundary Conditions== + [[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''']] - [[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]] + 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 - + - 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 +
${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$
${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$
Line 69: Line 58: which can be expressed as which can be expressed as -
${a_1}{T_1} = {a_E}{T_2} + b \qquad \qquad( )$
+
${a_1}{T_1} = {a_E}{T_2} + b \qquad \qquad(10)$
- (3.427) + where where -
${a_E} = \frac{{{k_1}{A_1}}}{{{{(\delta x)}_1}}} \qquad \qquad( )$
+
${a_E} = \frac{{{k_1}{A_1}}}{{{{(\delta x)}_1}}} \qquad \qquad(11)$
- (3.428) + -
${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} \qquad \qquad( )$
+
${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} \qquad \qquad(12)$
- (3.429) + -
$b = {q''_B}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad( )$
+
$b = {q''_B}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(13)$
- (3.430) + If the left boundary is subject to a convective condition, we have If the left boundary is subject to a convective condition, we have -
${q''_B} = h({T_f} - {T_1}) \qquad \qquad( )$
+
${q''_B} = h({T_f} - {T_1}) \qquad \qquad(14)$
- (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 + 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( )$
+
${a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} + h{A_1} \qquad \qquad(15)$
- (3.432) + -
$b = h{T_f}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad( )$
+
$b = h{T_f}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(16)$
- (3.433) + The boundary condition on the right hand side can be obtained by following a similar procedure (see Problem 3.46). 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 ${(\Delta 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 ${(\Delta x)_M}$ in the discretized equation for Practice A. + 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}$ in the discretized equation for Practice A. + + ==Solution of Discretized Equations== - [[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]] + 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 - + - ===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 [[#References|(Patankar, 1980)]].  The discretized equations for every point, eq. (3.423), 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)$
- + -
${a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}{T_{i - 1}} + {d_i},{\rm{ }}i = 1,2,...M \qquad \qquad( )$
+ - (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 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 Line 114: Line 93: To get the temperature, we hope that our temperature at the ${i_{th}}$ grid is related to the neighbor on the right, i.e., To get the temperature, we hope that our temperature at the ${i_{th}}$ grid is related to the neighbor on the right, i.e., -
${T_i} = {P_i}{T_{i + 1}} + {Q_i} \qquad \qquad( )$
+
${T_i} = {P_i}{T_{i + 1}} + {Q_i} \qquad \qquad(18)$
- (3.435) + - For the ${(i - 1)^{{\rm{th}}}}$ grid point, eq. (3.435) becomes + For the ${(i - 1)^{{\rm{th}}}}$ grid point, eq. (18) becomes -
${T_{i - 1}} = {P_{i - 1}}{T_i} + {Q_{i - 1}} \qquad \qquad( )$
+
${T_{i - 1}} = {P_{i - 1}}{T_i} + {Q_{i - 1}} \qquad \qquad(19)$
- (3.436) + - Substituting eq. (3.436) into eq. (3.434), we get + 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( )$
+
${a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}({P_{i - 1}}{T_i} + {Q_{i - 1}}) + {b_i} \qquad \qquad(20)$
- (3.437) + which can be rearranged to 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( )$
+
${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)$
- (3.438) + - Comparing eq. (3.435) and (3.438) yields + 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( )$
+
${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.439) + - For the first grid point at the left $(i = 1)$, ${c_1} = 0$ and eq. (3.439) becomes + For the first grid point at the left $(i = 1)$, ${c_1} = 0$ and eq. (22) becomes -
${P_1} = \frac{{{b_1}}}{{{a_1}}},{\rm{ }}{Q_1} = \frac{{{d_1}}}{{{a_1}}} \qquad \qquad( )$
+
${P_1} = \frac{{{b_1}}}{{{a_1}}},{\rm{ }}{Q_1} = \frac{{{d_1}}}{{{a_1}}} \qquad \qquad(23)$
- (3.440) + - For the last point at the right $(i = M)$, ${b_M} = 0$ and it follows that ${P_M} = 0$. Equation (3.435) for $i = M$ becomes + For the last point at the right $(i = M)$, ${b_M} = 0$ and it follows that ${P_M} = 0$. Equation (18) for $i = M$ becomes -
${T_N} = {Q_N} \qquad \qquad( )$
+
${T_N} = {Q_N} \qquad \qquad(24)$
- (3.441) + - The solution of the algebraic equations can thus be completed in the following steps: (1) obtain ${P_1}$ and ${Q_1}$ from eq. (3.440), (2) calculate ${P_i}$ and ${Q_i}$ 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). + The solution of the algebraic equations can thus be completed in the following steps: (1) obtain ${P_1}$ 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). 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. 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== ==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). ==Further Reading== ==Further Reading== ==External Links== ==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$

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

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