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

and Applied Analysis 3


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 with three physical properties or the equation with just one coefficient is solved to predict the temperature distribution during the heating or cooling process.Here  is the thermal conductivity,  is the density,   is the heat capacity, and  = /  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][2][3][4][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][7][8][9][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.Section 3 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][8][9][11][12][13], respectively, for an introduction to Lie symmetry method and homotopy perturbation method.

Formulation of Test Problems
The nonlinear PDE that describes the transient nonlinear heat conduction in a one-dimensional medium is where (, ) denotes the temperature at a point  at time  and () 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 () 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 secondorder 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 ( 2 ).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  with the thermal diffusivity functions of stainless steel and mild steel, respectively.

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  (, ,  (, ) ,   ,   ,   ,   ,   ) = 0. (5) 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][8][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  of PDE ( 5) and finding the conditions under which it leaves the boundaries invariant as well as imposing the restrictions on  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 The aim of the remaining steps is to find an approximate solution  Approx () of BVP of ODE (6) in function form and then use the similarity transformations of the symmetry to obtain approximate solution (, ) of the IBVP of PDE (5).
Step 2. Find numerical solution  Num of BVP of ODE (6) and use this as a benchmark for obtaining function form  Approx () of the approximate solution of BVP of ODE (6).
Step 3. Obtain an initial guess  Initial for the approximate analytic solution  Approx ().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  Initial but other kinds of approximations like upper/lower solutions can also be employed.
Step 4. Improve the initial approximation  Initial to get the approximate analytic solution  Approx () up to the desired level of accuracy.The improvement in the accuracy of initial approximation  Initial depends on reducing the difference between  Initial and the benchmark numerical solution curve  Num 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  Initial are utilized to generate accurate enough approximate analytic solution  Approx () of the BVP of ODE, whereas in Section 5, the error curve obtained through the difference of  Initial and  Num is treated to generate the approximate analytic solution  Approx () of the BVP of ODE.
Step 5. Use the similarity variables in  Approx () to get the approximate analytic solution (, ) of the IBVP of PDE (5).
Step 6. Validate by carrying out a comparative analysis of the approximate analytic solution (, ) at different times with the numerical solution the IBVP of PDE at corresponding times.
In Sections 4 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.

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 temperature   = 300 ∘ K and is subjected to a surface temperature   = 900 ∘ K at  = 0.The temperature-dependent diffusivity for stainless steel AISI 304, as estimated in Section 2, is given by  () =  + , where  = 2.0 × 10 −6 ,  = 0.0037.(7) The objective is to determine an approximate analytic expression for the temperature distribution (, ) satisfying the following IBVP: where  = 2.0 × 10 −6 and  = 0.0037.
It is first required to systematically find the symmetry that preserves the IBVP (( 8), ( 9)).Applying the Lie symmetry theory [7][8][9], the Lie symmetry algebra of the PDE ( 8) is determined to be four-dimensional and is generated by In order to obtain the symmetry that leaves the whole IBVP invariant, we consider the general symmetry operator of PDE (8) and search for the operator that preserves the boundary and the boundary conditions (9).
The invariance of the boundaries  = 0,  = 0 or equivalently implies so  must be In addition to the restrictions imposed by (13), the invariance of initial and boundary conditions, that is, implies that we must have Hence the IBVP (8), ( 9) is invariant under the symmetry where we have chosen  3 = 1.The similarity variables of  imply that the solution of IBVP ( 8), ( 9) is of the form  = (), where () satisfies the ODE with the boundary conditions 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][12][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: The homotopy parameter  has values 0 or 1; the value  = 0 corresponds to a linear equation which can easily be solved and the value  = 1 corresponds to the original equation ( 19).The solution () can be expanded in terms of the homotopy parameter  as follows: Substituting ( 22) in (21) and equating coefficients of the same powers of , we obtain a set of equations from which the first two consist of the equation with the boundary conditions and the equation with the boundary conditions The solution of the boundary value problem (23), (24) can easily be obtained as Substituting solution (27) in (25) and removing the secular terms, we obtain The solution of the boundary value problem ( 28) and ( 26) can be written in the following form: Substituting ( 27) and ( 29) in (30), we obtain an approximate closed form solution of BVP ( 19), (20) in the following form: We call this solution the initial guess  Initial for the intended approximate analytic solution of the BVP (19), (20).The numerical solution  Num of BVP ( 19), (20) is obtained using finite element method.Defining the initial relative error as it is found that Max      Initial     = 0.06119845971.
Figure 3 gives the plots of the initial approximation  Initial and the numerical solution  Num , while the initial relative error  Initial is shown in Figure 5.
In order to improve the accuracy level of the initial approximation (31), we introduce the parameters , , , , and  in (31) to obtain the expression Small variations of the parameters from the values in Ṽ() generate a sequence of curves in the neighborhood of the initial approximation curve  initial .Several careful numerical simulations are performed for small variation of the above parameters, leading to the values for which the relative error between Ṽ() and  Num is less than 0.4%.
Hence we obtain an approximate analytic solution  Approx of BVP ( 19), (20) as where the parameters , , , , and  are given by (36), with The plots of the numerical solution  Num and the approximate analytic solution  Approx are given in Figure 4, and Figure 5 shows the plot of Error() where  Finally the similarity variables (18) provide the approximate analytic solution of IBVP ( 8 8), (9).
where the parameters , , , , and  are given by (36).A surface plot of the approximate analytic solution (, ) is provided in Figure 6.
For validation purpose, we carry out a comparative analysis of the approximate analytic solution (, ) 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  = 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  =  at the end of the transient analysis.
Figure 7 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.

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 temperature   = 300 ∘ K and is subjected to a surface temperature   = 900 ∘ K at  = 0.The temperaturedependent diffusivity for mild steel, as estimated in Section 2, is given by where  = 1.0 × 10 −8 ,  = −3.0× 10 −5 ,  = 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 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 The similarity variables of  imply that the solution of IBVP (42), ( 43) is of the form  = (), where () satisfies the ODE with the boundary conditions 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  Initial for the approximate analytic solution  Approx of BVP (47), (48) whereas the numerical solution  Num will be used as the benchmark curve to improve the accuracy of the initially guessed homotopy solution.In order to find homotopy solution  Initial of the boundary value problem (47), (48), we construct a homotopy of the expression (47) in the following form: and apply the same procedure as in the previous section to obtain an approximate closed form solution of the BVP (47), (48) as We call this solution the initial guess  Initial for the intended approximate analytic solution of the BVP (47), (48).The numerical solution  Num of BVP (47), (48) is obtained using finite element method.Defining the initial relative error as it is found that Max      Initial     = 0.1500655840.
Figure 8 gives the plots of the initial approximation  Initial and the numerical solution  Num , while the initial relative error  Initial is shown in Figure 11.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  Num and  Initial in function form.Plotting the difference as shown in Figure 9, implies that the difference Δ() can be represented by a skewed function.As an ansatz for approximating Δ() we take the skewed function of the form where the parameter  mainly contributes to the location of the peak, the parameter  mainly affects the shrinking or expanding of the skewed curve, and  controls the height of the peak of the curve.Finally the similarity variables (46) provide the approximate analytic solution of IBVP (42), (43) as where the parameters , , and  are given by (55).A surface plot of the approximate analytic solution (, ) is provided in Figure 12.
For validation purpose, we carry out a comparative analysis of the approximate analytic solution (, ) 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.

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 surface (, ), practically involves approximating a curve () which is a tractable problem in comparison to increasingly complex and intractable problem of approximating the surface (, ) itself.It can become particularly efficient when a reasonably accurate initial approximation of the solution curve () 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].

Table 1 :
Thermal properties of AISI 304 stainless steel.

Table 2 :
Thermal properties of mild steel.

Table 3 :
Thermal diffusivity functions considered in the current study.
The plots of the numerical solution  Num and the approximate analytic solution  Approx are given in Figure10, and Figure11shows the plot of Error() where