Mathematical Methods and Models in the Natural to the Life Sciences
View this Special IssueResearch Article  Open Access
Approximate Analytic Solutions of Transient Nonlinear Heat Conduction with TemperatureDependent Thermal Diffusivity
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 semiinfinite 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.
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 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 materialspecific 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 temperaturedependent. 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 finitedifference method (FDM), the finiteelement method (FEM), and the boundaryelement 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 realtime 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 reallife 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 initialboundary 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–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 onedimensional medium is where denotes the temperature at a point at time and is the temperaturedependent 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 temperaturedependent 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 (). Using these thermal diffusivity functions, we will apply our method to investigate transient heat conduction in semiinfinite solid. Precisely, in Sections 4 and 5, approximate analytic solutions will be constructed for initialboundary value problem of the form with the thermal diffusivity functions of stainless steel and mild steel, respectively.



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
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 wellknown 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 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 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 of BVP of ODE (6) and use this as a benchmark for obtaining function form of the approximate solution of BVP of ODE (6).
Step 3. Obtain an initial guess for the approximate analytic solution . 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 but other kinds of approximations like upper/lower solutions can also be employed.
Step 4. Improve the initial approximation to get the approximate analytic solution up to the desired level of accuracy. The improvement in the accuracy of initial approximation depends on reducing the difference between and the benchmark numerical solution curve 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 are utilized to generate accurate enough approximate analytic solution of the BVP of ODE, whereas in Section 5, the error curve obtained through the difference of and is treated to generate the approximate analytic solution of the BVP of ODE.
Step 5. Use the similarity variables in 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.
4. Approximate Analytic Solution to IBVP for Stainless Steel AISI 304
As first test problem, we consider transient heat conduction in a semiinfinite solid bar made of AISI 304 stainless steel which is initially at temperature 00°K and is subjected to a surface temperature 00°K at . The temperaturedependent diffusivity for stainless steel AISI 304, as estimated in Section 2, is given by The objective is to determine an approximate analytic expression for the temperature distribution satisfying the following IBVP: where and .
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 fourdimensional 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 , 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 . 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–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 or ; the value corresponds to a linear equation which can easily be solved and the value 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: Since we are looking for an approximate closed form solution of the problem (19), setting in (22), we obtain a first order approximate solution 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 for the intended approximate analytic solution of the BVP (19), (20).
The numerical solution of BVP (19), (20) is obtained using finite element method. Defining the initial relative error as it is found that Figure 3 gives the plots of the initial approximation and the numerical solution , while the initial relative error 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 . Several careful numerical simulations are performed for small variation of the above parameters, leading to the values for which the relative error between and is less than .
Hence we obtain an approximate analytic solution of BVP (19), (20) as where the parameters , , , , and are given by (36), with The plots of the numerical solution and the approximate analytic solution are given in Figure 4, and Figure 5 shows the plot of where
Finally the similarity variables (18) provide the approximate analytic solution of IBVP (8), (9) as 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 steadystate 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 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.
5. Approximate Analytic Solution to IBVP for Mild Steel
In this test problem, we consider transient heat conduction in a semiinfinite solid bar made of mild steel which is initially at temperature 0°K and is subjected to a surface temperature 0°K at . The temperaturedependent diffusivity for mild steel, as estimated in Section 2, is given by The objective is to determine an approximate analytic expression for the temperature distribution satisfying the following IBVP: where , and .
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 threedimensional 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 for the approximate analytic solution of BVP (47), (48) whereas the numerical solution will be used as the benchmark curve to improve the accuracy of the initially guessed homotopy solution. In order to find homotopy solution 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 for the intended approximate analytic solution of the BVP (47), (48).
The numerical solution of BVP (47), (48) is obtained using finite element method. Defining the initial relative error as it is found that Figure 8 gives the plots of the initial approximation and the numerical solution , while the initial relative error 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 and 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. Performing several careful numerical simulations suggest the values of the parameters giving accurate enough function form representation of that makes the relative error between approximate analytic solution and numerical solution about , as shown below.
Defining where the parameters , , and are given by (55), provides an approximate analytic solution of BVP (47), (48) with The plots of the numerical solution and the approximate analytic solution are given in Figure 10, and Figure 11 shows the plot of where
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.
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 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 reducedviasimilarity 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].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
References
 K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equations: An Introduction, Cambridge University Press, New York, NY, USA, 1994. View at: MathSciNet
 J. Blobner, R. A. Białecki, and G. Kuhn, “Boundaryelement solution of coupled heat conductionradiation problems in the presence of shadow zones,” Numerical Heat Transfer B: Fundamentals, vol. 39, no. 5, pp. 451–478, 2001. View at: Publisher Site  Google Scholar
 B. Chen, Y. Gu, Z. Guan, and H. Zhang, “Nonlinear transient heat conduction analysis with precise time integration method,” Numerical Heat Transfer B, vol. 40, no. 4, pp. 325–341, 2001. View at: Publisher Site  Google Scholar
 F. Asllanaj, A. Milandri, G. Jeandel, and J. R. Roche, “A finite difference solution of nonlinear systems of radiativeconductive heat transfer equations,” International Journal for Numerical Methods in Engineering, vol. 54, no. 11, pp. 1649–1668, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 A. Fic, R. A. Białecki, and A. J. Kassab, “Solving transient nonlinear heat conduction problems by proper orthogonal decomposition and the finiteelement method,” Numerical Heat Transfer B: Fundamentals, vol. 48, no. 2, pp. 103–124, 2005. View at: Publisher Site  Google Scholar
 W. F. Ames, Nonlinear Partial Differential Equations in Engineering, vol. 1, Academic Press, New York, NY, USA, 1965.
 G. W. Bluman and S. Anco, Symmetry and Integration Methods for Differential Equations, Springer, New York, NY, USA, 2002. View at: MathSciNet
 N. H. Ibragimov, Elementary Lie Group Analysis and Ordinary Differential Equations, vol. 4, John Wiley & Sons, Chichester, UK, 1999. View at: MathSciNet
 P. J. Olver, Applications of Lie Groups to Differential Equations, Springer, New York, NY, USA, 1986.
 L. V. Ovsiannikov, Group Analysis of Differential Equations, Academic Press, New York, NY, USA, 1982. View at: MathSciNet
 J. H. He, “Homotopy perturbation technique,” Computer Methods in Applied Mechanics and Engineering, vol. 178, no. 34, pp. 257–262, 1999. View at: Publisher Site  Google Scholar  MathSciNet
 J. He, “Bookkeeping parameter in perturbation methods,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 2, no. 3, pp. 257–264, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J.H. He, “Addendum: new interpretation of homotopy perturbation method,” International Journal of Modern Physics B, vol. 20, no. 18, pp. 2561–2568, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 H. Azad, M. T. Mustafa, and A. F. M. Arif, “Analytic solutions of initialboundaryvalue problems of transient conduction using symmetries,” Applied Mathematics and Computation, vol. 215, no. 12, pp. 4132–4140, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 L. Dresner, Similarity Solutions of Nonlinear Partial Differential Equations, vol. 88 of Research Notes in Mathematics Series, Longman, 1983.
 A. Aziz, “A similarity solution for laminar thermal boundary layer over a flat plate with a convective surface boundary condition,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 4, pp. 1064–1068, 2009. View at: Publisher Site  Google Scholar
 S. Mukhopadhyay, K. Bhattacharyya, and G. C. Layek, “Steady boundary layer flow and heat transfer over a porous moving plate in presence of thermal radiation,” International Journal of Heat and Mass Transfer, vol. 54, no. 1314, pp. 2751–2757, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
Copyright
Copyright © 2014 M. T. Mustafa et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.