# Approximate Analytic Solutions of Transient Nonlinear Heat Conduction with Temperature-Dependent Thermal Diffusivity

Authors: M. T. Mustafa; A. F. M. Arif; Khalid Masood
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A new approach for generating approximate analytic solutions of transient nonlinear heat conduction problems is presented. It is based on an effective combination of Lie symmetry method, homotopy perturbation method, finite element method, and simulation based error reduction techniques. Implementation of the proposed approach is demonstrated by applying it to determine approximate analytic solutions of real life problems consisting of transient nonlinear heat conduction in semi-infinite bars made of stainless steel AISI 304 and mild steel. The results from the approximate analytical solutions and the numerical solution are compared indicating good agreement.


---

## Body

## 1. Introduction

Many physical problems involve transient or unsteady heat conduction. Examples include startup or shutdown of power plant, gas turbine, and heat exchanger. For a homogeneous material with constant properties, the Fourier differential equation(1)ρcp∂T∂t=k(∂2T∂x2+∂2T∂y2+∂2T∂z2)
with three physical properties or the equation
(2)∂T∂t=α(∂2T∂x2+∂2T∂y2+∂2T∂z2)
with just one coefficient is solved to predict the temperature distribution during the heating or cooling process. Here k is the thermal conductivity, ρ is the density, cp is the heat capacity, and α=k/ρcp is the thermal diffusivity. Thermal diffusivity α is a material-specific property for characterizing unsteady heat conduction and it describes how quickly a material reacts to a change in temperature. It should be noted that each of these quantities can vary with temperature. Any change in the thermal boundary conditions around a material results in heat flow in or out of the material until thermal equilibrium is achieved. Materials with a high thermal diffusivity will achieve thermal equilibrium faster than those with lower thermal diffusivity.Most metallic materials have thermal properties (thermal conductivity, specific heat, and density) that are usually temperature-dependent. This results in nonlinearities in the governing partial differential equations (PDEs) describing the temperature distribution through these materials. The analysis of nonlinear heat conduction problems is of increasing importance and interest in various engineering and scientific fields. However, no general theory exists for solving nonlinear partial differential equations, and because of the difficulties associated with the solution of nonlinear heat transfer problems, simplifying assumptions are usually made to linearize such problems. For example, constant thermal conductivity is generally assumed for materials having thermal conductivity which varies slightly with temperature. However, if temperature change is substantial or the thermal conductivity varies greatly with temperature, the assumption of constant thermal conductivity may lead to significant error in the solution. Therefore, such situations require solving nonlinear transient heat transfer problems. It is usually difficult and mathematically challenging to obtain exact closed form solutions of such nonlinear transient heat transfer problems. Therefore, a number of numerical methods such as time integration, the finite-difference method (FDM), the finite-element method (FEM), and the boundary-element method (BEM) have been proposed to solve such problems [1–5]. Though the solutions obtained through sophisticated numerical techniques are extremely accurate approximation of the analytic solution, yet the natural tabular form of numerical solutions hampers their real-time applicability unlike solutions in functional form. An alternate is to look for approximate solutions in function form which are more or less as accurate as numerical solutions. Such approximate analytic solutions can be used in real time and can be easily manipulated at different points in space and time, which may result in better insight into physical phenomenon governed by the PDE. With this motivation, the main objective of this work is to implement a new systematic procedure for generating approximate solutions of nonlinear heat conduction PDEs in function form which also meet the accuracy benchmarks of numerical solutions. As test cases, we develop mathematical models for heat conduction in stainless steel AISI 304 and mild steel and apply our procedure to find accurate enough approximate analytic solutions to PDEs governing these real-life problems.While there is no general theory for directly solving nonlinear PDEs, the approach of investigating nonlinear PDEs by transforming to ODEs or simpler PDEs has worked well in sufficient generality, compare [6]. This includes linearization transformations, methods of reduction of PDEs to ODEs, or transformations that result in reduction of the complexity in general. A powerful general technique for analyzing nonlinear PDEs by reducing them to ODEs is given by classical Lie symmetry method, also known as similarity method. PDEs that model physical phenomena naturally inherit symmetries from the underlying physical system. The similarity method takes advantage of these natural symmetries in a PDE and leads to determining special variables, namely, similarity variables, that give rise to the similarity reduction to ODEs. A large amount of literature about the classical Lie symmetry theory and its applications is available, for example, [6–10].The purpose of this paper is to apply a new approach for finding accurate enough approximate analytic solutions of nonlinear PDEs, particularly arising from heat conduction problems. The method utilizes an effective combination of Lie symmetry analysis, homotopy perturbation method, finite element method, and error reduction techniques. A summary of our approach is as follows. Lie symmetries are first utilized to reduce the initial-boundary value problem of PDE to a boundary value problem of ODE. The reduced ODE problem is then simultaneously solved by finite element method and homotopy perturbation method to, respectively, generate numerical solution and initial approximation of the approximate analytic solution. Next, using the numerically generated solution curve as the benchmark for accuracy, an error reduction of the initial approximation curve is carried out to generate explicit approximate solution of the ODE problem. Finally the similarity transformation provides the approximate analytic solution of the IBVP of original PDE.It should be noted that, in general, the solutions of PDEs are surfaces or hypersurfaces while the solutions of ODEs are represented by curves. Furthermore, it is well known that approximating a surface is an increasingly complex and intractable problem. The idea presented here for generating approximate solutions of PDE, that is, approximating the surface solution, practically involves only approximating a curve which is a tractable problem in comparison to the question of approximating surfaces. This makes it a promising approach especially when a reasonably accurate initial approximation of the solution curve of ODE can be obtained.In the next section we develop mathematical models for heat conduction in stainless steel AISI 304 and mild steel. Section3 presents the systematic procedure for generating approximate analytic solutions using a combination of Lie symmetry, finite element, homotopy perturbation, and error reduction methods. Sections 4 and 5 are devoted to application of the method for investigation of heat conduction models for stainless steel AISI 304 and mild steel, respectively.It is assumed that the reader is familiar with standard computations of Lie symmetry and homotopy perturbation method. The reader is referred to [7–9, 11–13], respectively, for an introduction to Lie symmetry method and homotopy perturbation method.


## 2. Formulation of Test Problems

The nonlinear PDE that describes the transient nonlinear heat conduction in a one-dimensional medium is(3)∂T(x,t)∂t=∂∂x(α(T)∂T(x,t)∂x),
where T(x,t) denotes the temperature at a point x at time t and α(T) is the temperature-dependent thermal diffusivity of the material.To apply our method, we will consider test problems for heat conduction in stainless steel AISI 304 and mild steel. For this purpose, the corresponding thermal diffusivity functionsα(T) need to be modeled first. Tables 1 and 2, respectively, give temperature-dependent thermal properties of AISI 394 stainless steel and mild steel considered in this work. Trend lines fitted on these datasets show that linear equation is the best fit for stainless steel as shown in Figure 1, and second-order polynomial is the most accurate fit for mild steel as illustrated in Figure 2. Table 3 gives the estimated thermal diffusivity functions for stainless steel and mild steel, along with the goodness of fit (R2). Using these thermal diffusivity functions, we will apply our method to investigate transient heat conduction in semi-infinite solid. Precisely, in Sections 4 and 5, approximate analytic solutions will be constructed for initial-boundary value problem of the form
(4)∂T(x,t)∂t=∂∂x(α(T)∂T(x,t)∂x),T|t=0=Ti,T|x=0=Ts,T|x=∞=Ti
with the thermal diffusivity functions of stainless steel and mild steel, respectively.Table 1
Thermal properties of AISI 304 stainless steel.
Temperature[K]
Specific heat  (c) [kJ/kg·K]
Conductivity   (k) [W/m K]
Density  (ρ) [kg/m3]
Thermal diffusivity  (α=k/ρc)  (m2/s)
100
0.272
9.2
7963
4.25
E
-
03
300
0.477
14.9
7900
3.95
E
-
03
500
0.536
18.2
7822
4.34
E
-
03
700
0.5695
21.2
7736
4.81
E
-
03
900
0.5965
24
7644
5.26
E
-
03
1100
0.6255
26.7
7550
5.65
E
-
03
1300
0.654
29.23
7455
6.00
E
-
03
1500
0.682
31.7
7362
6.31
E
-
03
1700
0.71
34.17
7269
6.62
E
-
03
1900
0.738
36.63
7175
6.92
E
-
03
2100
0.766
39.1
7082
7.21
E
-
03
2300
0.794
41.57
6989
7.49
E
-
03
2500
0.822
44.03
6895
7.77
E
-
03Table 2
Thermal properties of mild steel.
Temperature[K]
Specific heat  (c) [kJ/kg·K]
Conductivity   (k) [W/m K]
Density  (ρ) [kg/m3]
Thermal diffusivity  (α=k/ρc)  (m2/s)
100
0.346
68.1
7892
2.49
E
-
02
145.8
0.3662
66.36
7883
2.30
E
-
02
191.7
0.3863
64.62
7874
2.12
E
-
02
237.5
0.4065
62.88
7866
1.97
E
-
02
283.3
0.4267
61.13
7857
1.82
E
-
02
329.2
0.4468
59.39
7848
1.69
E
-
02
375
0.467
57.65
7840
1.57
E
-
02
420.8
0.4864
55.79
7827
1.47
E
-
02
466.7
0.505
53.8
7809
1.36
E
-
02
512.5
0.5236
51.81
7791
1.27
E
-
02
558.3
0.5421
49.81
7772
1.18
E
-
02
604.2
0.5616
47.82
7753
1.10
E
-
02
650
0.5905
45.8
7726
1.00
E
-
02
800
0.685
39.2
7726
7.41
E
-
03
1000
1.169
30
7726
3.32
E
-
03Table 3
Thermal diffusivity functions considered in the current study.
Material
Thermal diffusivity functionα(T), m2/s
R
2
Temperature range  (K)
Stainless steelAISI 304
α
(
T
)
=
a
T
+
b   a=2.0×10-6, b=0.0037
0.9968
100–2500
Mild steel
α
(
T
)
=
a
T
2
+
b
T
+
c    a=1.0×10-8, b=-3.0×10-5, c=0.0276
0.9957
100–1000Figure 1
Thermal diffusivity of stainless steel AISI 304.Figure 2
Thermal diffusivity of mild steel.


## 3. The Method for Generating Approximate Analytic Solutions

The main approach for finding approximate analytic solutions consists of implementing the steps highlighted below.Given an IBVP of a PDE of the form(5)F(x,t,T(x,t),Tx,Tt,Txx,Txt,Ttt)=0.Step 1.
Reduction of IVBP of PDE to a BVP of ODE: the Lie symmetries of PDE (5) can be found by employing the standard well-known procedure, compare [7–9]. It is further required to systematically identify the symmetry that leaves the boundaries and boundary conditions invariant [14]. This is done by taking the most general symmetry operator X of PDE (5) and finding the conditions under which it leaves the boundaries invariant as well as imposing the restrictions on X due to invariance of boundary conditions on the boundary. Further details about these steps are illustrated in Section 4 below. The similarity variables of the symmetry that leaves the whole IBVP invariant will lead to the reduction of IBVP of PDE (5) to a BVP of ODE of the form(6)G(z,V(z),Vz,Vzz)=0.The aim of the remaining steps is to find an approximate solutionVApprox(z) of BVP of ODE (6) in function form and then use the similarity transformations of the symmetry to obtain approximate solution T(x,t) of the IBVP of PDE (5).Step 2.
Find numerical solutionVNum of BVP of ODE (6) and use this as a benchmark for obtaining function form VApprox(z) of the approximate solution of BVP of ODE (6).Step 3.
Obtain an initial guess
V
Initial
for the approximate analytic solution
V
Approx
(
z
)
. This is a crucial step as it will be used as a basis to generate approximate analytic solution in function form. Here we use homotopy perturbation method to obtain initial guess VInitial but other kinds of approximations like upper/lower solutions can also be employed.Step 4.
Improve the initial approximation
V
Initial
to get the approximate analytic solution
V
Approx
(
z
)
up to the desired level of accuracy. The improvement in the accuracy of initial approximation VInitial depends on reducing the difference between VInitial and the benchmark numerical solution curve VNum of BVP of ODE (6). In case of good initial approximation, as obtained in the next sections through homotopy perturbation method, the simulation based techniques similar to noise reduction or smoothing techniques of image processing should be sufficient. This idea is implemented in two ways in subsequent sections to improve the level of accuracy of approximate analytic solutions. In Section 4, the simulations based on perturbing the initial approximation VInitial are utilized to generate accurate enough approximate analytic solution VApprox(z) of the BVP of ODE, whereas in Section 5, the error curve obtained through the difference of VInitial and VNum is treated to generate the approximate analytic solution VApprox(z) of the BVP of ODE.Step 5.
Use the similarity variables inVApprox(z) to get the approximate analytic solution T(x,t) of the IBVP of PDE (5).Step 6.
Validate by carrying out a comparative analysis of the approximate analytic solutionT(x,t) at different times with the numerical solution the IBVP of PDE at corresponding times.
In Sections4 and 5, we illustrate implementation of the above procedure and provide simulation results and approximate analytic solutions for nonlinear heat conduction in stainless steel and mild steel, respectively.


## 4. Approximate Analytic Solution to IBVP for Stainless Steel AISI 304

As first test problem, we consider transient heat conduction in a semi-infinite solid bar made of AISI 304 stainless steel which is initially at temperatureTi=300°K and is subjected to a surface temperature Ts=900°K at x=0. The temperature-dependent diffusivity for stainless steel AISI 304, as estimated in Section 2, is given by
(7)α(T)=aT+b,wherea=2.0×10-6,b=0.0037.
The objective is to determine an approximate analytic expression for the temperature distribution T(x,t) satisfying the following IBVP:
(8)∂T(x,t)∂t=∂∂x((aT(x,t)+b)∂T(x,t)∂x),(9)T|t=0=Ti,T|x=0=Ts,T|x=∞=Ti,
where a=2.0×10-6 and b=0.0037.It is first required to systematically find the symmetry that preserves the IBVP ((8), (9)). Applying the Lie symmetry theory [7–9], the Lie symmetry algebra of the PDE (8) is determined to be four-dimensional and is generated by
(10)X1=∂∂x,X2=∂∂t,X3=x∂∂x+2t∂∂t,X4=ax2∂∂x+(aT+b)∂∂T.
In order to obtain the symmetry that leaves the whole IBVP invariant, we consider the general symmetry operator
(11)X=k1X1+k2X2+k3X3+k4X4
of PDE (8) and search for the operator that preserves the boundary and the boundary conditions (9).The invariance of the boundariesx=0, t=0 or equivalently
(12)[X(x-0)]x=0=0,[X(t-0)]t=0=0
implies
(13)k1=k2=0,
so X must be
(14)X=k3X3+k4X4.
In addition to the restrictions imposed by (13), the invariance of initial and boundary conditions, that is,
(15)[X(T-Ti)]t=0=0,on⁡T=Ti,[X(T-Ts)]x=0=0,on⁡T=Ts,
implies that we must have
(16)k4=0.
Hence the IBVP (8), (9) is invariant under the symmetry
(17)X=x∂∂x+2t∂∂t,
where we have chosen k3=1. The similarity variables
(18)z(x,t)=xt,V(z)=T
of X imply that the solution of IBVP (8), (9) is of the form T=V(z), where V(z) satisfies the ODE
(19)(aV+b)d2Vdz2+a(dVdz)2+z2dVdz=0,
with the boundary conditions
(20)V(z=0)=Ts=900,V(z=∞)=Ti=300.Next the reduced BVP (19), (20) is to be simultaneously solved numerically and by employing homotopy perturbation method. The homotopy solution will be utilized later as initial guess for the approximate analytic solution of BVP (19), (20) whereas numerical solution will be used as the benchmark curve to improve the accuracy of the initially guessed homotopy solution.In order to apply homotopy perturbation method [11–13] for finding approximate closed form solution of the boundary value problem (19), (20), we construct a homotopy of the expression (19) in the following form:
(21)bd2Vdz2+z2dVdz+p(aVd2Vdz2+a(dVdz)2)=0.
The homotopy parameter p has values 0 or 1; the value p=0 corresponds to a linear equation which can easily be solved and the value p=1 corresponds to the original equation (19). The solution V(z) can be expanded in terms of the homotopy parameter p as follows:
(22)V=V0+pV1+p2V2+⋯.
Substituting (22) in (21) and equating coefficients of the same powers of p, we obtain a set of equations from which the first two consist of the equation
(23)bd2V0dz2+z2dV0dz=0
with the boundary conditions
(24)V0(z=0)=Ts,V0(z⟶∞)=Ti
and the equation
(25)bd2V1dz2+z2dV1dz+aV0d2V0dz2+a(dV0dz)2=0
with the boundary conditions
(26)V1(0)=0,V1(z⟶∞)=0.
The solution of the boundary value problem (23), (24) can easily be obtained as
(27)V0(z)=Ts-(Ts-Ti)erf⁡(z2b).
Substituting solution (27) in (25) and removing the secular terms, we obtain
(28)bd2V1dz2+z2dV1dz+ae-z2/2b(Ts-Ti)2bπ+zae-z2/2b(Ts-Ti)Ts2b3π=0.
The solution of the boundary value problem (28) and (26) can be written in the following form:
(29)V1(z)=a(Ts-Ti)2b×{(1-erf⁡(z2b))zb23e-z2/2bTsπ+(1-erf⁡(z2b))erf⁡(z2b)(Ts-Ti)}.
Since we are looking for an approximate closed form solution of the problem (19), setting p=1 in (22), we obtain a first order approximate solution
(30)V(z)=V0(z)+V1(z).
Substituting (27) and (29) in (30), we obtain an approximate closed form solution of BVP (19), (20) in the following form:
(31)V(z)=Ts-(Ts-Ti)erf⁡(z2b)+a(Ts-Ti)2b×{(1-erf⁡(z2b))zb23e-z2/2bTsπ+(1-erf⁡(z2b))erf⁡(z2b)(Ts-Ti)}.
We call this solution the initial guess VInitial for the intended approximate analytic solution of the BVP (19), (20).The numerical solutionVNum of BVP (19), (20) is obtained using finite element method. Defining the initial relative error as
(32)EInitial=VInitial-VNumVNum,
it is found that
(33)Max⁡|EInitial|=0.06119845971.
Figure 3 gives the plots of the initial approximation VInitial and the numerical solution VNum, while the initial relative error EInitial is shown in Figure 5.Figure 3
Initial approximationVInitial and numerical solution VNum of BVP (19), (20).In order to improve the accuracy level of the initial approximation (31), we introduce the parameters α, β, δ, ϵ, and λ in (31) to obtain the expression
(34)V~(z)=Ts-(Ts-Ti)erf⁡(ϵz2b)+a(Ts-Ti)2b×{(1-erf⁡(δz2b))βzb23e-(αz)2/2bTsπ+(1-erf⁡(δz2b))erf⁡(λz2b)(Ts-Ti)}.
Small variations of the parameters from the values
(35)α=1,β=1,δ=1,ϵ=1,λ=1
in V~(z) generate a sequence of curves in the neighborhood of the initial approximation curve Vinitial. Several careful numerical simulations are performed for small variation of the above parameters, leading to the values
(36)α=1.4,β=1,δ=0.93,ϵ=0.9101,λ=0.91
for which the relative error between V~(z) and VNum is less than 0.4%.Hence we obtain an approximate analytic solutionVApprox of BVP (19), (20) as
(37)VApprox(z)=Ts-(Ts-Ti)erf⁡(ϵz2b)+a(Ts-Ti)2b×{(1-erf⁡(δz2b))βzb23e-(αz)2/2bTsπ+(1-erf⁡(δz2b))erf⁡(λz2b)(Ts-Ti)},
where the parameters α, β, δ, ϵ, and λ are given by (36), with
(38)Max⁡|VApprox-VNumVNum|=0.003959810186.
The plots of the numerical solution VNum and the approximate analytic solution VApprox are given in Figure 4, and Figure 5 shows the plot of Error(z) where
(39)Error(z)=VApprox-VNumVNum.Figure 4
Approximate analytic solutionVApprox and numerical solution VNum of BVP (19), (20).Figure 5
Comparison of relative errors for the initial approximation and approximate analytic solution of BVP (19), (20).Finally the similarity variables (18) provide the approximate analytic solution of IBVP (8), (9) as
(40)T(x,t)=Ts-(Ts-Ti)erf⁡(ϵ(x/t)2b)+a(Ts-Ti)2b×{β(x/t)b23e-(α(x/t))2/2bTsπ+(1-erf⁡(δ(x/t)2b))×erf⁡(λ(x/t)2b)(Ts-Ti)β(x/t)b23e-(α(x/t))2/2bTsπ},
where the parameters α, β, δ, ϵ, and λ are given by (36). A surface plot of the approximate analytic solution T(x,t) is provided in Figure 6.Figure 6
Approximate analytic solutionT(x,t) of the IBVP (8), (9).For validation purpose, we carry out a comparative analysis of the approximate analytic solutionT(x,t) and the numerical solution of IBVP (8), (9). In order to numerically solve the IBVP (8), (9), a transient thermal analysis is performed using FEA software ANSYS. The given structure is modeled using 3D Conduction Bar Elements (LINK33). LINK33 is a uniaxial element with the ability to conduct heat between its nodes. The element has a single degree of freedom, temperature, at each node point. The conducting bar is applicable to a steady-state or transient thermal analysis. A refined uniform mesh is used to model the nonlinear thermal gradient through the solid. The length of the model is taken as L=2 m assuming that no significant temperature change occurs at the interior end point during the time period of interest. This assumption is validated by the temperature of node at x=L at the end of the transient analysis.Figure7 shows the comparison of temperature distribution at different times using the approximate analytic solution (40) and the numerical solutions at corresponding times. The figure shows good agreement between approximate analytic solution and the numerical results.Figure 7
Comparison of approximate analytic solutionT(x,t) and numerical solution of IBVP (8), (9) at different times.


## 5. Approximate Analytic Solution to IBVP for Mild Steel

In this test problem, we consider transient heat conduction in a semi-infinite solid bar made of mild steel which is initially at temperatureTi=300°K and is subjected to a surface temperature Ts=900°K at x=0. The temperature-dependent diffusivity for mild steel, as estimated in Section 2, is given by
(41)α(T)=aT2+bT+c,wherea=1.0×10-8,b=-3.0×10-5,c=0.0276.
The objective is to determine an approximate analytic expression for the temperature distribution T(x,t) satisfying the following IBVP:
(42)∂T(x,t)∂t=∂∂x((aT(x,t)2+bT(x,t)+c)∂T(x,t)∂x),(43)T|t=0=Ti,T|x=0=Ts,T|x=∞=Ti,
where a=1.0×10-8,b=-3.0×10-5, and b=0.0276.In order to identify the symmetry that preserves the IBVP (42), (43) we first find the Lie symmetry algebra of the PDE (42) which is determined to be three-dimensional and is generated by
(44)X1=∂∂x,X2=∂∂t,X3=x∂∂x+2t∂∂t.
Using the restrictions imposed by the boundaries and boundary conditions in the manner similar to Section 4, we find that the symmetry that leaves the whole IBVP (42), (43) invariant is
(45)X=x∂∂x+2t∂∂t.
The similarity variables
(46)z(x,t)=xt,V(z)=T
of X imply that the solution of IBVP (42), (43) is of the form T=V(z), where V(z) satisfies the ODE
(47)(aV2+bV+c)d2Vdz2+(2aV+b)(dVdz)2+z2dVdz=0
with the boundary conditions
(48)V(z=0)=Ts=900,V(z=∞)=Ti=300.
Next the reduced BVP (47), (48) is simultaneously solved numerically and by homotopy perturbation method. The homotopy solution will be utilized below as initial guess VInitial for the approximate analytic solution VApprox of BVP (47), (48) whereas the numerical solution VNum will be used as the benchmark curve to improve the accuracy of the initially guessed homotopy solution. In order to find homotopy solution VInitial of the boundary value problem (47), (48), we construct a homotopy of the expression (47) in the following form:
(49)cd2Vdz2+z2dVdz+p[aV2d2Vdz2+bVd2Vdz2+2aV(dVdz)2+b(dVdz)2]=0,
and apply the same procedure as in the previous section to obtain an approximate closed form solution of the BVP (47), (48) as
(50)V(z)=Ts-(Ts-Ti)erf⁡(z2c)+e-z2/2c(Ts-Ti)12cπ×{(z2c)2-6bTs+ez2/2cerf⁡c(z2c)×(4aπTs2-4aπTs2erf⁡c(z2c)2-6bπTierf⁡(z2c)+6bTs+9bπTserf⁡(z2c)-4aπTsTi+4aπTsTierf⁡c(z2c)2)}.
We call this solution the initial guess VInitial for the intended approximate analytic solution of the BVP (47), (48).The numerical solutionVNum of BVP (47), (48) is obtained using finite element method. Defining the initial relative error as
(51)EInitial=VInitial-VNumVNum,
it is found that
(52)Max⁡|EInitial|=0.1500655840.
Figure 8 gives the plots of the initial approximation VInitial and the numerical solution VNum, while the initial relative error EInitial is shown in Figure 11.Figure 8
Initial approximationVInitial and numerical solution VNum of BVP (47), (48).In order to improve the accuracy level of the initial approximation (31), we adopt an approach different from the previous section and focus on approximating the difference between VNum and VInitial in function form. Plotting the difference
(53)Δ(z)=VInitial-VNum,
as shown in Figure 9, implies that the difference Δ(z) can be represented by a skewed function. As an ansatz for approximating Δ(z) we take the skewed function of the form
(54)f(z)=βe-Kz2erf⁡(αz),
where the parameter α mainly contributes to the location of the peak, the parameter K mainly affects the shrinking or expanding of the skewed curve, and β controls the height of the peak of the curve. Performing several careful numerical simulations suggest the values of the parameters
(55)α=0.0105,K=0.000014,β=101.02
giving accurate enough function form representation of Δ(z) that makes the relative error between approximate analytic solution and numerical solution about 0.7%, as shown below.Figure 9
The differenceΔ between the initial approximation VInitial and numerical solution VNum of BVP (47), (48).Defining(56)VApprox(z)=VInitial-βe-Kz2erf⁡(αz),
where the parameters α, β, and K are given by (55), provides an approximate analytic solution VApprox of BVP (47), (48) with
(57)Max⁡|VApprox-VNumVNum|=0.00740878478.
The plots of the numerical solution VNum and the approximate analytic solution VApprox are given in Figure 10, and Figure 11 shows the plot of Error(z) where
(58)Error(z)=VApprox-VNumVNum.Figure 10
Approximate analytic solutionVApprox and numerical solution VNum of BVP (47), (48).Figure 11
Comparison of relative errors for the initial approximation and approximate analytic solution of BVP (47), (48).Finally the similarity variables (46) provide the approximate analytic solution of IBVP (42), (43) as
(59)T(x,t)=Ts-(Ts-Ti)erf⁡(x/t2c)+e-(x/t)2/2c(Ts-Ti)12cπ×{+4aπTsTierf⁡c(x/t2c)2)-6bTs+ez2/2cerf⁡c(x/t2c)×(4aπTs2-4aπTs2erf⁡c(x/t2c)2-6bπTierf⁡(x/t2c)+6bTs+9bπTserf⁡(x/t2c)-4aπTsTi+4aπTsTierf⁡c(x/t2c)2)}-βe-Kz2erf⁡(αz),
where the parameters α, β, and K are given by (55). A surface plot of the approximate analytic solution T(x,t) is provided in Figure 12.Figure 12
Approximate analytic solutionT(x,t) of the IBVP (42) and (43).For validation purpose, we carry out a comparative analysis of the approximate analytic solutionT(x,t) and the numerical solution of IBVP (42), (43). The numerical solutions of IBVP (42), (43) are obtained by performing a transient thermal analysis using FEA software ANSYS, as explained in the previous section. Figure 13 shows the comparison of temperature distribution at different times using the approximate analytic solution (59) and the numerical solutions at corresponding times. The figure shows good agreement between approximate analytic solution and the numerical results.Figure 13
Comparison of approximate analytic solutionT(x,t) and numerical solution of IBVP (42) and (43) at different times.


## 6. Conclusion

The mathematical modeling of most of the physical processes in fields like diffusion, chemical kinetics, fluid mechanics, wave mechanics, and general transport problems is governed by such nonlinear PDEs whose exact analytic solutions are hard to find. Therefore, the methods of finding approximate analytic solutions become important in investigating such nonlinear PDEs. In this paper, we present a new systematic procedure for generating approximate analytic solutions of nonlinear PDEs, particularly arising from heat conduction problems. The methodology followed here exploits an effective combination of Lie symmetry analysis, homotopy perturbation method, finite element method, and simulation based error reduction techniques. The similarity variables of an appropriate Lie symmetry are first utilized to reduce the IBVP of PDE to a BVP of ODE. The reduced ODE problem is then simultaneously solved numerically and by homotopy perturbation method to, respectively, generate numerical solution and initial approximation of the approximate analytic solution of BVP of ODE. Next, using the numerically generated solution curve as the accuracy benchmark, an error reduction of the initial approximation is carried out by utilizing simulations to generate an approximate analytic solution which meets the accuracy benchmark of numerical solution. Finally the similarity transformations provide the approximate analytic solution of the IBVP of the original PDE. The approach is applied to obtain approximate analytic solutions for test problems consisting of transient heat conduction in bars made of stainless steel AISI 304 and mild steel. The validity of solutions is verified by a comparison between the approximate analytic solutions and the numerical solutions of the test problems. The results indicate good agreement between the approximate analytical solutions and the corresponding numerical solutions. The idea presented here to get approximate analytic solution of PDE, that is, approximating the surfaceT(x,t), practically involves approximating a curve V(z) which is a tractable problem in comparison to increasingly complex and intractable problem of approximating the surface T(x,t) itself. It can become particularly efficient when a reasonably accurate initial approximation of the solution curve V(z) of the reduced ODE can be obtained, as was the case here in the form of homotopy solutions. Another advantage of the approximate analytic solutions obtained by our approach is that these can be used in real time or negligible CPU time to evaluate the solution of the PDE with the same accuracy as the numerical solution.Although, here, we have restricted applications of our approach to heat conduction problems, it can be adapted for obtaining approximate analytic solutions for the class of PDEs where the PDE can be reduced to an ODE through similarity variables. For instance, the approach can be directly applied to all the reduced-via-similarity BVPs of ODEs in [15]. For further application, the approach can be extended to obtain approximate solutions where the reduction is a system of ODEs like the reduced laminar boundary layer problems in [16, 17].



---

*Source: 423421-2014-07-13.xml*