# Numerical Solution of 1-D Steady Conduction

### From Thermal-FluidsPedia

Line 1: | Line 1: | ||

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

- | <center><math>\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = 0 \qquad \qquad( )</math></center> | + | <center><math>\frac{1}{{A(x)}}\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + S = 0 \qquad \qquad(1)</math></center> |

(3.418) | (3.418) | ||

where <math>A(x)</math> is the area of heat conduction. Equation (3.418) is valid for heat conduction in Cartesian (<math>A(x) = 1</math>), cylindrical (<math>x = r,{\rm{ }}A(r) = r</math>), and spherical (<math>x = r,</math> <math>A(r) = {r^2}</math>) coordinate systems, as well as heat transfer from extended surfaces (<math>A(x)</math> is the cross-sectional area of the fin). The source term, <math>S</math>, 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 <math>A(x)</math> is the area of heat conduction. Equation (3.418) is valid for heat conduction in Cartesian (<math>A(x) = 1</math>), cylindrical (<math>x = r,{\rm{ }}A(r) = r</math>), and spherical (<math>x = r,</math> <math>A(r) = {r^2}</math>) coordinate systems, as well as heat transfer from extended surfaces (<math>A(x)</math> is the cross-sectional area of the fin). The source term, <math>S</math>, 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 | ||

- | <center><math>S = {S_C} + {S_P}T \qquad \qquad( )</math></center> | + | <center><math>S = {S_C} + {S_P}T \qquad \qquad(2)</math></center> |

(3.419) | (3.419) | ||

where <math>{S_C}</math> is a constant and SP is a coefficient. | where <math>{S_C}</math> is a constant and SP is a coefficient. | ||

- | Substituting eq. ( | + | Substituting eq. (2) into eq. (1), and multiplying the resulting equation by <math>A(x)</math> yields |

<center><math>\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + ({S_C} + {S_P}T)A(x) = 0</math></center> | <center><math>\frac{d}{{dx}}\left( {kA(x)\frac{{dT}}{{dx}}} \right) + ({S_C} + {S_P}T)A(x) = 0</math></center> | ||

Line 17: | Line 17: | ||

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. 3.23), one obtains | ||

- | <center><math>{\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( )</math></center> | + | <center><math>{\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)</math></center> |

(3.420) | (3.420) | ||

- | Assuming the temperature <math>T</math> and conduction area <math>A</math> for each control volume are represented by their value at a grid point, the integral in eq. (3 | + | Assuming the temperature <math>T</math> and conduction area <math>A</math> for each control volume are represented by their value at a grid point, the integral in eq. (3) becomes |

<center><math>\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}</math></center> | <center><math>\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}</math></center> | ||

- | 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 | + | 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 |

<center><math>{\left( {kA\frac{{dT}}{{dx}}} \right)_e} = {k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}}</math></center> | <center><math>{\left( {kA\frac{{dT}}{{dx}}} \right)_e} = {k_e}{A_e}\frac{{{T_E} - {T_P}}}{{{{(\delta x)}_e}}}</math></center> | ||

Line 34: | Line 34: | ||

where <math>{k_e}</math> and <math>{k_w}</math> 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 <math>{k_e}</math> and <math>{k_w}</math> 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., | ||

- | <center><math>\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( )</math></center> | + | <center><math>\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)</math></center> |

(3.421) | (3.421) | ||

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

- | <center><math>{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( )</math></center> | + | <center><math>{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)</math></center> |

(3.422) | (3.422) | ||

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

- | <center><math>{a_P}{T_P} = {a_E}{T_E} + {a_W}{T_W} + b \qquad \qquad( )</math></center> | + | <center><math>{a_P}{T_P} = {a_E}{T_E} + {a_W}{T_W} + b \qquad \qquad(6)</math></center> |

(3.423) | (3.423) | ||

where | where | ||

- | <center><math>{a_E} = \frac{{{k_e}{A_e}}}{{{{(\delta x)}_e}}},{\rm{ }}{a_W} = \frac{{{k_w}{A_w}}}{{{{(\delta x)}_w}}} \qquad \qquad( )</math></center> | + | <center><math>{a_E} = \frac{{{k_e}{A_e}}}{{{{(\delta x)}_e}}},{\rm{ }}{a_W} = \frac{{{k_w}{A_w}}}{{{{(\delta x)}_w}}} \qquad \qquad(7)</math></center> |

(3.424) | (3.424) | ||

- | <center><math>{a_P} = {a_W} + {a_E} - {S_P}{A_P}{(\Delta x)_P} \qquad \qquad( )</math></center> | + | <center><math>{a_P} = {a_W} + {a_E} - {S_P}{A_P}{(\Delta x)_P} \qquad \qquad(8)</math></center> |

(3.425) | (3.425) | ||

- | <center><math>b = {S_C}{A_P}{(\Delta x)_P} \qquad \qquad( )</math></center> | + | <center><math>b = {S_C}{A_P}{(\Delta x)_P} \qquad \qquad(9)</math></center> |

(3.426) | (3.426) | ||

Line 69: | Line 69: | ||

which can be expressed as | which can be expressed as | ||

- | <center><math>{a_1}{T_1} = {a_E}{T_2} + b \qquad \qquad( )</math></center> | + | <center><math>{a_1}{T_1} = {a_E}{T_2} + b \qquad \qquad(10)</math></center> |

(3.427) | (3.427) | ||

where | where | ||

- | <center><math>{a_E} = \frac{{{k_1}{A_1}}}{{{{(\delta x)}_1}}} \qquad \qquad( )</math></center> | + | <center><math>{a_E} = \frac{{{k_1}{A_1}}}{{{{(\delta x)}_1}}} \qquad \qquad(11)</math></center> |

(3.428) | (3.428) | ||

- | <center><math>{a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} \qquad \qquad( )</math></center> | + | <center><math>{a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} \qquad \qquad(12)</math></center> |

(3.429) | (3.429) | ||

- | <center><math>b = {q''_B}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad( )</math></center> | + | <center><math>b = {q''_B}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(13)</math></center> |

(3.430) | (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 | ||

- | <center><math>{q''_B} = h({T_f} - {T_1}) \qquad \qquad( )</math></center> | + | <center><math>{q''_B} = h({T_f} - {T_1}) \qquad \qquad(14)</math></center> |

(3.431) | (3.431) | ||

- | It can be shown that eqs. (3.427) and (3.428) are still valid but eqs. ( | + | It can be shown that eqs. (3.427) and (3.428) are still valid but eqs. (12) and (13) should be modified to |

- | <center><math>{a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} + h{A_1} \qquad \qquad( )</math></center> | + | <center><math>{a_P} = {a_E} - {S_P}{A_1}{(\Delta x)_1} + h{A_1} \qquad \qquad(15)</math></center> |

(3.432) | (3.432) | ||

- | <center><math>b = h{T_f}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad( )</math></center> | + | <center><math>b = h{T_f}{A_1} + {S_C}{A_1}{(\Delta x)_1} \qquad \qquad(16)</math></center> |

(3.433) | (3.433) | ||

Line 105: | Line 105: | ||

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

- | <center><math>{a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}{T_{i - 1}} + {d_i},{\rm{ }}i = 1,2,...M \qquad \qquad( )</math></center> | + | <center><math>{a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}{T_{i - 1}} + {d_i},{\rm{ }}i = 1,2,...M \qquad \qquad(17)</math></center> |

(3.434) | (3.434) | ||

Line 114: | Line 114: | ||

To get the temperature, we hope that our temperature at the <math>{i_{th}}</math> grid is related to the neighbor on the right, i.e., | To get the temperature, we hope that our temperature at the <math>{i_{th}}</math> grid is related to the neighbor on the right, i.e., | ||

- | <center><math>{T_i} = {P_i}{T_{i + 1}} + {Q_i} \qquad \qquad( )</math></center> | + | <center><math>{T_i} = {P_i}{T_{i + 1}} + {Q_i} \qquad \qquad(18)</math></center> |

(3.435) | (3.435) | ||

- | For the <math>{(i - 1)^{{\rm{th}}}}</math> grid point, eq. ( | + | For the <math>{(i - 1)^{{\rm{th}}}}</math> grid point, eq. (18) becomes |

- | <center><math>{T_{i - 1}} = {P_{i - 1}}{T_i} + {Q_{i - 1}} \qquad \qquad( )</math></center> | + | <center><math>{T_{i - 1}} = {P_{i - 1}}{T_i} + {Q_{i - 1}} \qquad \qquad(19)</math></center> |

(3.436) | (3.436) | ||

- | Substituting eq. ( | + | Substituting eq. (19) into eq. (17), we get |

- | <center><math>{a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}({P_{i - 1}}{T_i} + {Q_{i - 1}}) + {b_i} \qquad \qquad( )</math></center> | + | <center><math>{a_i}{T_i} = {b_i}{T_{i + 1}} + {c_i}({P_{i - 1}}{T_i} + {Q_{i - 1}}) + {b_i} \qquad \qquad(20)</math></center> |

(3.437) | (3.437) | ||

which can be rearranged to | which can be rearranged to | ||

- | <center><math>{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( )</math></center> | + | <center><math>{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)</math></center> |

(3.438) | (3.438) | ||

- | Comparing eq. ( | + | Comparing eq. (18) and (21) yields |

- | <center><math>{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( )</math></center> | + | <center><math>{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)</math></center> |

(3.439) | (3.439) | ||

For the first grid point at the left <math>(i = 1)</math>, <math>{c_1} = 0</math> and eq. (3.439) becomes | For the first grid point at the left <math>(i = 1)</math>, <math>{c_1} = 0</math> and eq. (3.439) becomes | ||

- | <center><math>{P_1} = \frac{{{b_1}}}{{{a_1}}},{\rm{ }}{Q_1} = \frac{{{d_1}}}{{{a_1}}} \qquad \qquad( )</math></center> | + | <center><math>{P_1} = \frac{{{b_1}}}{{{a_1}}},{\rm{ }}{Q_1} = \frac{{{d_1}}}{{{a_1}}} \qquad \qquad(23)</math></center> |

(3.440) | (3.440) | ||

- | For the last point at the right <math>(i = M)</math>, <math>{b_M} = 0</math> and it follows that <math>{P_M} = 0</math>. Equation ( | + | For the last point at the right <math>(i = M)</math>, <math>{b_M} = 0</math> and it follows that <math>{P_M} = 0</math>. Equation (18) for <math>i = M</math> becomes |

- | <center><math>{T_N} = {Q_N} \qquad \qquad( )</math></center> | + | <center><math>{T_N} = {Q_N} \qquad \qquad(24)</math></center> |

(3.441) | (3.441) | ||

- | The solution of the algebraic equations can thus be completed in the following steps: (1) obtain <math>{P_1}</math> and <math>{Q_1}</math> from eq. (3.440), (2) calculate <math>{P_i}</math> and <math>{Q_i}</math> from eq. ( | + | The solution of the algebraic equations can thus be completed in the following steps: (1) obtain <math>{P_1}</math> and <math>{Q_1}</math> from eq. (3.440), (2) calculate <math>{P_i}</math> and <math>{Q_i}</math> 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. |

## Revision as of 00:28, 8 December 2009

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

(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*,*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

(3.419)

where *S*_{C} is a constant and SP is a coefficient.

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

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

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

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

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

(3.421)

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

(3.422)

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

(3.423)

where

(3.424)

(3.425)

(3.426)

## Contents |

### 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. 3.24). If the heat flux at the left boundary is known, the energy balance for the control volume at the left boundary is

which can be expressed as

(3.427)

where

(3.428)

(3.429)

(3.430)

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

(3.431)

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

(3.432)

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

### 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. (3.423), can be rewritten as

(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

*c*

_{1}= 0,

*b*

_{M}= 0

To get the temperature, we hope that our temperature at the *i*_{th} grid is related to the neighbor on the right, i.e.,

(3.435)

For the (*i* − 1)^{th} grid point, eq. (18) becomes

(3.436)

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

(3.437)

which can be rearranged to

(3.438)

Comparing eq. (18) and (21) yields

(3.439)

For the first grid point at the left (*i* = 1), *c*_{1} = 0 and eq. (3.439) becomes

(3.440)

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

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