Received Date: April 20, 2012; Accepted Date: June 20, 2012; Published Date: June 23, 2012
Citation: Nave O, Lehavi Y, Gol’dshtein V, Ajadi S (2012) Application of the Homotopy Perturbation Method (HPM) and the Homotopy Analysis Method (HAM) to the Problem of the Thermal Explosion in a Radiation Gas with Polydisperse Fuel Spray. J Appl Computat Math 1:115. doi: 10.4172/2168-9679.1000115
Copyright: © 2012 Nave O, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Applied & Computational Mathematics
The aim of this work is to apply the homotopy perturbation method and homotopy analysis method to the problem of thermal explosion in a flammable gas mixture with the addition of volatile fuel droplets. The system of equations that describes the effects of heating, evaporation, and combustion of fuel in a polydisperse spray is simplified. Both convective and radiative heating of droplets is taken into account in the model. The model for the radiative heating of droplets takes into account the semitransparency of the droplets. The results of the analysis have been applied to the modeling of the thermal explosion in diesel engines. We applied the Homotopy Perturbation Method and the Homotopy Analysis Method to the new model and we found the region of the convergence of the considered solutions of the relevant physical parameters. The results demonstrate that these methods are very effective for solving nonlinear problems in science and engineering.
Homotopy perturbation method (HPM); Homotopy analysis method (HAM); Nonlinear differential equations; Thermal explosion; Polydisperse fuel spray
• A pre-exponential rate factor ( s−1 )
• B universal gas constant ( Jkmol−1K−1 )
• C molar concentration ( kmolm−3 )
• c specific heat capacity ( Jkg−1K−1)
• E activation energy ( Jkmol−1 )
• F −M refer to the model (21)-(24) with b = 0
• k number of droplets
• L liquid evaporation energy (i.e., latent heat of evaporation, Enthalpy of evaporation) ( Jkg−1 )
• m droplet mass
• n number of droplets per unit volume (m−3 )
• Q combustion energy ( Jkg−1 )
• q heat flux (Wm−2 )
• R radius of droplet (m )
• r dimensionless radius
• T temperature ( K )
• t time ( s )
• treact characteristic reaction time ( s ) defined in Equation (25)
• α dimensionless volumetric phase content
• β dimensionless reduced initial temperature (with respect to the so-called activation temperature E / B )
• γ dimensionless parameter that represents the reciprocal of the final dimensionless adiabatic temperature of the thermally insulated system after the explosion has been completed
• ε1i i =1,..., k dimensionless parameters defined in Equation (25) and describes the competition between the combustion and the evaporation processes
• ε2i i =1,..., k dimensionless parameters defined in Equation (25). This parameter relates the heat released during combustion and energy that is needed to evaporate all the fuel droplets
• ε3i i =1,..., k dimensionless parameters defined in Equation (25). This parameter describes the ratio of treact and the characteristic droplet heating time
• ε4i i =1,..., k dimensionless parameters defined in Equation (25). This parameter is proportional to the ratio of radiative and convective fluxes
• η dimensionless fuel concentration
• θ dimensionless temperature
• λ thermal conductivity ( Wm−1K−1 )
• μ molar mass ( kgkmol−1 )
• ν dimensionless parameter defined in Equation (25)
• ρ density ( kgm−3 )
• σ Stefan-Boltzmann constant (W m−2K−4 )
• τ dimensionless time
• ψ represents the internal characteristics of the fuel (the ratio of the specific combustion energy and the latent heat of evaporation) defined in Equation (25) and for diesel fuel ψ >>1
• con convection
• d liquid fuel droplets
• f combustible gas component of the mixture
• g gas mixture
• l liquid phase
• p under constant pressure
• rad adiation
• 0 initial state
In this paper we investigated the problem of thermal explosion in a fuel mixture and gas. This problem, in most cases, has been studied based on the application of computational fluid dynamics (CFD) packages . This method could take into account the complicated geometry of the enclosure and the chemistry of the processes. Hence, this makes it particularly attractive for engineering applications including the modeling of combustion processes in diesel engines. Other approaches to this problem are based on a asymptotic analysis of equations describing the limiting cases of the processes. One of these approaches is based on the application of the zero-order approximation of the geometric version of the asymptotic method of integral manifolds (MIM), developed for combustion applications in [2,3]. For example, in  the method of integral manifold was applied to a specific problem of modeling of ignition process in a diesel engines by using the P-1 model. The chemical term was presented in the Arrhenius form with the pre-exponential factor calculated from the enthalpy equation, using the well known Shell autoignition model. The results predicted by the analytical solution were compared with those that predicted by the computational fluid dynamics package VECTIS. The effects of the thermal radiation were shown to be very significant especially at high temperatures.
Other asymptotic methods were proposed by  in 1992, known as the homotopy analysis method (HAM) and by  in 1998, known as the homotopy perturbation method (HPM), which is a special case of HAM. The HPM and the HAM methods are mathematical tools that are based on homotopy, a fundamental concept in topology and differential geometry. They are analytical approaches to formulate the series solution of linear and nonlinear partial differential equations. We refer the reader to  for an enlightening comparison between HAM and HPM.
HPM couples the homotopy technology and perturbation method including the modified Lindstedt-Poincare method . The authors of  modified the Multiple Scales method by incorporating the time transformation of Lindstedt Poincare method. In  the authors contrasted two different approaches of Lindstedt-Poincare methods using the duffing equation. The main deficiencies in applying perturbation methods is that a small parameter is needed in the equations.
The HPM was further developed and improved and applied to nonlinear oscillators with discontinuities , nonlinear wave equations , boundary value problem , limit cycle and bifurcation of nonlinear problems  and many other subjects. In recent years, the application of the homotopy perturbation method (HPM) in nonlinear problems has been developed by scientists and engineers, because this method deforms the difficult problem under study into a simple problem which is easy to solve [15-17]. Most perturbation methods assume a small parameter exists, but most nonlinear problems have no small parameter at all. Unlike analytical perturbation methods, the HPM and HAM do not depend on a small parameter which is difficult to find.
These two methods also provide a simple way to ensure the convergence of the series solution. Moreover, these methods provide a large degree of freedom to choose an appropriate base functions to approximate the linear and nonlinear problems . Another important advantage of this method is that one can construct a continuous mapping of an initial guess approximation to the exact solution of the given problem through an auxiliary linear operator. To ensure the convergence of the series solution an auxiliary parameter is used. In  Liao has substantiated that the HAM differs from the other analytical methods in that it ensures the convergence of the series solution by choosing a proper value for the convergence-control parameter.
In this paper we have rewritten the model that was proposed by  for polydisperse fuel spray and applied the HPM and HAM to the problem of thermal explosion in a fuel mixture and gas. Based on these two methods, we present an analytical solutions for various values of the relevant physical parameter and we discuss the convergence of these solutions. We also compare our results to numerical solutions.
To explain this method, let us consider the following equation:
with the boundary conditions of:
where A , B , f (r) are a general differential operator, a boundary operator, a known analytical function respectively. Ω is the domain. Generally, the operator A can be decomposed into a linear part L and a nonlinear part N(u) . Hence, Equation (1) can be written as:
By the homotopy technique, we construct a homotopy which satisfies:
where is an embedding parameter and α (r, p) is a function of r and p , and u0(r) denote the initial approximation of α (r) .
When p = 0 we have
and when p =1 we have
As we mention before, L denote an auxiliary linear operator. In addition L have the property:
Using (7), it is clear that for p = 0
is the solution of the equation:
And for p =1
is the solution of the equation:
When the embedding parameter p increase from 0 to1, the solution α (r, p) of the equation:
depends upon the embedding parameter p and the varies from the initial approximation u0 (r) to the solution u(r) of equation (1). In topology, such a kind of continuous variation is called deformation.
According to the HPM, we can first use the embedding parameter p as a “small parameter”, and assume that the solutions of Equation (4) can be written as a power series in p :
Setting p =1 yields in the approximate solution of (1) to
The combination of the perturbation method and the homotopy method is called the HPM, which eliminates the drawbacks of the traditional perturbation methods while keeping all their advantages. The rate of convergent of series (13) depends on the nonlinear operator A(u) .
The physical assumptions are as follows: The combustible gas mixture contains evaporating ideal spherical droplets of fuel. The liquid droplets form a polydisperse spray. The medium is assumed to be spatially homogeneous. The variations in pressure in the enclosure, and their influence on the combustion processes are ignored. The heat flux from the burning gas to the droplets is assumed to consist of two components: convection and radiation, and the form of these two components is as follows:
where for μm units the value of and are: and . The energy that is needed for heating fuel vapor from the droplet temperature to gas temperature is ignored. The thermal conductivity of the liquid phase is much greater than that of the gas phase. The volume fraction of the liquid phase is much less than that of the gas phase. The heat transfer coefficient in the liquidgas mixture is assumed to be controlled by the thermal properties of the gas phase. External heat losses are ignored. Fuel droplets are semitransparent. Combustion takes place in the gas phase only. Combustion is modeled as a one-step first-order exothermic reaction with gaseous fuel as a deficient reactant. The droplets are assumed to be stationary. Under these assumptions, we rewrite the model as in , which is in the form of monodisperse fuel spray, to a polydisperse fuel spray as follows:
In non-dimensional parameters and by applying the Frank- Kamenetskii approximation , the model has the form of:
where the following dimensionless parameters have been introduced:
The non-dimensional initial conditions are:
For simplicity, in our numerical simulations and when applying the HPM and the HAM we assume b = 0 .
By applying the HPM method to the system of Equations (21)-(24) we obtain the following HPM-system:
with the initial conditions:
According to HPM method the terms α1 , α2, α3i and α4i for i =1,..., k has the form of:
Suppose that the solution of (27)-(30) takes the form of
Substituting Equations (32) with the initial conditions (31) into Equations (27)-(30) for three different size of droplets i.e., k = 3 , using the Taylor expansion for the exponent [23,24] and finally collecting the terms in power of p up to order 3 we obtain:
equations for gas temperature θg
equations for concentration η
We derived a system of (24) ordinary differential equation with 24 with unknown functions: α1,1 , α1,2 , α1,3 , α2,1 , α2,2 , α2,3 , α13,1, α13,2, α13,3, α23,1, α23,2, α23,3. α33,1, α33,2, α33,3, α14,1, α14,2, α14,3, α24,1, α24,2, α24,3, α34,1, α34,2, α34,3 . The initial conditions are:
In this section, we discuss the convergence of the HAM solutions. The convergence depends on the so-called convergence - control - parameter , and so, we plot the -curve for θg(0) , η (0) , θd(0), and rd (0). The interval of convergence is determined by the flat portion of the -curve. In order to plot the -curve we applied the HAM as given in  to our new model (21)-(24).
An introduction to Homotopy Analysis Method (HAM)
Consider the following differential equation:
where N is a nonlinear operator, is a vector of spatial variables, t denotes time and u is an unknown function.
Zero order deformation of HAM: By means of generalizing the traditional concept of homotopy, Liao  constructs the so-called zero-order deformation equation:
where is a non-zero auxiliary parameter, H is an auxiliary function, is an auxiliary linear operator, is an initial guess of ; is a unknown function. The degree of freedom is to choose the initial guess, the auxiliary linear operator, the auxiliary parameter, and the auxiliary function H . Expanding in Taylor series with respect to the embedding parameter p , one has
If the auxiliary linear operator, the initial guess, the auxiliary parameter, and the auxiliary function are so properly chosen that the above series converges at p =1, one has
which must be one of the solutions of the original nonlinear equation, as proved in . If the same initial guess and the same auxiliary linear operator are chosen, the approximations given by the homotopy perturbation method are exactly a special case of those given by the homotopy analysis method when = −1 and H =1 . The series (54) itself is in principle a kind of Taylor series (at p = 1). Hence, mathematically, homotopy perturbation method itself is also a kind of generalized Taylor technique.
mth-order deformation: Define the vector:
Differentiating Equation (51) m -times with respect to the embedding parameter p and then setting p = 0 and finally dividing the terms by m! , we obtain the m th-order deformation equation in the form of:
and m χm is the unit step function. Applying the inverse operator on both side of Equation (56), we get
In this way, it is easy to obtain um for m ≥1 , at m th-order and finally get the solution as:
In our model we choose the initial guess to be , , , and which satisfied the initial conditions. The linear operator will be:
with the property , where c1 and c2 are constants of integration. According to the system (21)-(24) and the terms as in Equation (32) the nonlinear operators will be defined as follows:
By substituting the series (32) into Equations (61)-(65) correspondingly we get the terms for Rm accordingly to Equation (57) as follows:
Now, the solution of the m th-order deformation can be expressed according to Equation (58) which can be solved by symbolic software such as Mathematica 8.0, Maple, Matlab and so on. We obtain a family of solutions that depends on the auxiliary parameter . So, regarding as independent variable, it is easy to plot the curves. For example, we can plot the curves:
The curves (i =1, 2,3) are a function of and thus can be plotted by a curve . According to  there exists a horizontal line segment (flat portion of the -curve) in the figure of and called the valid region of which corresponds to a region of convergent of the solutions. Thus, if we choose any value in the valid region of we are sure that the corresponding solutions series are convergent. For given initial approximation , the auxiliary linear operator , and the auxiliary function , the valid region of for different special quantities are often nearly the same for a given problem. Hence, the so-called -curve provides us a convenient way to show the influence of on the convergence region of the solutions series.
We compared the system dynamics of the models (21)-(24) and (37)-(48) with and without the impact of the thermal radiation. The results are based on the following diesel engines parameter values:
Diesel − engines
We studied the problem of the the effect of fuel spray polydispersity on the ignition process in a fuel cloud by applying numerical simulation, the homotopy perturbation method and the homotopy analysis method for 30th order deformation. We compared between the homotopy perturbation method and by solving the full system of the model i.e., the system of Equations: (21)-(24) for the solution profiles of the of the gas temperature, droplet temperature, radius and concentration numerically. Although we take into account only three different size of droplets, our results show that the homotopy perturbation method provides an excellent approximation of the solutions of the system with high accuracy (Figures 1-4).
Figure 1: Solution profiles of the gas temperature θg −τ ; 1: Full model solved numerically with the impact of the thermal radiation, 2: HPM model with the impact of the thermal radiation for the = −1, (2): HAM model with the impact of the thermal radiation for the = 0.04 , 3: Full model solved numerically without the impact of the thermal radiation, 4: HPM model without the impact of the thermal radiation, (4): HAM model without the impact of the thermal radiation for = 0.04 .
Figure 2: Solution profiles of the droplet temperature θd −τ ; 1: Full model solved numerically with the impact of the thermal radiation, 2: HPM model with the impact of the thermal radiation for the = −1, (2): HAM model with the impact of the thermal radiation for the = 0.04 , 3: Full model solved numerically without the impact of the thermal radiation, 4: HPM model without the impact of the thermal radiation, (4): HAM model without the impact of the thermal radiation for = 0.04 .
Figure 3: Solution profiles of the radius r −τ ; 1: Full model solved numerically with the impact of the thermal radiation, 2: HPM model with the impact of the thermal radiation for the = −1 , (2): HAM model with the impact of the thermal radiation for the = 0.04 , 3: Full model solved numerically without the impact of the thermal radiation, 4: HPM model without the impact of the thermal radiation, (4): HAM model without the impact of the thermal radiation for = 0.04 .
Figure 4: Solution profiles of the concentration η −τ ; 1: Full model solved numerically with the impact of the thermal radiation, 2: HPM model with the impact of the thermal radiation for the = −1, (2): HAM model with the impact of the thermal radiation for the = 0.04 , 3: Full model solved numerically without the impact of the thermal radiation, 4: HPM model without the impact of the thermal radiation, (4): HAM model without the impact of the thermal radiation for = 0.04 .
The gas temperature trajectory, Figure 1, for all models with and without the impact of the thermal radiation, starts with a fast increase in the temperature until then the temperature decreases from , which means cooling before ignition, until . This continuous process of cooling before ignition is summarized in Table 1 and corresponding in figures 1-4. This dimensionless time, τ , refers to the ignition time and it is compatible for all the solution profiles for the gas and droplets temperature, radius and concentration. According to these results, the F −M has the smallest ignition time with and without the impact of the thermal radiation when comparing to HPM and HAM. The HPM results are closer to the F −M than the HAM results with and without the impact of the thermal radiation.
Table 1: The process of cooling before the ignition time for the different models. F −M refer to the full-model, HPM with = −1 and HAM with = 0.04 .
The presence of the small parameter γ in the gas temperature equations, such that Equations (37)-(39) form a singularly perturbed system, enables one to exploit the geometrical version of the method of the integral manifold and hence to separate the model into fast and slow subsystems. According to this method, the above results are also sustained with MIM .
In order to clarify the influence of the thermal radiation on the dimensionless time before the final explosion of the system we introduce the term impact of the thermal radiation as given in  which is measured in percent and defined as:
here the subscript F −M refers to the full model solved numerically, τ is the ignition time, and the superscript rad and no − rad refer to the model with and without the impact of the thermal radiation respectively. We also defined the parameters , , and which show the difference between the different models in percent. The results are as follows: , , and , , and . As we can see from these results, the HPM model is closer to the model solved with numerical simulations than the HAM model (the comparison between the F −M and the HPM model based on the value of as equal to −1, and the comparison between the F −M and the HAM model based on the value of as equal to 0.04 ).
As we mentioned in the previous section, the convergence depends on the convergence-control-parameter , and so we plot the -curve for and as shown in Figure 5 for m = 30 in Equation (59) i.e. 30 th order approximation. According to Figure 5 , the interval of convergence that agrees for all of the corresponding solutions is . In order to emphasize the impact of the convergence-control parameter on the solutions profiles we defined the terms:
which point out the difference (in dimensionless-time) between the profiles solutions of the HPM and HAM methods for different with the impact of the thermal radiation and without the impact of the thermal radiation respectively from the numerical results, and τ refer to the ignition time. The results are summarized in Table 2. According to these results, the HPM and the HAM methods are closed to the numerical results for both with and without the impact of the thermal radiation.
Table 2: The impact of the convergence-control parameter on the solutions profiles with and without the impact of the thermal radiation comparing to the numerical simulations.
We have shown the the solutions obtained by HPM and HAM are convergent and that they extremely well with numerical simulations. It has also been shown that the homotopy perturbation method, which is a special case of the homotopy analysis method when = −1 , yields convergent solutions for all of the cases considered. These results demonstrate that HPM and HAM are very effective analytical methods for solving nonlinear problems in science and engineering.
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals