# Numerical Solution of 1-D Steady Conduction

### From Thermal-FluidsPedia

Yuwen Zhang (Talk | contribs) |
|||

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

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

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

- | + | where <math>A(x)</math> is the area of heat conduction. Equation (1) 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. (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 | |

- | + | ||

- | + | ||

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

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

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

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

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

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

+ | |||

+ | 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)_w} = {k_w}{A_w}\frac{{{T_P} - {T_W}}}{{{{(\delta x)}_w}}}</math></center> | |

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

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

- | + | ||

- | + | ||

- | + | ||

- | + | ||

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

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

- | + | ||

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

- | + | ||

- | <center><math>{ | + | |

- | + | ||

- | where | + | where |

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

- | + | ||

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

- | + | ||

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

- | + | ==Boundary Conditions== | |

- | + | [[Image:Chapter3_(22).gif|thumb|400 px|alt=Control volumes at boundaries for Practice A|<center> (a) Left boundary (b) Right boundary </center> <br> '''Figure 2: Control volumes at boundaries for Practice A''']] | |

+ | [[Image:Chapter3_(23).gif|thumb|400 px|alt=Control volumes at boundaries for Practice A|<center> (a) Left boundary (b) Right boundary </center> <br> '''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 | |

- | + | ||

- | + | ||

- | + | ||

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

- | + | ||

- | + | ||

- | <center><math>{ | + | |

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | ||

- | + | which can be expressed as | |

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

- | + | ||

- | + | where | |

- | + | ||

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

- | + | ||

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

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

- | + | ||

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

- | + | ||

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

- | + | ||

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

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

- | + | ||

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

- | + | ||

- | + | 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 <math>{(\Delta x)_1}</math> 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 <math>{(\Delta x)_M}</math> 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 [[#References|(Patankar, 1980)]]. The discretized equations for every point, eq. (6), can be rewritten as | |

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

- | + | ||

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

- | + | ||

- | + | <center><math>{c_1} = 0,{\rm{ }}{b_M} = 0</math></center> | |

- | + | ||

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

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

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

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

- | + | ||

- | <center><math>\frac{{{ | + | |

- | + | ||

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

- | + | ||

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

- | + | ||

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

- | + | ||

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

- | + | ||

+ | 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. (23), (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. | ||

+ | ==References== | ||

- | + | Patankar, S.V., 1980, ''Numerical Heat Transfer and Fluid Flow'', Hemisphere, Washington, DC. | |

+ | Tao, W.Q., 2001, ''Numerical Heat Transfer'', 2<sup>nd</sup> Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese). | ||

+ | ==Further Reading== | ||

- | + | ==External Links== | |

- | + | ||

- | + | ||

- | + | ||

- | == | + |

## Current revision as of 08:42, 3 July 2010

## Contents |

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

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*) = *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

where *S*_{C} is a constant and *S*_{P} 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. 1), one obtains

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

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

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

where

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

which can be expressed as

where

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

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

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

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

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

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

which can be rearranged to

Comparing eq. (18) and (21) yields

For the first grid point at the left (*i* = 1), *c*_{1} = 0 and eq. (22) 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

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.

## References

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

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