Nonhomogeneous Problems

The transient one-dimensional conduction problems that we discussed so far are limited to the case that the problem is homogeneous and the method of separation of variables works. When the problem is not homogeneous due to a nonhomogeneous energy equation or boundary condition, the solution of a nonhomogeneous problem can be obtained by superposition of a particular solution of the nonhomogeneous problem and the general solution of the corresponding homogeneous problem.

Partial Solution

Let us consider a finite slab with thickness of L and a uniform initial temperature of Ti as shown in Fig. 3.17. At time t= 0, the temperature on the left side of the slab is suddenly increased to T0 while the temperature on the right side of the slab is maintained at Ti. Assuming that there is no internal heat generation in the slab and the thermophysical properties of the slab are constants, the energy equation is

$\frac{{{\partial ^2}T}}{{\partial {x^2}}} = \frac{1}{\alpha }\frac{{\partial T}}{{\partial t}},{\rm{ }}0 < x < L,{\rm{ }}t > 0 \qquad \qquad( )$

(3.268)

Insert Image Figure 3.17 Heat conduction under boundary condition of the first kind

subject to the following boundary and initial conditions

$T = {T_0},{\rm{ }}x = 0 \qquad \qquad( )$

(3.269)

$T = {T_i},{\rm{ }}x = L \qquad \qquad( )$

(3.270)

$T = {T_i},{\rm{ }}0 < x < L,{\rm{ }}t = 0 \qquad \qquad( )$

(3.271)

By defining the following dimensionless variables

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

(3.272)

eqs. (3.268) – (3.271) will be nondimensionalized as

$\frac{{{\partial ^2}\theta }}{{\partial {X^2}}} = \frac{{\partial \theta }}{{\partial {\rm{Fo}}}},{\rm{ }}0 < X < 1,{\rm{ Fo}} > 0 \qquad \qquad( )$

(3.273)

$\theta = 1,{\rm{ }}X = 0 \qquad \qquad( )$

(3.274)

$\theta = 0,{\rm{ X}} = 1 \qquad \qquad( )$

(3.275)

$\theta = 0,{\rm{ }}0 < X < 1,{\rm{ Fo}} = 0 \qquad \qquad( )$

(3.276)

It can be seen that the problem is still nonhomogeneous after nondimensionalization because eq. (3.274) is not homogeneous. When $t \to \infty$(i.e., $Fo \to \infty$), the temperature distribution can reach to steady state. If the steady state temperature is represented by θs, it must satisfy the following equations:

$\frac{{{\partial ^2}{\theta _s}}}{{\partial {X^2}}} = 0,{\rm{ }}0 < X < 1 \qquad \qquad( )$

(3.277)

${\theta _s} = 1,{\rm{ }}X = 0 \qquad \qquad( )$

(3.278)

${\theta _s} = 0,{\rm{ X}} = 1 \qquad \qquad( )$

(3.279)

which have the following solution:

${\theta _s} = 1 - X \qquad \qquad( )$

(3.280)

To obtain the generation of the problem described by eqs. (3.273) – (3.276), a method of partial solution (Myers, 1987) will be employed. In this methodology, it is assumed that the solution of a nonhomogeneous problem can be expressed as

$\theta (X,{\rm{Fo}}) = {\theta _s}(X) + {\theta _h}(X,{\rm{Fo}}) \qquad \qquad( )$

(3.281)

where θh represent the solution of a homogeneous problem. Substituting eqs. (3.273) – (3.276) and considering eqs. (3.277) – (3.279), we have

$\frac{{{\partial ^2}{\theta _h}}}{{\partial {X^2}}} = \frac{{\partial {\theta _h}}}{{\partial {\rm{Fo}}}},{\rm{ }}0 < X < 1,{\rm{ Fo}} > 0 \qquad \qquad( )$

(3.282)

${\theta _h} = 0,{\rm{ }}X = 0 \qquad \qquad( )$

(3.283)

${\theta _h} = 0,{\rm{ X}} = 1 \qquad \qquad( )$

(3.284)

${\theta _h} = X - 1,{\rm{ }}0 < X < 1,{\rm{ Fo}} = 0 \qquad \qquad( )$

(3.285)

which represent a new homogeneous problem. This problem can be solved using the method of separation of variables (see Problem 3.31) and the result is

${\theta _h} = - \frac{2}{\pi }\sum\limits_{n = 1}^\infty {\frac{{\sin (n\pi X)}}{n}{e^{ - {{(n\pi )}^2}{\rm{Fo}}}}} \qquad \qquad( )$

(3.286)

The solution of the nonhomogeneous problem thus becomes

$\theta = 1 - X - \frac{2}{\pi }\sum\limits_{n = 1}^\infty {\frac{{\sin (n\pi X)}}{n}{e^{ - {{(n\pi )}^2}{\rm{Fo}}}}} \qquad \qquad( )$

(3.287)

Variation of Parameter

The partial solution only works if the steady-state solution exists. If the steady-state solution does not exist, we can use the method of variation of parameters to solve the problem. Let us consider a finite slab with thickness of L and a uniform initial temperature of Ti. At time t= 0, the left side is subject to a constant heat flux while the right side of the slab is adiabatic (see Fig. 3.18). Assuming that there is no internal heat generation in the slab and the thermophysical properties of the slab are constants, the energy equation is

$\frac{{{\partial ^2}T}}{{\partial {x^2}}} = \frac{1}{\alpha }\frac{{\partial T}}{{\partial t}},{\rm{ }}0 < x < L,{\rm{ }}t > 0 \qquad \qquad( )$

(3.288)

subject to the following boundary and initial conditions

$- k\frac{{\partial T}}{{\partial x}} = {q''_0},{\rm{ }}x = 0 \qquad \qquad( )$

(3.289)

Insert Image Figure 3.18 Heat conduction under boundary condition of the second kind

$\frac{{\partial T}}{{\partial x}} = 0,{\rm{ }}x = L \qquad \qquad( )$

(3.290)

$T = {T_i},{\rm{ }}0 < x < L,{\rm{ }}t = 0 \qquad \qquad( )$

(3.291)

By defining the following dimensionless variables

$\theta = \frac{{T - {T_i}}}{{{{q''}_0}L/k}},{\rm{ }}X = \frac{x}{L},{\rm{ Fo}} = \frac{{\alpha t}}{{{L^2}}} \qquad \qquad( )$

(3.292)

eqs. (3.288) – (3.291) will be nondimensionalized as

$\frac{{{\partial ^2}\theta }}{{\partial {X^2}}} = \frac{{\partial \theta }}{{\partial {\rm{Fo}}}},{\rm{ }}0 < X < 1,{\rm{ Fo}} > 0 \qquad \qquad( )$

(3.293)

$\frac{{\partial \theta }}{{\partial X}} = - 1,{\rm{ }}X = 0 \qquad \qquad( )$

(3.294)

$\frac{{\partial \theta }}{{\partial X}} = 0,{\rm{ }}X = 1 \qquad \qquad( )$

(3.295)

$\theta = 0,{\rm{ }}0 < X < 1,{\rm{ Fo}} = 0 \qquad \qquad( )$

(3.296)

This nonhomogeneous problem does not have a steady-state solution, and therefore the partial solution cannot be applied. We will use the method of variation of parameters (Myers, 1987) to solve this problem. This method requires the following steps:

1. Set up a homogeneous problem by dropping the nonhomogeneous terms.

2. Solve the homogeneous problem to get eigenvalue λn and eigenfunctions Θn(X)

3. Assuming the solution of the original nonhomogeneous problem has the form of $\theta (X,{\rm{Fo}}) = \sum\limits_{n = 1}^n {{A_n}({\rm{Fo}}){\Theta _n}(X)}$

4. Solve for An(Fo) using orthogonal property of Θn

5. Obtain an ordinary differential equation (ODE) for An(Fo) and solve for An(Fo) from the ODE

6. Put together the final solution.

We will solve this nonhomogeneous problem by following the above procedure. The corresponding homogeneous problem is:

$\frac{{{\partial ^2}{\theta _h}}}{{\partial {X^2}}} = \frac{{\partial {\theta _h}}}{{\partial {\rm{Fo}}}},{\rm{ }}0 < X < 1,{\rm{ Fo}} > 0 \qquad \qquad( )$

(3.297)

$\frac{{\partial {\theta _h}}}{{\partial X}} = 0,{\rm{ }}X = 0 \qquad \qquad( )$

(3.298)

$\frac{{\partial {\theta _h}}}{{\partial X}} = 0,{\rm{ }}X = 1 \qquad \qquad( )$

(3.299)

${\theta _h} = 0,{\rm{ }}0 < X < 1,{\rm{ Fo}} = 0 \qquad \qquad( )$

(3.300)

Assuming the solution of the above homogeneous problem is

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

(3.301)

eq. (3.297) becomes

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

(3.302)

Since the objective here is to get the eigenvalue and eigen functions, we do not need to solve for Γ and only need to solve for Θ. The eigenvalue problem is

$\Theta '' + {\lambda ^2}\Theta = 0 \qquad \qquad( )$

(3.303)

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

(3.304)

$\Theta '(1) = 0 \qquad \qquad( )$

(3.305)

Solving eqs. (3.303) – (3.305) yields the following eigenvalues and eigen functions

${\lambda _n} = n\pi \qquad \qquad( )$

(3.306)

${\Theta _n}(X) = \cos (n\pi X),{\rm{ }}n = 0,1,2,... \qquad \qquad( )$

(3.307)

Now, let us assume that the solution of the original nonhomogeneous problem is

$\theta (X,{\rm{Fo}}) = \sum\limits_{n = 0}^\infty {{A_n}({\rm{Fo}})\cos \left( {n\pi X} \right)} \qquad \qquad( )$

(3.308)

Multiplying eq. (3.308) by $\cos \left( {m\pi X} \right)$ and integrating the resulting equation in the interval of (0, 1), one obtains

$\int_0^1 {\theta (X,{\rm{Fo}})\cos (m\pi X)dX} = \sum\limits_{n = 1}^\infty {{A_n}\int_0^1 {\cos (m\pi X)\cos (n\pi X)dX} } \qquad \qquad( )$

(3.309)

The integral on the right-hand side of eq. (3.309) can be evaluated as

$\int_0^1 {\cos (m\pi X)\cos (n\pi X)dX} = \left\{ {\begin{array}{*{20}{c}} {0,{\rm{ }}m \ne n{\rm{ }}} \\ \begin{array}{l} 1/2,{\rm{ }}m = n \ne 0 \\ 1,{\rm{ }}m = n = 0 \\ \end{array} \\ \end{array}} \right. \qquad \qquad( )$

(3.310)

thus, eq. (3.309) becomes

${A_{\rm{0}}}({\rm{Fo}}) = \int_0^1 {\theta (X,{\rm{Fo}})dX} ,{\rm{ }}m = 0 \qquad \qquad( )$

(3.311)

${A_{\rm{m}}}({\rm{Fo}}) = 2\int_0^1 {\theta (X,{\rm{Fo}})\cos (m\pi X)dX} ,{\rm{ }}m \ne 0 \qquad \qquad( )$

(3.312)

Differentiating eq. (3.311) with respect to Fo, one obtains:

$\frac{{d{A_{\rm{0}}}}}{{{\rm{dFo}}}} = \int_0^1 {\frac{{\partial \theta }}{{\partial {\rm{Fo}}}}dX} \qquad \qquad( )$

(3.313)

Substituting eq. (3.293) into eq. (3.313) and integrating with respect to X yields

$\frac{{d{A_{\rm{0}}}}}{{{\rm{dFo}}}} = \int_0^1 {\frac{{{\partial ^2}\theta }}{{\partial {X^2}}}dX = } {\left. {\frac{{\partial \theta }}{{\partial X}}} \right|_{X = 1}} - {\left. {\frac{{\partial \theta }}{{\partial X}}} \right|_{X = 0}} = 1 \qquad \qquad( )$

(3.314)

Integrating eq. (3.314) with respect to Fo, we have

${A_0}({\rm{Fo}}) = {\rm{Fo}} + {C_1} = \int_0^1 {\theta (X,{\rm{Fo}})dX} \qquad \qquad( )$

(3.315)

When Fo = 0, eq. (3.315) becomes

${A_0}(0) = {C_1} = \int_0^1 {\theta (X,0)dX} = 0 \qquad \qquad( )$

(3.316)

thus, we have

${A_0}({\rm{Fo}}) = {\rm{Fo}} \qquad \qquad( )$

(3.317)

Differentiating eq. (3.312) and considering eq. (3.293) yield

$\frac{{d{A_{\rm{m}}}}}{{d{\rm{Fo}}}} = 2\int_0^1 {\frac{{\partial \theta }}{{\partial {\rm{Fo}}}}\cos (m\pi X)dX} = 2\int_0^1 {\frac{{{\partial ^2}\theta }}{{\partial {X^2}}}\cos (m\pi X)dX} \qquad \qquad( )$

(3.318)

Using integration by parts twice, the following ODE is obtained:

$\frac{{d{A_{\rm{m}}}}}{{d{\rm{Fo}}}} = 2 - {(m\pi )^2}{A_m} \qquad \qquad( )$

(3.319)

Multiplying eq. (3.319) by an integrating factor ${e^{{{(m\pi )}^2}{\rm{Fo}}}}$, we have

$\frac{d}{{d{\rm{Fo}}}}\left[ {{A_{\rm{m}}}{e^{{{(m\pi )}^2}{\rm{Fo}}}}} \right] = 2{e^{{{(m\pi )}^2}{\rm{Fo}}}} \qquad \qquad( )$

(3.320)

which can be integrated to get

${A_m} = \frac{2}{{{{(m\pi )}^2}}} + {C_2}{e^{ - {{(m\pi )}^2}{\rm{Fo}}}} \qquad \qquad( )$

(3.321)

where C2 is an integral constant that needs to be determined by an initial condition. For Fo = 0, eq. (3.312) becomes

${A_{\rm{m}}}({\rm{0}}) = 2\int_0^1 {\theta (X,{\rm{0}})\cos (m\pi X)dX} = 2\int_0^1 {\cos (m\pi X)dX} = 0 \qquad \qquad( )$

(3.322)

Substituting eq. (3.321) into eq. (3.322), one obtains ${C_2} = - \frac{2}{{{{(m\pi )}^2}}}$

therefore, we have ${A_m} = \frac{2}{{{{(m\pi )}^2}}} - \frac{2}{{{{(m\pi )}^2}}}{e^{ - {{(m\pi )}^2}{\rm{Fo}}}}$

Changing m back to n for notation,

${A_n} = \frac{2}{{{{(n\pi )}^2}}} - \frac{2}{{{{(n\pi )}^2}}}{e^{ - {{(n\pi )}^2}{\rm{Fo}}}} \qquad \qquad( )$

(3.323)

Substituting eqs. (3.317) and (3.323) into eq. (3.308), the solution becomes

$\theta (X,{\rm{Fo}}) = {\rm{Fo + }}\frac{{\rm{2}}}{{{\pi ^{\rm{2}}}}}\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {n\pi X} \right)}}{{{n^2}}}} - \frac{{\rm{2}}}{{{\pi ^{\rm{2}}}}}\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {n\pi X} \right)}}{{{n^2}}}{e^{ - {{(n\pi )}^2}{\rm{Fo}}}}} \qquad \qquad( )$

(3.324)

When the time (Fourier number) becomes large, the last term on the right-hand side will become zero and the solution is represented by the first two terms only. To simplify eq. (3.324), let us assume the solution at large Fo can be expressed as

$\theta (X,{\rm{Fo}}) = {\rm{Fo + }}f(X) \qquad \qquad( )$

(3.325)

which is referred to as an asymptotic solution and it must satisfy eqs. (3.293) – (3.295). Substituting eq. (3.325) into eqs. (3.293) – (3.295), we have

$f''(X) = 1\qquad \qquad( )$

(3.326)

$f'(0) = - 1,{\rm{ }}f'(1) = 0\qquad \qquad( )$

(3.327)

Integrating eq. (3.326) and considering eq. (3.327), we obtain

$f(X) = \frac{{{X^2}}}{2} - X + C\qquad \qquad( )$

(3.328)

where C cannot be determined from eq. (3.327) because both boundary conditions are for the first order derivative. To determine C, we can expand f(X) defined in eq. (3.328) into cosine Fourier series, i.e.

$f(X) = \frac{{{X^2}}}{2} - X + C = {a_0} + \sum\limits_{n = 1}^\infty {{a_n}\cos (n\pi X)} \qquad \qquad( )$

(3.329)

After determining a0 and an, and considering that f(X) is identical to the second term on the right-hand side of eq. (3.324), we have

$\frac{{{X^2}}}{2} - X + \frac{1}{3} = \frac{{\rm{2}}}{{{\pi ^{\rm{2}}}}}\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {n\pi X} \right)}}{{{n^2}}}} \qquad \qquad( )$

(3.330)

Substituting eq. (3.330) into eq. (3.324), the final solution becomes

$\theta (X,{\rm{Fo}}) = {\rm{Fo + }}\frac{{{X^2}}}{2} - X + \frac{1}{3} - \frac{{\rm{2}}}{{{\pi ^{\rm{2}}}}}\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {n\pi X} \right)}}{{{n^2}}}{e^{ - {{(n\pi )}^2}{\rm{Fo}}}}} \qquad \qquad( )$

(3.331)

Transient 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. 3.19, 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( )$

(3.332)

Insert Image Figure 3.19 Heat conduction in a semi-infinite body

$T(x,t) = {T_0}\quad \quad x = 0,\quad t > 0 \qquad \qquad( )$

(3.333)

$T(x,t) = {T_i}\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad( )$

(3.334)

$T(x,t) = {T_i}\quad \quad x > 0,\quad t = 0 \qquad \qquad( )$

(3.335)

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

(3.336)

where L is a characteristic length, eqs. (3.332) – (3.335) 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( )$

(3.337)

$\theta = 0,{\rm{ }}X = 0 \qquad \qquad( )$

(3.338)

$\theta = 1,{\rm{ X}} \to \infty \qquad \qquad( )$

(3.339)

$\theta = 1,{\rm{ }}X > 0,{\rm{ Fo}} = 0 \qquad \qquad( )$

(3.340)

Assuming that the temperature can be expressed as

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

(3.341)

and substituting eq. (3.341) into eq. (3.337), one obtains

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

(3.342)

whose general solutions are

$\Theta = {C_1}\cos \lambda X + {C_2}\sin \lambda X \qquad \qquad( )$

(3.343)

$\Gamma = {C_3}{e^{ - {\lambda ^2}{\rm{Fo}}}} \qquad \qquad( )$

(3.344)

where C1,C2, and C3 are integral constants. Substituting eq. (3.341) into eq. (3.338), the following boundary condition of eq. (3.343) is obtained

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

(3.345)

Substituting eq. (3.343) into eq. (3.345) yields C1 = 0 and eq. (3.343) becomes

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

(3.346)

The boundary condition at $X \to \infty$, eq. (3.339) 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. (3.346) and (3.344) into eq. (3.341), the solution becomes

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

(3.347)

where C = C2C3. The general solution for the problem can be obtained by using linear combination of eq. (3.347) 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( )$

(3.348)

Substituting eq. (3.348) into eq. (3.340), 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. (3.348), 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( )$

(3.349) Let us define a new variable

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

The first integral in eq. (3.349), 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. (3.349), 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( )$

(3.350)

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

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

(3.351)

Equation (3.350) 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( )$

(3.352)

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

(3.353)

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. (3.332) – (3.335) but replace eq. (3.333) by

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

(3.354)

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. (3.354) 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( )$

(3.355)

$\vartheta (x,t) = f(t) = A\cos (\omega t - \beta )\quad \quad x = 0,\quad t > 0 \qquad \qquad( )$

(3.356)

$\vartheta (x,t) = 0\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad( )$

(3.357)

$\vartheta (x,t) = 0\quad \quad x > 0,\quad t = 0 \qquad \qquad( )$

(3.358)

Duhamel’s theorem can be used to handle the time-dependent boundary condition. Instead of solving eqs. (3.355) – (3.358) 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( )$

(3.359)

$\Phi (x,t) = f(\tau ) = A\cos (\omega \tau - \beta )\quad \quad x = 0,\quad t > 0 \qquad \qquad( )$

(3.360)

$\Phi (x,t) = 0\quad \quad x \to \infty ,\quad t > 0 \qquad \qquad( )$

(3.361)

$\Phi (x,t) = 0\quad \quad x > 0,\quad t = 0 \qquad \qquad( )$

(3.362)

where τ in eq. (3.360) 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( )$

(3.363)

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

(3.364)

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

(3.365)

therefore, eq. (3.364) becomes

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

(3.366)

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

(3.367)

The partial derivative appearing in eq. (3.366) 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( )$

(3.368)

Substituting eq. (3.368) into eq. (3.366), 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( )$

(3.369)

Introducing a new independent variable

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

eq. (3.369) becomes

<center>$\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( )$

(3.370)

For the periodic boundary condition specified in eq. (3.356), 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( )$

(3.371)

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

(3.372)

Evaluating the first integral on the right hand side of eq. (3.372) 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( )$

(3.373)

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

(3.374)

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

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. 3.20). 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( )$

(3.375)

$T(x,t) = {T_i}\quad \quad x = \delta (t) \qquad \qquad( )$

(3.376)

Integrating eq. (3.332) 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( )$

(3.377)

Insert Image Figure 3.20 Heat conduction in a semi-infinite body with constant wall temperature.

The right-hand side of eq. (3.377) 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( )$

(3.378)

which represents the energy balance within the thermal penetration depth. Substituting eqs. (3.375) and (3.376) into eq. (3.378) yields

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

(3.379)

where

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

(3.380)

Equation (3.379) 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. (3.379) does not necessarily satisfy differential eq. (3.332), which describes the energy balance at any and all points in the domain of the problem.

In order to use the integral equation (3.379) 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( )$

(3.381)

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. (3.333), (3.375) and (3.376) – one more condition must be identified so that all four constants in eq. (3.381) 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( )$

(3.382)

Substituting eq. (3.332) into eq. (3.382) yields

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

(3.383)

Substituting eq. (3.381) into eqs. (3.333), (3.375), (3.376) and (3.383) yields four equations for the constants in eq. (3.381). Solving for the four constants and substituting the results into eq. (3.381), 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( )$

(3.384)

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

$4\alpha = \delta \frac{{d\delta }}{{dt}}\quad \quad t > 0$

(3.385)

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

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

(3.386)

The solution of eqs. (3.385) and (3.386) is

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

(3.387)

which is consistent with the result of scale analysis, $\delta \sim \sqrt {\alpha t}$ [see eq. (1.186)]. The temperature distribution in the thermal penetration depth can be obtained by eq. (3.384), 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.