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

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.

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) where 𝜌 is density, 𝑐_𝑝 is the specific heat at constant pressure, and 𝑘 is thermal conductivity. In this paper, we assume all of them depend only on temperature 𝑇(𝑥,𝑡). Following Zel'dovich and Raizer [3], we introduce internal energy per unit volume 𝐸(𝑇)=^𝑇₀𝜌(𝜃)𝑐_𝑝(𝜃)𝑑𝜃, (1.2) then the temperature can be solved as a function of the energy density 𝐸(𝑥,𝑡)𝑇(𝑥,𝑡)=𝑇(𝐸), (1.3) and the nonlinear heat conduction equation can be rewritten in terms of 𝐸,𝜕𝐸(𝑥,𝑡)=𝜕𝜕𝑡𝜕𝑥(𝑘_𝑒(𝐸)𝜕𝐸𝜕𝑥), (1.4) where 𝑘_𝑒(𝐸)=𝑘(𝐸)𝑑𝑇𝑑𝐸 (1.5) 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 𝑇 and the energy density 𝐸.

Figure 1 

Schematic of a semi-infinite body.

With the heat flux 𝐽(𝑥,𝑡) defined as 𝐽(𝑥,𝑡)=−𝑘_𝑒(𝐸)𝜕𝐸,𝜕𝑥 (1.6) equation (1.4) can be equivalently written as 𝜕𝐸𝜕𝑡=−𝜕𝐽(𝑥,𝑡).𝜕𝑥 (1.7) We choose the initial and boundary conditions as 𝑇(𝑥,0)=0for𝑥>0,𝑇(0,𝑡)=𝑇₀(𝑡) (1.8) or correspondingly𝐸(𝑥,0)=0,for𝑥>0,𝐸(0,𝑡)=𝐸₀(𝑡). (1.9)(1.10) Also, we assume that the effective heat conductivity 𝑘_𝑒(𝐸) tends to zero with the energy density 𝐸; therefore, the thermal wave has a finite extension and a definite front [3]. With 𝑥=ℎ(𝑡) as the front coordinate of thermal wave at time 𝑡, we have 𝐸(𝑥,𝑡)|_{𝑥≥ℎ(𝑡)}=𝑘_𝑒(𝐸)|_{𝑥≥ℎ(𝑡)}=𝐽(𝑥,𝑡)|_{𝑥≥ℎ(𝑡)}=0. (1.11)

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 density 𝐸 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., 𝑘_𝑒=𝑘₀(𝐸+𝑏𝐸²+⋯)). 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 𝑥, one finds 𝜕[𝜕𝑡^𝑥₀𝐸(𝜁,𝑡)𝑑𝜁]=𝐽₀(𝑡)−𝐽(𝑥,𝑡), (2.1) where 𝐽₀(𝑡)=𝐽(0,𝑡)=−𝑘_𝑒(𝐸₀)𝐸₁𝐸(𝑡),₁(𝑡)=𝜕𝐸(𝑥,𝑡)|𝜕𝑥_𝑥=0. (2.2)(2.3) We define 𝑆(𝐸)=^𝐸₀𝑘_𝑒(𝜁)𝑑𝜁, (2.4) and from (1.6) we may write the heat flux as 𝐽(𝑥,𝑡)=−𝜕𝑆(𝐸).𝜕𝑥 (2.5) Also it can be verified that 𝑆(𝑥,𝑡)|_{𝑥=ℎ(𝑡)}=𝑆(𝐸)|_𝐸=0=0. (2.6) Substituting (2.5) into (2.1) and integrating (2.1) from 0 to 𝑥, we obtain 𝜕[𝜕𝑡^𝑥₀𝑑𝜂^𝜂₀𝐸(𝜁,𝑡)𝑑𝜁]=𝑆(𝑥,𝑡)−𝑆(𝐸₀)+𝐽₀(𝑡)𝑥. (2.7) Noting that𝜕[𝜕𝑡^𝑥₀𝑑𝜂^𝜂₀𝜕𝐸(𝜁,𝑡)𝑑𝜁]=[𝜕𝑡^𝑥₀𝜕(𝑥−𝜁)𝐸(𝜁,𝑡)𝑑𝜁]=𝑥[𝜕𝑡^𝑥₀𝜕𝐸(𝜁,𝑡)𝑑𝜁]−[𝜕𝑡^𝑥₀𝜁𝐸(𝜁,𝑡)𝑑𝜁], (2.8) then substituting (2.8) and (2.1) into (2.7) results in 𝜕[𝜕𝑡^𝑥₀𝜁𝐸(𝜁,𝑡)𝑑𝜁]=𝑆(𝐸₀)−𝑆(𝑥,𝑡)−𝐽(𝑥,𝑡)𝑥. (2.9) Notice that (2.9) is an equivalent integral form of the nonlinear heat conduction equation (1.4).

Substituting 𝑥=ℎ(𝑡) and the boundary condition (1.11) and (2.6) into (2.9) and (2.1) gives 𝑑[𝑑𝑡^ℎ₀𝜁𝐸(𝜁,𝑡)𝑑𝜁]=𝑆(𝐸₀𝑑),[𝑑𝑡^ℎ₀𝐸(𝜁,𝑡)𝑑𝜁]=𝐽₀(𝑡). (2.10)(2.11) 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 𝐸(𝑥,𝑡) about the boundary point 𝑥=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 𝐸(𝑥,𝑡). 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 𝑘(𝑇) is a strong nonlinear function of temperature 𝑇, or correspondingly, when the effective heat conductivity 𝑘_𝑒(𝐸) is a strong nonlinear function of energy density 𝐸, 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 (𝑥=0) to the wave front (𝑥=ℎ(𝑡)). 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 𝐸(𝑥,𝑡) as a linear function of 𝑥: 𝐸(𝑥,𝑡)≈𝐸₀(𝑡)+𝐸₁(𝑡)𝑥,0≤𝑥≤ℎ(𝑡), (2.12) where ℎ(𝑡) is the wave front, and 𝐸₀(𝑡) and 𝐸₁(𝑡) are defined in (1.10) and (2.3). Notice that 𝐸₀(𝑡) is known from the boundary condition (1.10) so that we have two unknown functions of time, 𝐸₁(𝑡) and ℎ(𝑡) 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𝐸₁ℎ³)=𝑆(𝐸₀𝑑),𝑑𝑡(𝐸₀𝐸ℎ+¹2ℎ²)=−𝑘_𝑒(𝐸₀)𝐸₁. (2.13)(2.14) Integrating (2.13) and using the initial condition ℎ(𝑡)|_𝑡=0=𝐸₁ℎ³|_𝑡=0=0, (2.15) we may solve for ℎ² in the form ℎ²𝐸=𝐺(𝑡)(⁰2+𝐸₁ℎ3)⁻¹, (2.16) where 𝐺(𝑡)=^𝑡₀𝑆(𝐸₀(𝑡₁))𝑑𝑡₁. (2.17) From (2.14), we haveℎ𝑑𝑑𝑡[ℎ(𝐸₀+𝐸₁ℎ2)]=−𝑘_𝑒(𝐸₀)𝐸₁ℎ (2.18) orℎ²𝑑𝑑𝑡(𝐸₀+𝐸₁ℎ21)+2(𝐸₀+𝐸₁ℎ2)𝑑(ℎ²)𝑑𝑡=−𝑘_𝑒(𝐸₀)𝐸₁ℎ. (2.19) Using (2.16), we get 𝐸𝐺(𝑡)(⁰2+𝐸₁ℎ3)⁻¹𝑑𝑑𝑡(𝐸₀+𝐸₁ℎ21)+2(𝐸₀+𝐸₁ℎ2)𝑑𝐸𝑑𝑡[𝐺(𝑡)(⁰2+𝐸₁ℎ3)⁻¹]=−𝑘_𝑒(𝐸₀)𝐸₁ℎ (2.20) or𝐸𝐺(𝑡)(⁰2+𝐸₁ℎ3)⁻¹[12−16(𝐸₀+𝐸₁ℎ2𝐸)(⁰2+𝐸₁ℎ3)⁻¹]𝑑𝑑𝑡(𝐸₁ℎ)=−𝑘_𝑒(𝐸₀)𝐸₁1ℎ−2(𝐸₀+𝐸₁ℎ2𝐸)(⁰2+𝐸₁ℎ3)⁻¹[𝑑𝐺(𝑡)−𝑑𝑡𝐺(𝑡)2(𝐸₀2+𝐸₁ℎ3)⁻¹𝑑𝐸₀]𝐸𝑑𝑡−𝐺(𝑡)(⁰2+𝐸₁ℎ3)⁻¹𝑑𝐸₀.𝑑𝑡 (2.21) Substituting (2.17) into (2.21) yields[12−16(𝐸₀+𝐸₁ℎ2𝐸)(⁰2+𝐸₁ℎ3)⁻¹]𝑑𝑑𝑡(𝐸₁𝑘ℎ)=−^𝑒(𝐸₀)𝐸₁ℎ(𝐸𝐺(𝑡)2⁰+𝐸₁ℎ31)−2𝐺(𝑡)(𝐸₀+𝐸₁ℎ2)[𝑆(𝐸₀)−𝐺(𝑡)2(𝐸₀2+𝐸₁ℎ3)⁻¹𝑑𝐸₀𝑑𝑡]−𝑑𝐸₀𝑑𝑡 (2.22) or𝑑𝑑𝑡(𝐸₁1ℎ)=[2−16(𝐸₀+𝐸₁ℎ2𝐸)(⁰2+𝐸₁ℎ3)⁻¹]⁻¹𝑘×{−^𝑒(𝐸₀)𝐸₁ℎ(𝐸𝐺(𝑡)2⁰+𝐸₁ℎ31)−2𝐺(𝑡)(𝐸₀+𝐸₁ℎ2)𝑆(𝐸₀)1+[4(𝐸₀+𝐸₁ℎ2𝐸)(⁰2+𝐸₁ℎ3)⁻¹−1]𝑑𝐸₀𝑑𝑡}. (2.23)

Equation (2.23) is a nonlinear ordinary differential equation for 𝐸₁ℎ. It is easy to understand that the approximate expansion (2.12) is suitable only for the case of ℎ(𝑡)>0. From (2.15), (2.23) is valid only for 𝑡>0. For the initial time, we may simply assume that 𝑑𝑑𝑡(𝐸₁ℎ)|_𝑡=0=𝐸₁ℎ|_𝑡=0=0. (2.24)

With 𝐸₀ given in the boundary condition (1.10), 𝐸₁ℎ can be solved from (2.23), (2.24), and (2.17). Then the heat wave front ℎ will be available from (2.16) and 𝐸₁ from 𝐸₁ℎ and ℎ. Finally, the approximate values of the energy density 𝐸 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 Section 2.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 𝑥 and we can obtain more accurate wave fronts and energy density profiles. In general, we may expand the energy density 𝐸(𝑥,𝑡) as a Taylor polynomial of degree up to 𝑀 (𝑀≥2) in 𝑥: 𝐸(𝑥,𝑡)=^𝑀_𝑗=0𝐸_𝑗(𝑡)𝑥^𝑗𝐸,0≤𝑥≤ℎ(𝑡),_𝑗1(𝑡)=𝜕𝑗!^𝑗𝐸(𝑥,𝑡)𝜕𝑥^𝑗|_𝑥=0. (2.25)(2.26) In this case, we have 𝑀+1 unknowns, 𝐸₁(𝑡),𝐸₂(𝑡),…,𝐸_𝑀(𝑡), and ℎ(𝑡). (𝐸₀(𝑡) is known from the boundary condition (1.10).) Substituting (2.25) and (2.2) into (2.10) and (2.11) provides two equations:𝑑[𝑑𝑡^𝑀_𝑗=0𝐸_𝑗(𝑡)ℎ𝑗+2^𝑗+2]=𝑆(𝐸₀𝑑),[𝑑𝑡^𝑀_𝑗=0𝐸_𝑗(𝑡)ℎ𝑗+1^𝑗+1]=−𝑘_𝑒(𝐸₀)𝐸₁(𝑡). (2.27)(2.28) In order to close the system, we need (𝑀−1) more equations. Using (2.4), the heat conduction equation (1.4) can be written as 𝜕𝐸(𝑥,𝑡)=𝜕𝜕𝑡²𝑆(𝑥,𝑡)𝜕𝑥². (2.29) For positive integer 𝑙=2,…,𝑀, we study 𝑑[𝑑𝑡^ℎ₀𝜁^𝑙𝐸(𝜁,𝑡)𝑑𝜁]. (2.30) We show in Appendix A that, together with (2.27) and (2.28), this approach leads to the following 𝑀+1 equations: ℎ²=𝐺(𝑡)(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻¹_𝑀𝑗=1,[1−1𝑗+1(2(𝑗+2)^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸_𝑗ℎ^𝑗)𝑑𝑡=−𝑆₁1(𝑡)ℎ+[4(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1]ℎ²𝑑𝐸₀−1𝑑𝑡2(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀),_𝑀𝑗=1[1−𝑗+𝑙+1𝑙+1(2(𝑗+2)^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+𝑙+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸_𝑗ℎ^𝑗)𝑑𝑡=[𝑙+14(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+𝑙+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1𝑙+1]ℎ²𝑑𝐸₀−𝑑𝑡𝑙+12(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+𝑙+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀)+𝑙(𝑙−1)^𝑀+1_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗𝑗+𝑙−1,𝑙=2,…,𝑀, (2.31)(2.32)(2.33) where 𝑆_𝑗1(𝑡)=𝜕𝑗!^𝑗𝑆(𝑥,𝑡)𝜕𝑥^𝑗|_𝑥=0=𝑆_𝑗(𝐸₀,𝐸₁⋯𝐸_𝑗),0≤𝑗≤𝑀, (2.34) and 𝑆_𝑀+1(𝑡) is determined by ^𝑀+1_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗=0. (2.35)

Equations (2.31), (2.32), and (2.33) form a close set for the 𝑀+1 unknowns 𝐸₁(𝑡)⋯𝐸_𝑀(𝑡), and ℎ(𝑡). Once 𝐸₁(𝑡)⋯𝐸_𝑀(𝑡), and ℎ(𝑡) are available, the energy density profile can be determined from (2.25) and (1.11).

For the case of a quadratic approximation, 𝑀=2, as shown in Appendix B, we may getℎ²𝐸=𝐺(𝑡)(⁰2+𝐸₁ℎ3+𝐸₂ℎ²)𝑗+2⁻¹,𝑑𝑑𝑡(𝐸₁1ℎ)=ℎ²𝐶₀𝐵₂−𝐶₂𝐵₀𝐴₀𝐵₂−𝐴₂𝐵₀=1(𝐺(𝑡)²_𝑗=0𝐸_𝑗ℎ^𝑗)𝐶𝑗+2⁰𝐵₂−𝐶₂𝐵₀𝐴₀𝐵₂−𝐴₂𝐵₀𝑑,𝑑𝑡(𝐸₂ℎ²1)=ℎ²𝐴₀𝐶₂−𝐴₂𝐶₀𝐴₀𝐵₂−𝐴₂𝐵₀=1(𝐺(𝑡)²_𝑗=0𝐸_𝑗ℎ^𝑗)𝐴𝑗+2⁰𝐶₂−𝐴₂𝐶₀𝐴₀𝐵₂−𝐴₂𝐵₀, (2.36)(2.37)(2.38) where𝐴₀=12−16(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹,𝐵₀=13−18(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹,𝐶₀=−𝑆₁1ℎ+[4(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1]ℎ²𝑑𝐸₀−1𝑑𝑡2(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀𝐴),₂=14−12(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹,𝐵₂=15−38(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹,𝐶₂3=[4(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−13]ℎ²𝑑𝐸₀−3𝑑𝑡2(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀)+2³_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗,𝑗+1 (2.39)(2.40)(2.41)(2.42) and the 𝑆₁ℎ, 𝑆₂ℎ², and 𝑆₃ℎ³ terms that appeared in (2.40) and (2.42) are given by 𝑆₁ℎ=ℎ𝜕𝑆(𝑥,𝑡)|𝜕𝑥_𝑥=0=𝑘_𝑒(𝐸₀)𝐸₁𝑆ℎ,₂ℎ²=ℎ²2𝜕²𝑆(𝑥,𝑡)𝜕𝑥²|_𝑥=0=12𝑘^′_𝑒(𝐸₀)(𝐸₁ℎ)²+𝑘_𝑒(𝐸₀)𝐸₂ℎ²,𝑆₃ℎ³=−[𝑆₂ℎ²+𝑆₁ℎ+𝑆(𝐸₀)]. (2.43) Here (2.36), (2.37), and (2.38), forming a close set of equations for ℎ², 𝐸₁, and 𝐸₂, 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 Section 2, 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 𝑘_𝑒(𝐸)=𝑘_𝑒0𝐸^𝑛,𝐸₀=Λ₀exp(2𝛼𝑡). (3.1)(3.2) 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 𝑘𝑆(𝐸)=^𝑒0𝐸^𝑛+1,𝑛+1 (3.3) and from (2.17) and (3.2) that 𝐺(𝑡)=^𝑡₀𝑘_𝑒0Λ₀^𝑛+1𝑛+1exp[2𝛼(𝑛+1)𝑡₁]𝑑𝑡₁=𝑘_𝑒0Λ₀^𝑛+12𝛼(𝑛+1)²{exp[2𝛼(𝑛+1)𝑡]−1}. (3.4) When 𝑡 is large, we have from (3.2) 𝑘𝐺≈^𝑒0𝐸₀^𝑛+12𝛼(𝑛+1)². (3.5) Substituting (3.1)–(3.3) and (3.5) into (2.23), we obtain𝑑𝑑𝑡(𝛾_𝑛𝐸₀1)=[2−16𝛾(1+^𝑛21)(2+𝛾_𝑛3)⁻¹]⁻¹×{−2𝛼(𝑛+1)²𝐸₀(12+𝛾_𝑛3)−𝛼(𝑛+1)𝐸₀𝛾(1+^𝑛2)1+2𝛼[4𝛾(1+^𝑛21)(2+𝛾_𝑛3)⁻¹−1]𝐸₀}, (3.6) where 𝛾_𝑛𝐸(𝑡)=¹ℎ𝐸₀. (3.7) A constant solution to (3.6) can be determined by𝛾_𝑛1=[2−16𝛾(1+^𝑛21)(2+𝛾_𝑛3)⁻¹]⁻¹×{−(𝑛+1)²(12+𝛾_𝑛3)−𝑛+12𝛾(1+^𝑛21)+[4𝛾(1+^𝑛21)(2+𝛾_𝑛3)⁻¹−1]}. (3.8) Solving for 𝛾_𝑛, we have 𝛾_𝑛=3𝑛(2𝑛+3)+9𝑛²(2𝑛+3)²−24(𝑛+1)²(𝑛+2)8(𝑛+1)²−32. (3.9) Substituting (3.9) and (3.5) into (2.16) yields ℎ²=123𝑛(2𝑛+3)+9𝑛²(2𝑛+3)²−24(𝑛+1)²𝑘(𝑛+2)^𝑒0𝐸^𝑛₀𝛼 (3.10) or ℎ=[12𝑛²3𝑛(2𝑛+3)+9𝑛²(2𝑛+3)²−24(𝑛+1)²](𝑛+2)^1/2ℎ^Mar≡𝐶_𝑛ℎ^Mar, (3.11) where the Marshak's solution, ℎ^Mar, is given byℎ^Mar=1𝑛𝑘_𝑒0𝐸^𝑛₀𝛼. (3.12) Since we have [12𝑛²3𝑛(2𝑛+3)+9𝑛²(2𝑛+3)²−24(𝑛+1)²](𝑛+2)^1/2→1,as𝑛→∞, (3.13) equation (3.11) indicates that our linear approximate solution tends to Marshak's result when 𝑛 becomes large and in general we have ℎ1>ℎ^Mar>0.81,for𝑛≥1.2. (3.14) In Figure 2, we plot 𝐶_𝑛 against 𝑛 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 of $n$.

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), 𝑘_𝑒0=Λ₀=1, 𝑛=4, and 𝛼=0.9 so that 𝑘_𝑒(𝐸)=𝐸⁴ and 𝐸₀(𝑡)=exp(1.8𝑡). 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 ℎ 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 𝐸₀(𝑡)=1 and an effective conductivity𝑘_𝑒(𝐸)=10𝐸^2.5(1+10𝐸^2.57)exp(2𝐸^2.5)+40𝐸^4.5(1+9𝐸^4.511)exp(𝐸^4.5). (3.15) 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, 𝑘_𝑒(𝐸)=400𝐸⁴, and boundary energy density, 𝐸₀(𝑡)=√√𝑡−sin(6𝜋𝑡).6𝜋 (3.16) 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 where $E_0(t)=\exp \hspace{0.17em}(1.8t)$ and $k_e(E)=E^4$.

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 where $E_0(t)=1$ and effective conductivity is given by (3.15).

Figure 5 

Same as Figure 4 for the case where $k_e(E)=400E^4$ and boundary energy density is given by (3.16).

Figure 6 

Same as Figure 4 for the case boundary energy density is given by (3.16) and effective conductivity given by (3.15).

Next, in Figure 7 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 𝑥 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 Figure 5.

Figure 8 

Boundary energy density as a function of time for the case in Figure 7.

Figure 9 

Evolution of heat wave over the entire field for the case in Figure 5.

Figure 10 

Evolution of heat wave over the entire field for the case in Figure 6.

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 Figures 3–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 Section 2.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.

Appendices

A. Derivation of Governing Equations for Higher-Order Approximation

We show here that we can obtain 𝑀−1 equations from (2.30). Together with (2.27) and (2.28), we obtain 𝑀+1 equations.

For positive integers 𝑙=2,…,𝑀, we do integration by part in (2.30) and use (2.29) to obtain𝑑[𝑑𝑡^ℎ₀𝑥^𝑙𝐸(𝑥,𝑡)𝑑𝑥]=ℎ^𝑙𝐸(ℎ,𝑡)𝑑ℎ+𝑑𝑡^ℎ₀𝑥^𝑙𝜕𝐸(𝑥,𝑡)𝜕𝑡𝑑𝑥=ℎ^𝑙𝐸(ℎ,𝑡)𝑑ℎ+𝑑𝑡^ℎ₀𝑥^𝑙𝜕²𝑆(𝑥,𝑡)𝜕𝑥²𝑑𝑥=ℎ^𝑙𝐸(ℎ,𝑡)𝑑ℎ𝑑𝑡+𝑥^𝑙𝜕𝑆(𝑥,𝑡)|𝜕𝑥^𝑥=ℎ_𝑥=0−𝑙^ℎ₀𝑥^𝑙−1𝜕𝑆(𝑥,𝑡)𝜕𝑥𝑑𝑥=ℎ^𝑙𝐸(ℎ,𝑡)𝑑ℎ𝑑𝑡+𝑥^𝑙𝐽(𝑥,𝑡)|^𝑥=ℎ_𝑥=0−𝑙𝑥^𝑙−1𝑆(𝑥,𝑡)|^𝑥=ℎ_𝑥=0+𝑙(𝑙−1)^ℎ₀𝑥^𝑙−2𝑆(𝑥,𝑡)𝑑𝑥=ℎ^𝑙𝐸(ℎ,𝑡)𝑑ℎ𝑑𝑡+ℎ^𝑙𝐽(𝑥,𝑡)|_𝑥=ℎ−𝑙ℎ^𝑙−1𝑆(𝑥,𝑡)|_𝑥=ℎ+𝑙(𝑙−1)^ℎ₀𝑥^𝑙−2𝑆(𝑥,𝑡)𝑑𝑥. (A.1) Substituting (1.11), (2.6), and (2.25) into (A.1) yields 𝑑[𝑑𝑡^𝑀_𝑗=0𝐸_𝑗(𝑡)ℎ^𝑗+𝑙+1𝑗+𝑙+1]=𝑙(𝑙−1)^ℎ₀𝑥^𝑙−2𝑆(𝑥,𝑡)𝑑𝑥,𝑙=2,…,𝑀. (A.2) Noting that from (1.9) 𝑆(𝑥,𝑡) is a function of energy density 𝐸, we may expand 𝑆(𝑥,𝑡) as a polynomial of 𝑥. Since𝑘𝑆(𝑥,𝑡)=𝑆(𝐸),_𝑒(𝐸)=𝑑𝑆(𝐸),𝑑𝐸𝜕𝑆(𝑥,𝑡)𝜕𝑥=𝑘_𝑒(𝐸)𝜕𝐸(𝑥,𝑡),𝜕𝜕𝑥₂𝑆(𝑥,𝑡)𝜕𝑥²=𝜕𝜕𝑥[𝑘_𝑒(𝐸)𝜕𝐸(𝑥,𝑡)𝜕𝑥]=𝜕𝑘_𝑒(𝐸)𝜕𝑥𝜕𝐸(𝑥,𝑡)𝜕𝑥+𝑘_𝑒𝜕(𝐸)²𝐸(𝑥,𝑡)𝜕𝑥²=𝑘^′_𝑒(𝐸)[𝜕𝐸(𝑥,𝑡)]𝜕𝑥²+𝑘_𝑒𝜕(𝐸)²𝐸(𝑥,𝑡)𝜕𝑥², (A.3) where 𝑘^′_𝑒(𝐸)≡𝑑𝑘_𝑒(𝐸),𝑑𝐸 (A.4) so 𝑆(𝑥,𝑡)|_𝑥=0=𝑆(𝐸(0,𝑡))=𝑆(𝐸₀),𝜕𝑆(𝑥,𝑡)|𝜕𝑥_𝑥=0=𝑘_𝑒(𝐸₀)𝐸₁1(𝑡),2𝜕²𝑆(𝑥,𝑡)𝜕𝑥²|_𝑥=0=12𝑘^′_𝑒(𝐸₀)𝐸²₁(𝑡)+𝑘_𝑒(𝐸₀)𝐸₂(𝑡). (A.5) Here we have used (2.26). It is easy to show that any of 𝜕^𝑗𝑆(𝑥,𝑡)/𝜕𝑥^𝑗|_𝑥=0 can be represented as a function of 𝐸₀(𝑡), 𝐸₁(𝑡)⋯𝐸_𝑗(𝑡). With 𝐸₀(𝑡), 𝐸₁(𝑡)⋯𝐸_𝑀(𝑡), and boundary condition (2.6) given, we may approach 𝑆(𝑥,𝑡) as 𝑆(𝑥,𝑡)=^𝑀+1_𝑗=0𝑆_𝑗(𝑡)𝑥^𝑗,0≤𝑥≤ℎ(𝑡), (A.6) where 𝑆_𝑗1(𝑡)=𝜕𝑗!^𝑗𝑆(𝑥,𝑡)𝜕𝑥^𝑗|_𝑥=0=𝑆_𝑗(𝐸₀,𝐸₁⋯𝐸_𝑗),for0≤𝑗≤𝑀, (A.7) and 𝑆_𝑀+1(𝑡) is determined by ^𝑀+1_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗=0. (A.8)

Substituting (A.6) into (A.2) gives𝑑[𝑑𝑡^𝑀_𝑗=0𝐸_𝑗(𝑡)ℎ^𝑗+𝑙+1𝑗+𝑙+1]=𝑙(𝑙−1)^𝑀+1_𝑗=0𝑆_𝑗(𝑡)ℎ^{𝑗+𝑙−1}𝑗+𝑙−1,𝑙=2,…,𝑀, (A.9) or ℎ²𝑑(𝑑𝑡^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗𝑗+𝑙+1)+𝑙+12(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+𝑙+1𝑑(ℎ²)𝑑𝑡=𝑙(𝑙−1)^𝑀+2_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗𝑗+𝑙−1,𝑙=2,…,𝑀. (A.10) Integrating (2.27) and using (2.17) yields ℎ²=𝐺(𝑡)(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻¹, (A.11) so𝑑(ℎ²)𝑑𝑡=𝑆(𝐸₀)(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻¹−𝐺(𝑡)(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻²(^𝑀_𝑗=01𝑗+2𝑑(𝐸_𝑗ℎ^𝑗))𝑑𝑡=𝑆(𝐸₀)(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻¹−ℎ²2(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻¹𝑑𝑇₀𝑑𝑡−ℎ²(^𝑀_𝑗=0𝐸_𝑗ℎ^𝑗)𝑗+2⁻¹(^𝑀_𝑗=11𝑗+2𝑑(𝐸_𝑗ℎ^𝑗)𝑑𝑡). (A.12) Substituting (A.11) and (A.12) into (A.10), we get^𝑀_𝑗=1[1−𝑗+𝑙+1𝑙+1(2(𝑗+2)^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+𝑙+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸_𝑗ℎ^𝑗)𝑑𝑡=[𝑙+14(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+𝑙+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1𝑙+1]ℎ²𝑑𝐸₀−𝑑𝑡𝑙+12(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+𝑙+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀)+𝑙(𝑙−1)^𝑀+1_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗𝑗+𝑙−1,𝑙=2,…,𝑀. (A.13)

Similarly, substituting (A.11) and (A.12) into (2.28) yields ^𝑀_𝑗=1[1−1𝑗+1(2(𝑗+2)^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸_𝑗ℎ^𝑗)𝑑𝑡=−𝑆₁1(𝑡)ℎ+[4(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1]ℎ²𝑑𝐸₀−1𝑑𝑡2(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(^𝑀_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀). (A.14) Equations (A.11), (A.14), and (A.13) are equations (2.31), (2.32), and (2.33), respectively.

B. Derivation of Governing Equations for Quadratic Approximation

For the case of quadratic approximation, 𝑀=2. Therefore, 𝑙=2 only in (2.33). In particular, (2.33) and (2.32) reduce, respectively, to[14−12(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸₁ℎ)1𝑑𝑡+[5−38(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸₂ℎ²)𝑑𝑡=2⁴_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗3𝑗+1+[4(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−13]ℎ²𝑑𝐸₀−3𝑑𝑡2(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀[1),2−16(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸₁ℎ)1𝑑𝑡+[3−18(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹]ℎ²𝑑(𝐸₂ℎ²)𝑑𝑡=−𝑆₁1(𝑡)ℎ+[4(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1]ℎ²𝑑𝐸₀−1𝑑𝑡2(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀). (B.1)(B.2) Setting 𝑀=2 in (2.31), we find that ℎ²𝐸=𝐺(𝑡)(⁰2+𝐸₁ℎ3+𝐸₂ℎ²)𝑗+2⁻¹. (B.3) Define𝐴₀=12−16(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝐵₀,=13−18(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝐶₀,=−𝑆₁1ℎ+[4(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−1]ℎ²𝑑𝐸₀−1𝑑𝑡2(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+1)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀𝐴),₂=14−12(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝐵₂=15−38(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝐶₂,3=[4(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹−13]ℎ²𝑑𝐸₀−3𝑑𝑡2(²_𝑖=0𝐸_𝑖ℎ^𝑖𝑖+3)(²_𝑖=0𝐸_𝑖ℎ^𝑖)𝑖+2⁻¹𝑆(𝐸₀)+2³_𝑗=0𝑆_𝑗(𝑡)ℎ^𝑗𝑗+1, (B.4)(1)(2)(3)(4)(5)(6)(7) where 𝑆₁ℎ=ℎ𝜕𝑆(𝑥,𝑡)|𝜕𝑥_𝑥=0=𝑘_𝑒(𝐸₀)𝐸₁𝑆ℎ,₂ℎ²=ℎ²2𝜕²𝑆(𝑥,𝑡)𝜕𝑥²|_𝑥=0=12𝑘^′_𝑒(𝐸₀)(𝐸₁ℎ)²+𝑘_𝑒(𝐸₀)𝐸₂ℎ²,𝑆₃ℎ³=−(𝑆₂ℎ²+𝑆₁ℎ+𝑆(𝐸₀)). (B.5) Solving (𝑑/𝑑𝑡)(𝐸₁ℎ) between (B.1) and (B.2) yields 𝑑𝑑𝑡(𝐸₁1ℎ)=ℎ²𝐶₀𝐵₂−𝐶₂𝐵₀𝐴₀𝐵₂−𝐴₂𝐵₀=1(𝐺(𝑡)²_𝑗=0𝐸_𝑗ℎ^𝑗)𝐶𝑗+2⁰𝐵₂−𝐶₂𝐵₀𝐴₀𝐵₂−𝐴₂𝐵₀. (B.6) Similarly for (𝑑/𝑑𝑡)(𝐸₂ℎ²), 𝑑𝑑𝑡(𝐸₂ℎ²1)=ℎ²𝐴₀𝐶₂−𝐴₂𝐶₀𝐴₀𝐵₂−𝐴₂𝐵₀=1(𝐺(𝑡)²_𝑗=0𝐸_𝑗ℎ^𝑗)𝐴𝑗+2⁰𝐶₂−𝐴₂𝐶₀𝐴₀𝐵₂−𝐴₂𝐵₀. (B.7) Equations (B.3), (B.6), and (B.7) are (2.36), (2.37), and (2.38), respectively.