# One-dimensional transient heat conduction in semi-infinite body

Figure 1: Heat conduction in a semi-infinite body

Let us consider heat conduction in a semi-infinite body (x > 0) with an initial temperature of Ti. At time t = 0, the surface temperature of the semi-infinite body is suddenly increased to a temperature T0. The temperature near the surface of the semi-infinite body will increase because of the surface temperature change, while the temperature far from the surface of the semi-infinite body is not affected and remains at the initial temperature Ti. The physical model of the problem is illustrated in Fig. 1, and the governing equation of the heat conduction problem and the corresponding initial and boundary conditions are:

$\frac{{{\partial ^2}T(x,t)}}{{\partial {x^2}}} = \frac{1}{\alpha }\frac{{\partial T(x,t)}}{{\partial t}}\quad \quad x > 0,\quad t > 0 \qquad \qquad(1)$
$T(x,t) = {T_0}\quad \quad x = 0,\quad t > 0 \qquad \qquad(2)$
$T(x,t) = {T_i}\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad(3)$
$T(x,t) = {T_i}\quad \quad x > 0,\quad t = 0 \qquad \qquad(4)$

which can be solved by using the method of separation of variables or integral approximate solution.

## Separation of Variables

Defining the following dimensionless variables

$\theta = \frac{{T - {T_0}}}{{{T_i} - {T_0}}},{\rm{ }}X = \frac{x}{L},{\rm{ Fo}} = \frac{{\alpha t}}{{{L^2}}} \qquad \qquad(5)$

where L is a characteristic length, eqs. (1) – (4) will be nondimensionalized as

$\frac{{{\partial ^2}\theta }}{{\partial {X^2}}} = \frac{{\partial \theta }}{{\partial {\rm{Fo}}}},{\rm{ }}X > 0,{\rm{ Fo}} > 0 \qquad \qquad(6)$
$\theta = 0,{\rm{ }}X = 0 \qquad \qquad(7)$
$\theta = 1,{\rm{ X}} \to \infty \qquad \qquad(8)$
$\theta = 1,{\rm{ }}X > 0,{\rm{ Fo}} = 0 \qquad \qquad(9)$

Assuming that the temperature can be expressed as

$\theta (X,{\rm{Fo}}) = \Theta (X)\Gamma ({\rm{Fo}}) \qquad \qquad(10)$

and substituting eq. (10) into eq. (6), one obtains

$\frac{{\Theta ''(X)}}{{\Theta (X)}} = \frac{{\Gamma '({\rm{Fo}})}}{{\Gamma ({\rm{Fo}})}} = - {\lambda ^2} \qquad \qquad(11)$

whose general solutions are

$\Theta = {C_1}\cos \lambda X + {C_2}\sin \lambda X \qquad \qquad(12)$
$\Gamma = {C_3}{e^{ - {\lambda ^2}{\rm{Fo}}}} \qquad \qquad(13)$

where C1,C2, and C3 are integral constants. Substituting eq. (10) into eq. (7), the following boundary condition of eq. (12) is obtained

$\Theta (0) = 0 \qquad \qquad(14)$

Substituting eq. (12) into eq. (14) yields C1 = 0 and eq. (12) becomes

$\Theta = {C_2}\sin \lambda X \qquad \qquad(15)$

The boundary condition at $X \to \infty$, eq. (8) cannot be used to get the eigen equation for λ and therefore, there will be no restriction on the value of λ for heat conduction in a semi-infinite body. Substituting eqs. (15) and (13) into eq. (10), the solution becomes

${\theta _\lambda } = C(\lambda )\sin \left( {\lambda X} \right){e^{ - {\lambda ^2}{\rm{Fo}}}} \qquad \qquad(16)$

where C = C2C3. The general solution for the problem can be obtained by using linear combination of eq. (16) for all possible λ, i.e.,

$\theta = \int_0^\infty {C(\lambda )\sin \left( {\lambda X} \right){e^{ - {\lambda ^2}{\rm{Fo}}}}d\lambda } \qquad \qquad(17)$

Substituting eq. (17) into eq. (9), one obtains

$1 = \int_{\lambda = 0}^\infty {C(\lambda )\sin \left( {\lambda X} \right)d\lambda }$

If we solve the problem by using Laplace transformation, we have

$1 = \int_{\lambda = 0}^\infty {\sin \left( {\lambda X} \right)\left[ {\frac{2}{\pi }\int_{X' = 0}^\infty {\sin \left( {\lambda X'} \right)dX'} } \right]d\lambda }$

Comparing the above two equations, an expression of C is obtained:

$C(\lambda ) = \frac{2}{\pi }\int_{X' = 0}^\infty {\sin \left( {\lambda X'} \right)dX'}$

The temperature distribution, eq. (17), becomes

$\theta = \frac{2}{\pi }\int_{X' = 0}^\infty {\int_{\lambda = 0}^\infty {{e^{ - {\lambda ^2}{\rm{Fo}}}}\sin \left( {\lambda X'} \right)dX'} \sin \left( {\lambda X} \right)d\lambda } dX'$

which can be rewritten as

$\theta = \frac{2}{\pi }\int_{X' = 0}^\infty {\int_{\lambda = 0}^\infty {{e^{ - {\lambda ^2}{\rm{Fo}}}}\left[ {\cos \lambda (X - X') - \cos \lambda (X + X')} \right]} d\lambda } dX'$

Evaluating the integral with respect toλ yields

$\int_{\lambda = 0}^\infty {{e^{ - {\lambda ^2}{\rm{Fo}}}}\cos \lambda (X - X')} d\lambda = \sqrt {\frac{\pi }{{4{\rm{Fo}}}}} \exp \left[ { - \frac{{{{(X - X')}^2}}}{{4{\rm{Fo}}}}} \right]$
$\int_{\lambda = 0}^\infty {{e^{ - {\lambda ^2}{\rm{Fo}}}}\cos \lambda (X + X')} d\lambda = \sqrt {\frac{\pi }{{4{\rm{Fo}}}}} \exp \left[ { - \frac{{{{(X + X')}^2}}}{{4{\rm{Fo}}}}} \right]$

thus

$\theta = \frac{1}{{\sqrt {4\pi {\rm{Fo}}} }}\left\{ {\int_{X' = 0}^\infty {\exp \left[ { - \frac{{{{(X - X')}^2}}}{{4{\rm{Fo}}}}} \right]dX'} - \int_{X' = 0}^\infty {\exp \left[ { - \frac{{{{(X + X')}^2}}}{{4{\rm{Fo}}}}} \right]dX'} } \right\} \qquad \qquad(18)$

Let us define a new variable

$\eta = \frac{{X - X'}}{{\sqrt {4{\rm{Fo}}} }}$

The first integral in eq. (18), becomes

$\int_{X' = 0}^\infty {\exp \left[ { - \frac{{{{(X - X')}^2}}}{{4{\rm{Fo}}}}} \right]dX'} = \sqrt {4{\rm{Fo}}} \int_{ - X/\sqrt {4{\rm{Fo}}} }^\infty {\exp ( - {\eta ^2})d\eta }$

Similarly, the second integral can be evaluated by following a similar procedure:

$\int_{X' = 0}^\infty {\exp \left[ { - \frac{{{{(X + X')}^2}}}{{4{\rm{Fo}}}}} \right]dX'} = \sqrt {4{\rm{Fo}}} \int_{X/\sqrt {4{\rm{Fo}}} }^\infty {\exp ( - {\eta ^2})d\eta }$

Substituting the above two equations into eq. (18), we have

$\theta = \frac{1}{{\sqrt \pi }}\left[ {\int_{ - X/\sqrt {4{\rm{Fo}}} }^\infty {\exp ( - {\eta ^2})d\eta } - \int_{X/\sqrt {4{\rm{Fo}}} }^\infty {\exp ( - {\eta ^2})d\eta } } \right] = \frac{2}{{\sqrt \pi }}\int_0^{X/\sqrt {4{\rm{Fo}}} } {\exp ( - {\eta ^2})d\eta }$

which can be rewritten as

$\theta = {\rm{erf}}\left( {\frac{X}{{\sqrt {4{\rm{Fo}}} }}} \right) \qquad \qquad(19)$

where erf in eq. (19) is the error function defined as:

${\rm{erf}}(z) = \frac{2}{{\sqrt \pi }}\int_0^z \,{e^{ - {z^2}}}dz \qquad \qquad(20)$

Equation (19) can also be rewritten as dimensional form:

$\frac{{T - {T_0}}}{{{T_i} - {T_0}}} = {\rm{erf}}\left( {\frac{x}{{\sqrt {4\alpha t} }}} \right) \qquad \qquad(21)$

The surface heat flux can be obtained by applying Fourier’s law

${q''_0}(t) = - k{\left. {\frac{{\partial T}}{{\partial x}}} \right|_{x = 0}} = \frac{{k({T_0} - {T_i})}}{{\sqrt {\pi \alpha t} }} \qquad \qquad(22)$

The solution of heat conduction in a semi-infinite body under the boundary conditions of the second and third kinds can also be obtained by using the method of separation of variables (Ozisik, 1993).

## Time-Dependent Boundary Condition

Our discussions thus far have been limited to the case that the boundary condition is not a function of time. Periodic boundary conditions can be encountered in various applications ranging from heat conduction in a building during day and night to emerging technologies such as pulsed laser processing of materials. Let us reconsider the problem described by eqs. (1) – (4) but replace eq. (2) by

$T = {T_i} + f(t) = {T_i} + A\cos (\omega t - \beta ),{\rm{ }}x = 0,\quad t > 0 \qquad \qquad(23)$

where A is the amplitude of oscillation, ω is the angular frequency, and β is the phase delay. The method of separation of variables will not work since eq. (23) is not homogeneous. Introducing excess temperature $\vartheta = T - {T_i}$, the governing equation and corresponding boundary and initial conditions become

$\frac{{{\partial ^2}\vartheta (x,t)}}{{\partial {x^2}}} = \frac{1}{\alpha }\frac{{\partial \vartheta (x,t)}}{{\partial t}}\quad \quad x > 0,\quad t > 0 \qquad \qquad(24)$
$\vartheta (x,t) = f(t) = A\cos (\omega t - \beta )\quad \quad x = 0,\quad t > 0 \qquad \qquad(25)$
$\vartheta (x,t) = 0\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad(26)$
$\vartheta (x,t) = 0\quad \quad x > 0,\quad t = 0 \qquad \qquad(27)$

Duhamel’s theorem can be used to handle the time-dependent boundary condition. Instead of solving eqs. (24) – (27) directly, we will start with a simpler auxiliary problem defined below:

$\frac{{{\partial ^2}\Phi (x,t)}}{{\partial {x^2}}} = \frac{1}{\alpha }\frac{{\partial \Phi (x,t)}}{{\partial t}}\quad \quad x > 0,\quad t > 0 \qquad \qquad(28)$
$\Phi (x,t) = f(\tau ) = A\cos (\omega \tau - \beta )\quad \quad x = 0,\quad t > 0 \qquad \qquad(29)$
$\Phi (x,t) = 0\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad(30)$
$\Phi (x,t) = 0\quad \quad x > 0,\quad t = 0 \qquad \qquad(31)$

where τ in eq. (29) is treated as a parameter, rather than time. Duhamel’s theorem (Ozisik, 1993) stated that the solution to the original problem is related to the solution of the auxiliary problem by

$\vartheta (x,t) = \frac{\partial }{{\partial t}}\int_{\tau = 0}^t {\Phi (x,t - \tau ,\tau )d\tau } \qquad \qquad(32)$

which can rewritten using Leibniz’s rule

$\vartheta (x,t) = \int_{\tau = 0}^t {\frac{\partial }{{\partial t}}\Phi (x,t - \tau ,\tau )d\tau } + \Phi {\left. {(x,t - \tau ,\tau )} \right|_{\tau = t}} \qquad \qquad(33)$

The second term on the right hand side is

$\Phi {\left. {(x,t - \tau ,\tau )} \right|_{\tau = t}} = \Phi (x,0,\tau ) = 0 \qquad \qquad(34)$

therefore, eq. (33) becomes

$\vartheta (x,t) = \int_{\tau = 0}^t {\frac{\partial }{{\partial t}}\Phi (x,t - \tau ,\tau )d\tau } \qquad \qquad(35)$

The solution of the auxiliary problem can be expressed as

$\Phi (x,t,\tau ) = f(\tau )\left[ {1 - {\rm{erf}}\left( {\frac{x}{{\sqrt {4\alpha t} }}} \right)} \right] = \frac{{2f(\tau )}}{{\sqrt \pi }}\int_{x/\sqrt {4\alpha t} }^\infty {{e^{ - {\eta ^2}}}d\eta } \qquad \qquad(36)$

The partial derivative appearing in eq. (35) can be evaluated as

$\frac{\partial }{{\partial t}}\Phi (x,t - \tau ,\tau ) = f(\tau )\frac{x}{{\sqrt {4\pi \alpha } {{(t - \tau )}^{3/2}}}}\exp \left[ { - \frac{{{x^2}}}{{4\alpha (t - \tau )}}} \right] \qquad \qquad(37)$

Substituting eq. (37) into eq. (35), the solution of the original problem becomes

$\vartheta (x,t) = \frac{x}{{\sqrt {4\pi \alpha } }}\int_{\tau = 0}^t {\frac{{f(\tau )}}{{{{(t - \tau )}^{3/2}}}}\exp \left[ { - \frac{{{x^2}}}{{4\alpha (t - \tau )}}} \right]d\tau } \qquad \qquad(38)$

Introducing a new independent variable

$\xi = \frac{x}{{\sqrt {4\alpha (t - \tau )} }}$

eq. (38) becomes

$\vartheta (x,t) = \frac{2}{{\sqrt \pi }}\int_{x/\sqrt {4\alpha t} }^\infty {f\left( {t - \frac{{{x^2}}}{{4\alpha {\xi ^2}}}} \right)\exp ( - {\xi ^2})d\xi } \qquad \qquad(39)$

For the periodic boundary condition specified in eq. (25), we have

$\vartheta (x,t) = \frac{{2A}}{{\sqrt \pi }}\int_{x/\sqrt {4\alpha t} }^\infty {\cos \left[ {\omega \left( {t - \frac{{{x^2}}}{{4\alpha {\xi ^2}}}} \right) - \beta } \right]\exp ( - {\xi ^2})d\xi } \qquad \qquad(40)$

which can be rewritten as

$\begin{array}{l} \vartheta (x,t) = \frac{{2A}}{{\sqrt \pi }}\int_0^\infty {\cos \left[ {\omega \left( {t - \frac{{{x^2}}}{{4\alpha {\xi ^2}}}} \right) - \beta } \right]\exp ( - {\xi ^2})d\xi } \\ {\rm{ }} - \frac{{2A}}{{\sqrt \pi }}\int_0^{x/\sqrt {4\alpha t} } {\cos \left[ {\omega \left( {t - \frac{{{x^2}}}{{4\alpha {\xi ^2}}}} \right) - \beta } \right]\exp ( - {\xi ^2})d\xi } \\ \end{array} \qquad \qquad(41)$

Evaluating the first integral on the right hand side of eq. (41) yields

$\begin{array}{l} \vartheta (x,t) = A\exp \left[ { - x{{\left( {\frac{\omega }{{2\alpha }}} \right)}^{1/2}}} \right]\cos \left[ {\omega t - x{{\left( {\frac{\omega }{{2\alpha }}} \right)}^{1/2}} - \beta } \right] \\ {\rm{ }} - \frac{{2A}}{{\sqrt \pi }}\int_0^{x/\sqrt {4\alpha t} } {\cos \left[ {\omega \left( {t - \frac{{{x^2}}}{{4\alpha {\xi ^2}}}} \right) - \beta } \right]\exp ( - {\xi ^2})d\xi } \\ \end{array} \qquad \qquad(42)$

It can be seen that as $t \to \infty$, the second term will become zero and the first term represents the steady oscillation.

${\vartheta _s}(x,t) = A\exp \left[ { - x{{\left( {\frac{\omega }{{2\alpha }}} \right)}^{1/2}}} \right]\cos \left[ {\omega t - x{{\left( {\frac{\omega }{{2\alpha }}} \right)}^{1/2}} - \beta } \right] \qquad \qquad(43)$

where $A\exp \left[ { - x{{\left( {\omega /(2\alpha )} \right)}^{1/2}}} \right]$ represents the amplitude of oscillation at point x, and $- x{\left( {\omega /(2\alpha )} \right)^{1/2}}$ in the cosine function represents the phase delay of oscillation at the point x relative to the oscillation of the surface temperature. It is apparent that the amplitude of the oscillation decreases as x increases, and the phase of oscillation delays with increasing x. The oscillation at x is not as strong as that at the surface and there is a delay from the time that the surface temperature changes to the time that the temperature at x responds to such change.

## Integral Approximate Solution

Figure 2: Heat conduction in a semi-infinite body with constant wall temperature

The early applications of integral approximate solutions to heat transfer problems included integral approximate solutions of boundary-layer momentum and energy equations. The method can also be used to solve linear or nonlinear transient conduction problems. The effect of sudden surface temperature change gradually propagates into the semi-infinite body. It is useful here to introduce a concept similar to the thermal boundary layer for convective heat transfer – thermal penetration depth. Assuming the thickness of the thermal penetration depth at time t is δ, the temperature of the semi-infinite body at x < δ will be affected but the temperature at x > δ will remain unchanged (see Fig. 2). It should be pointed out that the thermal penetration depth, δ, increases with increasing time.

According to the definition of the thermal penetration depth, the temperature at the thermal penetration depth should satisfy

$\frac{{\partial T(x,t)}}{{\partial x}} = 0\quad \quad x = \delta (t) \qquad \qquad(44)$
$T(x,t) = {T_i}\quad \quad x = \delta (t) \qquad \qquad(45)$

Integrating eq. (1) in the interval (0, δ), one obtains

$\mathop {\left. {\frac{{\partial T}}{{\partial x}}} \right|}{x = \delta (t)} - \mathop {\left. {\frac{{\partial T}}{{\partial x}}} \right|}{x = 0} = \frac{1}{\alpha }\int_0^{\delta (t)} \,\frac{{\partial T(x,t)}}{{\partial t}}dx \qquad \qquad(46)$

The right-hand side of eq. (46) can be rewritten using Leibnitz’s rule, i.e.,

$\mathop {\left. {\frac{{\partial T}}{{\partial x}}} \right|}{x = \delta (t)} - \mathop {\left. {\frac{{\partial T}}{{\partial x}}} \right|}{x = 0} = \frac{1}{\alpha }\left[ {\frac{d}{{dt}}\left( {\int_0^\delta \,Tdx} \right) - \mathop {\left. T \right|}{x = \delta } \frac{{d\delta }}{{dt}}} \right] \qquad \qquad(47)$

which represents the energy balance within the thermal penetration depth. Substituting eqs. (44) and (45) into eq. (47) yields

$\mathop {\left. { - \alpha \frac{{\partial T}}{{\partial x}}} \right|}{x = 0} = \frac{d}{{dt}}(\Theta - {T_i}\delta ) \qquad \qquad(48)$

where

$\Theta (t) = \int_0^{\delta (t)} \,T(x,t)dx \qquad \qquad(49)$

Equation (48) is the integral energy equation of the conduction problem, and this equation pertains to the entire thermal penetration depth. It follows that a temperature distribution that satisfies eq. (48) does not necessarily satisfy differential eq. (1), which describes the energy balance at any and all points in the domain of the problem.

In order to use the integral equation (48) to solve the conduction problem, the temperature profile in the thermal penetration depth is needed. The next step of the integral approximate solution is to assume a temperature distribution in the thermal penetration depth. The assumed temperature distribution can be any arbitrary function provided that the boundary conditions at x = 0 and x = δ are satisfied. Let us assume that the temperature distribution in the thermal penetration depth is a third-order polynomial function of x, i.e.,

$T(x,t) = {A_0} + {A_1}x + {A_2}{x^2} + {A_3}{x^3} \qquad \qquad(50)$

where A0, A1, A2, and A3 are four constants to be determined using the boundary conditions. Since there are only three boundary conditions available – eqs. (2), (44) and (45) – one more condition must be identified so that all four constants in eq. (50) can be determined. The surface temperature of the semi-infinite body, T0, is not a function of time t, so

$\frac{{\partial T(x,t)}}{{\partial t}} = 0\quad \quad x = 0 \qquad \qquad(51)$

Substituting eq. (1) into eq. (51) yields

$\frac{{{\partial ^2}T(x,t)}}{{\partial {x^2}}} = 0\quad \quad x = 0 \qquad \qquad(52)$

Substituting eq. (50) into eqs. (2), (44), (45) and (52) yields four equations for the constants in eq. (50). Solving for the four constants and substituting the results into eq. (50), the temperature distribution in the thermal penetration depth becomes

$\frac{{T(x,t) - {T_i}}}{{{T_0} - {T_i}}} = 1 - \frac{3}{2}\left( {\frac{x}{\delta }} \right) + \frac{1}{2}\mathop {\left( {\frac{x}{\delta }} \right)}^3 \qquad \qquad(53)$

where the thermal penetration depth, δ, is still unknown. Substituting eq. (53) into eq. (48), an ordinary differential equation for δis obtained:

$4\alpha = \delta \frac{{d\delta }}{{dt}}\quad \quad t > 0 \qquad \qquad(54)$

Since the thermal penetration depth equals zero at the beginning of the heat conduction, eq. (54) is subject to the following initial condition:

$\delta = 0\quad \quad t = 0 \qquad \qquad(55)$

The solution of eqs. (54) and (55) is

$\delta = \sqrt {8\alpha t} \qquad \qquad(56)$

which is consistent with the result of scale analysis, $\delta \sim \sqrt {\alpha t}$. The temperature distribution in the thermal penetration depth can be obtained by eq. (53), and the temperature in the semi-infinite body beyond the thermal penetration depth equals the initial temperature, Ti. The temperature distribution in the thermal penetration depth obtained here, as well as the thermal penetration depth thickness, depend on the assumed temperature distribution. The degree of the polynomial function for the temperature distribution in the thermal penetration depth should not be higher than four, or the result obtained will oscillate around the actual temperature profile, giving erroneous results. From the above analysis, we can summarize the procedure of the integral approximate solution as follows:

1. Obtain the integral equation of the problem by integrating the partial differential equation over the thermal penetration depth.

2. Assume an appropriate temperature distribution – usually a polynomial function – in the thermal penetration depth, and determine the unknown constants in the polynomial by using the boundary conditions at x=0 and x=δ.

3. Additional boundary conditions, if needed, can be obtained by further analysis of the boundary conditions and the conduction equation.

4. Obtain an ordinary differential equation of thermal penetration depth by substituting the temperature distribution into the integral equation. The thermal penetration depth thickness can be obtained by solving this ordinary differential equation.

The temperature distribution in the thermal penetration depth can be obtained by combining the thermal penetration depth thickness from step 3 into the temperature distribution from step 2.

## References

Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.

Ozisik, M.N., 1993, Heat Conduction, 2nd ed., Wiley-Interscience, New York.