# Approximate Solution of the Nonlinear Heat Conduction Equation in a Semi-Infinite Domain

Authors: Jun Yu; Yi Yang; Antonio Campo
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

We use an approximation method to study the solution to a nonlinear
heat conduction equation in a semi-infinite domain. By expanding an energy density
function (defined as the internal energy per unit volume) as a Taylor polynomial in
a spatial domain, we reduce the partial differential equation to a set of first-order
ordinary differential equations in time. We describe a systematic approach to derive
approximate solutions using Taylor polynomials of a different degree. For a special case,
we derive an analytical solution and compare it with the result of a self-similar analysis.
A comparison with the numerically integrated results demonstrates good accuracy of our
approximate solutions. We also show that our approximation method can be applied
to cases where boundary energy density and the corresponding effective conductivity
are more general than those that are suitable for the self-similar method. Propagation
of nonlinear heat waves is studied for different boundary energy density and the
conductivity functions.


---

## Body

## 1. Introduction

The nonlinear heat equation, as given in (1.1), has applications in various branches of science and engineering, including thermal processing of materials [1], liquid movement in porous media [2], and radiation heat wave [3]. In particular, radiation heat wave plays an important role in indirect drive inertial confinement fusion (ICF) [4] and has a potential application in measurement of opacity [5] and equation of state [6], of hot and dense materials. In this paper, we focus on analysis of nonlinear heat waves governed by a nonlinear heat equation in a semi-infinite domain, with zero initial condition and a prescribed boundary condition at the origin. In principle, a semi-infinite solid extends to infinity in all but one direction. As a result, a single identifiable surface characterizes a semi-infinite solid as illustrated in Figure 1. The semi-infinite solid provides a useful idealization for two types of practical problems in transient heat conduction. For instance, it may be used to determine transient heat transfer near the surface of the earth or to approximate the transient response of a finite solid, such as a thick slab. For the slab application, the approximation would be reasonable for the early portion of the transient near the surface. In many applications, if a thermal change is imposed at the surface, a one-dimensional temperature wave will be propagated by heat conduction within the semi-infinite solid. Specifically, we study the 1D nonlinear heat conduction equation:
(1.1)ρ(T)cp(T)∂T(x,t)∂t=∂∂x(k(T)∂T∂x),
where ρ is density, cp is the specific heat at constant pressure, and k is thermal conductivity. In this paper, we assume all of them depend only on temperature T(x,t). Following Zel'dovich and Raizer [3], we introduce internal energy per unit volume (1.2)E(T)=∫0Tρ(θ)cp(θ)dθ,
then the temperature can be solved as a function of the energy density E(x,t)(1.3)T(x,t)=T(E),
and the nonlinear heat conduction equation can be rewritten in terms of E,(1.4)∂E(x,t)∂t=∂∂x(ke(E)∂E∂x),
where
(1.5)ke(E)=k(E)dTdE
is considered as effective heat conductivity in this paper. By virtue of (1.2) and (1.3), we have a one-to-one function correspondence between the temperature T and the energy density E.Figure 1
Schematic of a semi-infinite body.With the heat fluxJ(x,t) defined as (1.6)J(x,t)=-ke(E)∂E∂x,
equation (1.4) can be equivalently written as (1.7)∂E∂t=-∂J(x,t)∂x.
We choose the initial and boundary conditions as (1.8)T(x,0)=0forx>0,T(0,t)=T0(t)
or correspondingly(1.9)E(x,0)=0,forx>0,(1.10)E(0,t)=E0(t).
Also, we assume that the effective heat conductivity ke(E) tends to zero with the energy density E; therefore, the thermal wave has a finite extension and a definite front [3]. With x=h(t) as the front coordinate of thermal wave at time t, we have (1.11)E(x,t)|x≥h(t)=ke(E)|x≥h(t)=J(x,t)|x≥h(t)=0.Analytical solutions can be obtained for materials that possess constant thermophysical properties. However, when the thermophysical properties are affected by temperature, finite-difference techniques or elaborate analytical procedures need to be employed. For example, Barbaro et al. [7] applied the Kirchhoff transform to the enthalpy formulation of the heat conduction equation to obtain approximate solutions for temperature-dependent thermal properties. The solutions compared favorably with those produced by numerical techniques. Singh et al. [8] used a meshless element-free Galerkin method to obtain numerical solutions of a semi-infinite solid with temperature-dependent thermal conductivity. A quasilinearization scheme is adopted to avoid iterations, and for the time integration the backward-difference procedure was employed.Further, a variety of analytic techniques [9] have been applied to the nonlinear heat conduction problems with temperature-dependent thermal conductivity. One of main approaches is the similarity analysis. Marshak [10] was the first to obtain a self-similar solution to the nonlinear radiation heat conduction equation. Zel'dovich and Razier [3] gave a thorough and clear account of the propagation of thermal wave using self-similar method. Recently, Sigel et al. [11–14] extended the self-similar approach and applied it to experimental study of the ablative and nonablative radiation thermal wave. However, all self-similar schemes impose restriction on the choices of boundary temperature and heat conductivity. Specifically, in order to obtain a self-similar solution to the parabolic differential equation, boundary temperature was assumed to be either an exponential [10] or a power function [3] of time and conductivity was assumed to be a power function of temperature. This limited the application of self-similar method to more general and realistic cases. Another approach, initiated by Parlange, is based on the observation that even though the conductivity may be a complicated function, its integral is far easier to handle. By assuming such a chosen integral a quadratic polynomial in spatial coordinate, Parlange et al. [2] obtained approximate analytical solution of the nonlinear diffusion equation for arbitrary boundary conditions.In this paper, we propose to analyze an integral with respect to the spatial variable involving energy densityE directly. First an integral form of the nonlinear diffusion equation is derived. Due to the strong nonlinear nature of the problem considered here, the profile of the energy density has a sharp wave front. This allows us to approximate the integrals in the new formulation by expanding the energy density as Taylor polynomial of spatial variable from the boundary point to the wave front. Then we derive a system of ordinary differential equations governing the evolution of wave front and energy density profile approximately. These will be illustrated in Section 2. Approximation using a linear Taylor expansion will be described in Section 2.1. One advantage of our approach is that it can be improved by simply expanding the energy density as a higher-order polynomial, as discussed in Section 2.2. In addition, our method has no requirements on the form of boundary temperature and conductivity, or correspondingly, no requirements on the form of boundary energy density and effective conductivity. (As mentioned earlier, all self-similar methods do.) In general, our method is suitable for many different function forms describing the dependency of heat conductivity on energy density, as long as it equals zero when energy density is zero (e.g., ke=k0(E+bE2+⋯)). Therefore, our method is more general as compared with the self-similar method. To illustrate the idea, several examples are studied in Section 3, demonstrating good agreement between numerical and our approximate solutions of the heat conduction equation. Finally, we make concluding remarks in Section 4.


## 2. Approximation Method

We start with a derivation of an integral formulation of our nonlinear heat wave problem described by (1.4) and (1.9)–(1.11). Integrating (1.7) from 0 to x, one finds (2.1)∂∂t[∫0xE(ζ,t)dζ]=J0(t)-J(x,t),
where (2.2)J0(t)=J(0,t)=-ke(E0)E1(t),(2.3)E1(t)=∂E(x,t)∂x|x=0.
We define (2.4)S(E)=∫0Eke(ζ)dζ,
and from (1.6) we may write the heat flux as (2.5)J(x,t)=-∂S(E)∂x.
Also it can be verified that (2.6)S(x,t)|x=h(t)=S(E)|E=0=0.
Substituting (2.5) into (2.1) and integrating (2.1) from 0 to x, we obtain (2.7)∂∂t[∫0xdη∫0ηE(ζ,t)dζ]=S(x,t)-S(E0)+J0(t)x.
Noting that(2.8)∂∂t[∫0xdη∫0ηE(ζ,t)dζ]=∂∂t[∫0x(x-ζ)E(ζ,t)dζ]=x∂∂t[∫0xE(ζ,t)dζ]-∂∂t[∫0xζE(ζ,t)dζ],
then substituting (2.8) and (2.1) into (2.7) results in (2.9)∂∂t[∫0xζE(ζ,t)dζ]=S(E0)-S(x,t)-J(x,t)x.
Notice that (2.9) is an equivalent integral form of the nonlinear heat conduction equation (1.4).Substitutingx=h(t) and the boundary condition (1.11) and (2.6) into (2.9) and (2.1) gives (2.10)ddt[∫0hζE(ζ,t)dζ]=S(E0),(2.11)ddt[∫0hE(ζ,t)dζ]=J0(t).
It should be noted that up to this point no approximation has been made. We now show that we can systematically approximate the integrals in (2.10) and (2.11) using Taylor polynomials of a different degree for E(x,t) about the boundary point x=0. This approach leads to an approximation method for solving the nonlinear heat conduction equation. What follows next is a derivation of the linear approximation, using the Taylor polynomial of degree 1 for E(x,t). Higher-order approximation will be described in Section 2.2.


### 2.1. Linear Approximation

We first study the integral term in (2.11), which represents the area under the energy density curve. It is well known that when the heat conductivity k(T) is a strong nonlinear function of temperature T, or correspondingly, when the effective heat conductivity ke(E) is a strong nonlinear function of energy density E, the solution to the nonlinear heat conduction equation exhibits a sharp front. So the integral represented as an area under the energy density function curve can be approximated by the area under the linear expansion of the energy density from the boundary point (x=0) to the wave front (x=h(t)). A similar analysis may be applied to the integral in (2.10). Therefore, to approximate the integral terms in (2.10) and (2.11), we expand the energy density E(x,t) as a linear function of x:
(2.12)E(x,t)≈E0(t)+E1(t)x,0≤x≤h(t),
where h(t) is the wave front, and E0(t) and E1(t) are defined in (1.10) and (2.3). Notice that E0(t) is known from the boundary condition (1.10) so that we have two unknown functions of time, E1(t) and h(t) only. In order to obtain two equations to close the system, we substitute (2.12) and (2.2) into (2.10) and (2.11). This yields (2.13)ddt(E02h2+13E1h3)=S(E0),(2.14)ddt(E0h+E12h2)=-ke(E0)E1.
Integrating (2.13) and using the initial condition (2.15)h(t)|t=0=E1h3|t=0=0,
we may solve for h2 in the form (2.16)h2=G(t)(E02+E1h3)-1,
where (2.17)G(t)=∫0tS(E0(t1))dt1.
From (2.14), we have(2.18)hddt[h(E0+E1h2)]=-ke(E0)E1h
or(2.19)h2ddt(E0+E1h2)+12(E0+E1h2)d(h2)dt=-ke(E0)E1h.
Using (2.16), we get
(2.20)G(t)(E02+E1h3)-1ddt(E0+E1h2)+12(E0+E1h2)ddt[G(t)(E02+E1h3)-1]=-ke(E0)E1h
or(2.21)G(t)(E02+E1h3)-1[12-16(E0+E1h2)(E02+E1h3)-1]ddt(E1h)=-ke(E0)E1h-12(E0+E1h2)(E02+E1h3)-1[dG(t)dt-G(t)2(E02+E1h3)-1dE0dt]-G(t)(E02+E1h3)-1dE0dt.
Substituting (2.17) into (2.21) yields(2.22)[12-16(E0+E1h2)(E02+E1h3)-1]ddt(E1h)=-ke(E0)E1hG(t)(E02+E1h3)-12G(t)(E0+E1h2)[S(E0)-G(t)2(E02+E1h3)-1dE0dt]-dE0dt
or(2.23)ddt(E1h)=[12-16(E0+E1h2)(E02+E1h3)-1]-1×{-ke(E0)E1hG(t)(E02+E1h3)-12G(t)(E0+E1h2)S(E0)+[14(E0+E1h2)(E02+E1h3)-1-1]dE0dt}.Equation (2.23) is a nonlinear ordinary differential equation for E1h. It is easy to understand that the approximate expansion (2.12) is suitable only for the case of h(t)>0. From (2.15), (2.23) is valid only for t>0. For the initial time, we may simply assume that (2.24)ddt(E1h)|t=0=E1h|t=0=0.WithE0 given in the boundary condition (1.10), E1h can be solved from (2.23), (2.24), and (2.17). Then the heat wave front h will be available from (2.16) and E1 from E1h and h. Finally, the approximate values of the energy density E can be estimated using (2.12) and (1.11).


### 2.2. Higher-Order Approximation

To improve our approximation to the solution of the nonlinear heat conduction equation, we extend our linear approximation method, as described in Section2.1, to higher order. Following a procedure similar to that in Section 2.1, it is easy to see that the integral term in (2.11) can be calculated more accurately by expanding the energy density as a quadratic polynomial of x and we can obtain more accurate wave fronts and energy density profiles. In general, we may expand the energy density E(x,t) as a Taylor polynomial of degree up to M (M≥2) in x: (2.25)E(x,t)=∑j=0MEj(t)xj,0≤x≤h(t),(2.26)Ej(t)=1j!∂jE(x,t)∂xj|x=0.
In this case, we have M+1 unknowns, E1(t),E2(t),…,EM(t), and h(t). (E0(t) is known from the boundary condition (1.10).) Substituting (2.25) and (2.2) into (2.10) and (2.11) provides two equations:(2.27)ddt[∑j=0MEj(t)j+2hj+2]=S(E0),(2.28)ddt[∑j=0MEj(t)j+1hj+1]=-ke(E0)E1(t).
In order to close the system, we need (M-1) more equations. Using (2.4), the heat conduction equation (1.4) can be written as (2.29)∂E(x,t)∂t=∂2S(x,t)∂x2.
For positive integer l=2,…,M, we study (2.30)ddt[∫0hζlE(ζ,t)dζ].
We show in Appendix A that, together with (2.27) and (2.28), this approach leads to the following M+1 equations: (2.31)h2=G(t)(∑j=0MEjhjj+2)-1,(2.32)∑j=1M[1j+1-12(j+2)(∑i=0MEihii+1)(∑i=0MEihii+2)-1]h2d(Ejhj)dt=-S1(t)h+[14(∑i=0MEihii+1)(∑i=0MEihii+2)-1-1]h2dE0dt-12(∑i=0MEihii+1)(∑i=0MEihii+2)-1S(E0),(2.33)∑j=1M[1j+l+1-l+12(j+2)(∑i=0MEihii+l+1)(∑i=0MEihii+2)-1]h2d(Ejhj)dt=[l+14(∑i=0MEihii+l+1)(∑i=0MEihii+2)-1-1l+1]h2dE0dt-l+12(∑i=0MEihii+l+1)(∑i=0MEihii+2)-1S(E0)+l(l-1)∑j=0M+1Sj(t)hjj+l-1,l=2,…,M,
where
(2.34)Sj(t)=1j!∂jS(x,t)∂xj|x=0=Sj(E0,E1⋯Ej),0≤j≤M,
and SM+1(t) is determined by
(2.35)∑j=0M+1Sj(t)hj=0.Equations (2.31), (2.32), and (2.33) form a close set for the M+1 unknowns E1(t)⋯EM(t), and h(t). Once E1(t)⋯EM(t), and h(t) are available, the energy density profile can be determined from (2.25) and (1.11).For the case of a quadratic approximation,M=2, as shown in Appendix B, we may get(2.36)h2=G(t)(E02+E1h3+E2h2j+2)-1,(2.37)ddt(E1h)=1h2C0B2-C2B0A0B2-A2B0=1G(t)(∑j=02Ejhjj+2)C0B2-C2B0A0B2-A2B0,(2.38)ddt(E2h2)=1h2A0C2-A2C0A0B2-A2B0=1G(t)(∑j=02Ejhjj+2)A0C2-A2C0A0B2-A2B0,
where(2.39)A0=12-16(∑i=02Eihii+1)(∑i=02Eihii+2)-1,B0=13-18(∑i=02Eihii+1)(∑i=02Eihii+2)-1,(2.40)C0=-S1h+[14(∑i=02Eihii+1)(∑i=02Eihii+2)-1-1]h2dE0dt-12(∑i=02Eihii+1)(∑i=02Eihii+2)-1S(E0),(2.41)A2=14-12(∑i=02Eihii+3)(∑i=02Eihii+2)-1,B2=15-38(∑i=02Eihii+3)(∑i=02Eihii+2)-1,(2.42)C2=[34(∑i=02Eihii+3)(∑i=02Eihii+2)-1-13]h2dE0dt-32(∑i=02Eihii+3)(∑i=02Eihii+2)-1S(E0)+2∑j=03Sj(t)hjj+1,
and the S1h, S2h2, and S3h3 terms that appeared in (2.40) and (2.42) are given by (2.43)S1h=h∂S(x,t)∂x|x=0=ke(E0)E1h,S2h2=h22∂2S(x,t)∂x2|x=0=12ke′(E0)(E1h)2+ke(E0)E2h2,S3h3=-[S2h2+S1h+S(E0)].
Here (2.36), (2.37), and (2.38), forming a close set of equations for h2, E1, and E2, constitutes the formulas of the quadratic approximation.


## 2.1. Linear Approximation

We first study the integral term in (2.11), which represents the area under the energy density curve. It is well known that when the heat conductivity k(T) is a strong nonlinear function of temperature T, or correspondingly, when the effective heat conductivity ke(E) is a strong nonlinear function of energy density E, the solution to the nonlinear heat conduction equation exhibits a sharp front. So the integral represented as an area under the energy density function curve can be approximated by the area under the linear expansion of the energy density from the boundary point (x=0) to the wave front (x=h(t)). A similar analysis may be applied to the integral in (2.10). Therefore, to approximate the integral terms in (2.10) and (2.11), we expand the energy density E(x,t) as a linear function of x:
(2.12)E(x,t)≈E0(t)+E1(t)x,0≤x≤h(t),
where h(t) is the wave front, and E0(t) and E1(t) are defined in (1.10) and (2.3). Notice that E0(t) is known from the boundary condition (1.10) so that we have two unknown functions of time, E1(t) and h(t) only. In order to obtain two equations to close the system, we substitute (2.12) and (2.2) into (2.10) and (2.11). This yields (2.13)ddt(E02h2+13E1h3)=S(E0),(2.14)ddt(E0h+E12h2)=-ke(E0)E1.
Integrating (2.13) and using the initial condition (2.15)h(t)|t=0=E1h3|t=0=0,
we may solve for h2 in the form (2.16)h2=G(t)(E02+E1h3)-1,
where (2.17)G(t)=∫0tS(E0(t1))dt1.
From (2.14), we have(2.18)hddt[h(E0+E1h2)]=-ke(E0)E1h
or(2.19)h2ddt(E0+E1h2)+12(E0+E1h2)d(h2)dt=-ke(E0)E1h.
Using (2.16), we get
(2.20)G(t)(E02+E1h3)-1ddt(E0+E1h2)+12(E0+E1h2)ddt[G(t)(E02+E1h3)-1]=-ke(E0)E1h
or(2.21)G(t)(E02+E1h3)-1[12-16(E0+E1h2)(E02+E1h3)-1]ddt(E1h)=-ke(E0)E1h-12(E0+E1h2)(E02+E1h3)-1[dG(t)dt-G(t)2(E02+E1h3)-1dE0dt]-G(t)(E02+E1h3)-1dE0dt.
Substituting (2.17) into (2.21) yields(2.22)[12-16(E0+E1h2)(E02+E1h3)-1]ddt(E1h)=-ke(E0)E1hG(t)(E02+E1h3)-12G(t)(E0+E1h2)[S(E0)-G(t)2(E02+E1h3)-1dE0dt]-dE0dt
or(2.23)ddt(E1h)=[12-16(E0+E1h2)(E02+E1h3)-1]-1×{-ke(E0)E1hG(t)(E02+E1h3)-12G(t)(E0+E1h2)S(E0)+[14(E0+E1h2)(E02+E1h3)-1-1]dE0dt}.Equation (2.23) is a nonlinear ordinary differential equation for E1h. It is easy to understand that the approximate expansion (2.12) is suitable only for the case of h(t)>0. From (2.15), (2.23) is valid only for t>0. For the initial time, we may simply assume that (2.24)ddt(E1h)|t=0=E1h|t=0=0.WithE0 given in the boundary condition (1.10), E1h can be solved from (2.23), (2.24), and (2.17). Then the heat wave front h will be available from (2.16) and E1 from E1h and h. Finally, the approximate values of the energy density E can be estimated using (2.12) and (1.11).


## 2.2. Higher-Order Approximation

To improve our approximation to the solution of the nonlinear heat conduction equation, we extend our linear approximation method, as described in Section2.1, to higher order. Following a procedure similar to that in Section 2.1, it is easy to see that the integral term in (2.11) can be calculated more accurately by expanding the energy density as a quadratic polynomial of x and we can obtain more accurate wave fronts and energy density profiles. In general, we may expand the energy density E(x,t) as a Taylor polynomial of degree up to M (M≥2) in x: (2.25)E(x,t)=∑j=0MEj(t)xj,0≤x≤h(t),(2.26)Ej(t)=1j!∂jE(x,t)∂xj|x=0.
In this case, we have M+1 unknowns, E1(t),E2(t),…,EM(t), and h(t). (E0(t) is known from the boundary condition (1.10).) Substituting (2.25) and (2.2) into (2.10) and (2.11) provides two equations:(2.27)ddt[∑j=0MEj(t)j+2hj+2]=S(E0),(2.28)ddt[∑j=0MEj(t)j+1hj+1]=-ke(E0)E1(t).
In order to close the system, we need (M-1) more equations. Using (2.4), the heat conduction equation (1.4) can be written as (2.29)∂E(x,t)∂t=∂2S(x,t)∂x2.
For positive integer l=2,…,M, we study (2.30)ddt[∫0hζlE(ζ,t)dζ].
We show in Appendix A that, together with (2.27) and (2.28), this approach leads to the following M+1 equations: (2.31)h2=G(t)(∑j=0MEjhjj+2)-1,(2.32)∑j=1M[1j+1-12(j+2)(∑i=0MEihii+1)(∑i=0MEihii+2)-1]h2d(Ejhj)dt=-S1(t)h+[14(∑i=0MEihii+1)(∑i=0MEihii+2)-1-1]h2dE0dt-12(∑i=0MEihii+1)(∑i=0MEihii+2)-1S(E0),(2.33)∑j=1M[1j+l+1-l+12(j+2)(∑i=0MEihii+l+1)(∑i=0MEihii+2)-1]h2d(Ejhj)dt=[l+14(∑i=0MEihii+l+1)(∑i=0MEihii+2)-1-1l+1]h2dE0dt-l+12(∑i=0MEihii+l+1)(∑i=0MEihii+2)-1S(E0)+l(l-1)∑j=0M+1Sj(t)hjj+l-1,l=2,…,M,
where
(2.34)Sj(t)=1j!∂jS(x,t)∂xj|x=0=Sj(E0,E1⋯Ej),0≤j≤M,
and SM+1(t) is determined by
(2.35)∑j=0M+1Sj(t)hj=0.Equations (2.31), (2.32), and (2.33) form a close set for the M+1 unknowns E1(t)⋯EM(t), and h(t). Once E1(t)⋯EM(t), and h(t) are available, the energy density profile can be determined from (2.25) and (1.11).For the case of a quadratic approximation,M=2, as shown in Appendix B, we may get(2.36)h2=G(t)(E02+E1h3+E2h2j+2)-1,(2.37)ddt(E1h)=1h2C0B2-C2B0A0B2-A2B0=1G(t)(∑j=02Ejhjj+2)C0B2-C2B0A0B2-A2B0,(2.38)ddt(E2h2)=1h2A0C2-A2C0A0B2-A2B0=1G(t)(∑j=02Ejhjj+2)A0C2-A2C0A0B2-A2B0,
where(2.39)A0=12-16(∑i=02Eihii+1)(∑i=02Eihii+2)-1,B0=13-18(∑i=02Eihii+1)(∑i=02Eihii+2)-1,(2.40)C0=-S1h+[14(∑i=02Eihii+1)(∑i=02Eihii+2)-1-1]h2dE0dt-12(∑i=02Eihii+1)(∑i=02Eihii+2)-1S(E0),(2.41)A2=14-12(∑i=02Eihii+3)(∑i=02Eihii+2)-1,B2=15-38(∑i=02Eihii+3)(∑i=02Eihii+2)-1,(2.42)C2=[34(∑i=02Eihii+3)(∑i=02Eihii+2)-1-13]h2dE0dt-32(∑i=02Eihii+3)(∑i=02Eihii+2)-1S(E0)+2∑j=03Sj(t)hjj+1,
and the S1h, S2h2, and S3h3 terms that appeared in (2.40) and (2.42) are given by (2.43)S1h=h∂S(x,t)∂x|x=0=ke(E0)E1h,S2h2=h22∂2S(x,t)∂x2|x=0=12ke′(E0)(E1h)2+ke(E0)E2h2,S3h3=-[S2h2+S1h+S(E0)].
Here (2.36), (2.37), and (2.38), forming a close set of equations for h2, E1, and E2, constitutes the formulas of the quadratic approximation.


## 3. Comparison with Self-Similar Solution and Numerical Results

In this section, we do comparison of our approximation, described in Section2, with both the self-similar method and a direct numerical integration of (1.4), (1.9), and (1.10). First we consider the case where (3.1)ke(E)=ke0En,(3.2)E0=Λ0exp(2αt).
This is a well-known problem studied by Marshak using self-similar method in [10]. Here, we derive and compare analytically our linear approximate solution described in Section 2.1 with the Marshak's self-similar solution. Notice from (2.4) that (3.3)S(E)=ke0En+1n+1,
and from (2.17) and (3.2) that (3.4)G(t)=∫0tke0Λ0n+1n+1exp[2α(n+1)t1]dt1=ke0Λ0n+12α(n+1)2{exp[2α(n+1)t]-1}.
When t is large, we have from (3.2) (3.5)G≈ke0E0n+12α(n+1)2.
Substituting (3.1)–(3.3) and (3.5) into (2.23), we obtain(3.6)ddt(γnE0)=[12-16(1+γn2)(12+γn3)-1]-1×{-2α(n+1)2E0(12+γn3)-α(n+1)E0(1+γn2)+2α[14(1+γn2)(12+γn3)-1-1]E0},
where
(3.7)γn(t)=E1hE0.
A constant solution to (3.6) can be determined by(3.8)γn=[12-16(1+γn2)(12+γn3)-1]-1×{-(n+1)2(12+γn3)-n+12(1+γn2)+[14(1+γn2)(12+γn3)-1-1]}.
Solving for γn, we have (3.9)γn=3n(2n+3)+9n2(2n+3)2-24(n+1)2(n+2)8(n+1)2-32.
Substituting (3.9) and (3.5) into (2.16) yields (3.10)h2=123n(2n+3)+9n2(2n+3)2-24(n+1)2(n+2)ke0E0nα
or (3.11)h=[12n23n(2n+3)+9n2(2n+3)2-24(n+1)2(n+2)]1/2hMar≡CnhMar,
where the Marshak's solution, hMar, is given by(3.12)hMar=1nke0E0nα.
Since we have (3.13)[12n23n(2n+3)+9n2(2n+3)2-24(n+1)2(n+2)]1/2→1,asn→∞,
equation (3.11) indicates that our linear approximate solution tends to Marshak's result when n becomes large and in general we have (3.14)1>hhMar>0.81,forn≥1.2.
In Figure 2, we plot Cn against n to show the limiting behavior as given in (3.13).Figure 2
Ratio between the linear approximation and the Marshak's solution, as a function ofn.Now we perform numerical comparisons. These include comparisons of our approximate solutions with self-similar solution and with the numerically integrated solutions. Four examples will be given. In the first example, we choose, in (3.1) and (3.2), ke0=Λ0=1, n=4, and α=0.9 so that ke(E)=E4 and E0(t)=exp(1.8t). For this example, both our linear approximation and Marshak's self-similar solution are suitable and are given in (3.11) and (3.12), respectively. In Figure 3, we show the heat wave front h as a function of time from numerically integrated (solid line), linear approximate (dashed line), quadratic approximate (dotted line), and Marshak's (dotted-dashed line) solutions. Both Marshak's and our two approximate solutions agree well with the numerically integrated result. Marshak's solution is perhaps slightly more accurate than our linear approximation. However, we see a clear improvement of our quadratic approximation over the linear case and it becomes comparable in accuracy with the Marshak's solution. Here we point out that our approximation method has advantage in two aspects. First, our method has already been systematically extended to higher-order approximation (see Section 2.2) and we have shown clear improvement in accuracy from linear to quadratic approximations. Second, notice that the self-similar method that Marshak used to derive his solution is only suitable for the conductivity being a power function of energy density and boundary energy density being either an exponential [10] or a power function [3] of time. Our approximation method described in Section 2 does not have these restrictions. In particular, as mentioned earlier, our method is suitable for many different function forms describing the dependency of heat conductivity on energy density, as long as it equals zero when energy density is zero. To demonstrate this idea, we choose in the second example, a constant boundary energy density E0(t)=1 and an effective conductivity(3.15)ke(E)=10E2.5(1+10E2.57)exp(2E2.5)+40E4.5(1+9E4.511)exp(E4.5).
This effective conductivity is remarkably different from pure power law of energy density. Therefore, self-similar method is not suitable (and Marshak's solution is not available). However, our approximate solutions are still available. The third example corresponds to effective conductivity, ke(E)=400E4, and boundary energy density, (3.16)E0(t)=t-sin(6πt)6π.
This time, effective conductivity is a pure power law of energy density, but the boundary energy density is neither a pure exponential nor a power law function of time. Finally, in order to construct a stringent test case, we combine, in the fourth example, the effective conductivity from the second example and boundary energy density from the third example as given in (3.15) and (3.16), respectively. Again, self-similar method is not suitable for both of these last two examples; while our approximation method is. In Figures 4, 5, and 6, we show the heat wave front as a function of time from the numerically integrated (solid line), the linear approximate (dashed line), and the quadratic approximate (dotted line) solutions for examples 2, 3, and 4, respectively. One observes a good agreement of our two approximate solutions to the numerically integrated result for all three examples, including the stringent test case of the example 4. This demonstrates that our approach can provide approximate solutions of nonlinear heat conduction equation for more general conductivity and boundary energy density.Figure 3
Heat wave front for numerical integrated (solid), linear approximate (dashed), quadratic approximate (dotted), and Marshak's (dotted-dashed) solutions as a function of time for the case whereE0(t)=exp(1.8t) and ke(E)=E4.Figure 4
Heat wave front for numerical integrated (solid), linear approximate (dashed), and quadratic approximate (dotted) solutions as a function of time for the case whereE0(t)=1 and effective conductivity is given by (3.15).Figure 5
Same as Figure4 for the case where ke(E)=400E4 and boundary energy density is given by (3.16).Figure 6
Same as Figure4 for the case boundary energy density is given by (3.16) and effective conductivity given by (3.15).Next, in Figure7 we show velocity values of the heat wave as a function of time from numerically integrated (solid line), the linear approximate (dashed line), and the quadratic approximate (dotted line) solutions for the third example. The boundary energy density for this example, given by (3.16), is demonstrated in Figure 8. It is interesting to note that the behavior of velocity of heat wave is strongly influenced by the boundary condition. As shown in Figures 7 and 8, when the boundary energy density increases rapidly with time, the heat wave will accelerate; when the boundary energy density increases slowly with time, the velocity of heat wave will decrease. Finally, in order to view the evolution of the nonlinear heat wave over the entire field for this example, we plot, in Figure 9, the numerical solution of the energy density versus x at different times. Figure 10 is a similar plot for the solution of energy density in example 4. In both cases, the nonlinear heat wave accelerates and decelerates according to the rapid and slow increases of the boundary energy density at the origin, respectively.Figure 7
Wave-front velocity as a function of time for the case in Figure5.Figure 8
Boundary energy density as a function of time for the case in Figure7.Figure 9
Evolution of heat wave over the entire field for the case in Figure5.Figure 10
Evolution of heat wave over the entire field for the case in Figure6.


## 4. Conclusion

We have developed an approximation method for solving nonlinear heat conduction problem. Comparison with both self-similar and numerical integrated solutions confirms the accuracy of our approximate results and it also indicates that our approximation method is suitable and efficient for more general and realistic cases since it assumes no restriction on the form of boundary energy density and heat conductivity. The consistent improvement of the quadratic approximation over the linear one observed in all of the Figures3–7 indicates that our method can be systematically extended to achieve higher order of accuracy. Finally, our approximation results also reveal a strong dependence of the heat wave velocity on the boundary energy density.Through our method, we have reduced a mathematical problem with a partial differential equation, which may be considered as an infinite number of ordinary differential equations, to a set of finite number of ordinary differential equations. In particular, for the linear approximation given in Section2.1, the problem is reduced to only one nonlinear ordinary differential equation (2.23) and an algebraic equation (2.16). These provide opportunity for further study of the nonlinear dynamics of the problem. The mathematical formulas, such as (2.16) resulted from our approximation method, could also useful in engineering practices.The idea in our approximation method is quite general, as it can be applied to any nonlinear physical processes in which the solution exhibits wave front behavior. This makes our approximation methodology possibly be useful in many application fields, including fluid dynamics, combustion, and environmental and material sciences.



---

*Source: 421657-2010-08-10.xml*