Analysis combustion of fuel droplets release to the atmosphere using homotopy analysis method

Abstract: In this paper, we investigate the model of combustion of polydisperse fuel spray that release to the atmosphere. The mathematical model using the wellknown Eulerian–Lagrangian approach for transient flows of the fuel vapor–droplets mixture. The eddy breakup model is used to describe the combustion of the spray droplets as well as the k– model is used to describe the turbulence of the process. In order to investigate the mathematical/physical model, we applied the well-known analytic approximate method, the homotopy analysis method (HAM). According to the theory of HAM, the convergence and the rate of solution series are dependent on the artificial convergent control parameter ħ.


Introduction
Many chemical industries stored a large quantity of fuel at high pressure in liquid state. This phenomenon is dangerous for the environment and hence for people when a flammable substance released into the atmosphere. In many cases, two-phase outflow occurs so that the fuel cloud, which builds up in the atmosphere, contains a mixture of fuel vapor and liquid fuel droplets. A combustion/ ignition of a mixture of fuel droplet-vapor-air cloud can cause shock-free combustion with the formation of a fireball (Makhviladze, Roberts, & Yakush, 1999a, 1999b. Up to now, there are no many researches on the release of liquid fuel to the atmosphere. This paper is a natural continuous of the model presented in (Utyuzhnikov, 2002). In comparison to (Utyuzhnikov, 2002), we used with the operator ⟨⋅⟩ which is an operation of averaging over the droplet sizes, i.e. the droplets size are presented using with a probability density function (PDF) n(r). The governing equations are divided to gas-phase equations and liquid-phase equations. The gas-phase equations are described by a system of Favre averaged Navier-Stokes equations completed by the k-model of turbulence and eddy breakup model for turbulent combustion (Launder & Spalding, 1974). The liquid phase consists of ABOUT THE AUTHOR My research area is focused on all aspects of combustion of polydisperse and monodisperse fuel spray.
The models that describe the physical phenomena of these processes are systems of nonlinear partial differential equations. In order to investigate these models, we applied a semianalytical and asymptotic methods such as the homotopy analysis method (HAM), the method of integral invariant method (MIM), and singular perturbad homotopy analysis method (SPHAM).

PUBLIC INTEREST STATEMENT
In this paper, we applied the well-known semianalytical method the homotopy analysis method in order to investigate the problem of combustion of polydisperse fuel spray that release to the atmosphere.
interaction between the droplets. The continuity and energy equations include evaporation of droplets and inter-phase energy exchange using the source terms. The reaction rate for individual component are w i = ± i w, where i is the mass stoichiometric coefficient ( F = 1), the sign should be taken as negative for fuel and oxidizer and positive for combustion products. The radiative flux associated with kth gray gas is a function of the U k , i.e. the kth gas radiative energy density and U k , i.e. the black-body radiative energy density. The influence of turbulent fluctuation was neglected. The source terms in the gas phase determined by integrating over all droplets size.
Soot formation and oxidation were modeled using the two-step global scheme with large-scale flame. In order to apply the eddy breakup model, we use for calculation the mean flow characteristics. The radiation consists of infrared banded molecular radiation. The species having the strongest radiation bands being CO 2 and H 2 O and continuous spectra radiation from the soot particles. The optical thickness that corresponding to the k-th gray gas is determined by integrating the absorption coefficient along the vertical and horizontal lines of sight and taking the maximum of the resulting values overall all the computational domain. The model of volume emission is used for the gray gas. The radiative heat transfer calculations were based on the mean distribution of the temperature and concentration. In this way, one can neglect the influence of the turbulent fluctuations. The droplets are described by a Lagrangian approach. The volume fraction of the dispersed phase is negligible and there is no interaction between the droplets. This assumption can be made due to the spray dilute approximation. We take into account the effect of the turbulent fluctuations on the droplets motion in such a way that the gas velocity is represented as a sum of its average and fluctuating components. The evaporation rate of the droplets is described by the steady-state model in which the droplet temperature is determined by heat and mass balances with a Ranz-Marshal correction for the effect of the droplet motion. According to the assumptions above, the governing equations of the model are as follows (Makhviladze et al., 1999a(Makhviladze et al., , 1999bUtyuzhnikov, 2002  The following expressions are used in the above model: where ⃗ I is the unit tensor, = l + t , the total viscosity, is the sum of laminar and turbulent viscosities, respectively, C is a constant in the turbulence model, Q v is the heat of evaporation, is an absorption coefficient, a is a weight coefficient, ⃗ U ′ is the turbulent fluctuating velocity corresponding to the Gaussian probability distribution with the mean 2/3k.

The initial and boundary conditions
The initial and boundary conditions are as follows (Utyuzhnikov, 2002): The problem is considered as an axisymmetric one in the cylindrical coordinate system (r, z). The fuel and the droplets are injected with constant temperature and pressure. The injection velocity is represented by a piecewise function of time. The spatial distribution of the velocity follows the Gaussian distribution with a maximum value of U in (0) = √ (2(p 0 − p a )∕ l ) (using the incompressible Bernoulli equation, where p 0 is the pressure in the tank and l is the liquid propane density at the pressure p 0 ) on the z-axis and decrease in the velocity at the source boundary by the factor of √ 2.
The initial size distribution of droplet follows the gamma distribution. The heat capacity and the heat of evaporation are assumed to be constant. The droplets do not reach the right and the upper boundaries since we assumed that these boundaries are located far enough. On the solid wall the following boundary conditions for k and are k = 0, ∇ ⋅ ⃗ n = 0, where ⃗ n is the internal unit normal vector. For the radiative fluxes for optically thick gray gases, the following boundary conditions are Finally, we assumed that the ignition is spark like and takes place near the source surface in a small volume.

Preliminaries
In this section, we apply the HAM as presented in (Liao, 1998).

The basic idea of homotopy analysis method
Consider the following differential equation: where N is a nonlinear operator, ⃗ r is a vector of spatial variables, t denotes time, and u is an unknown function.

Zero-order deformation equation
By means of generalizing the traditional concept of homotopy, Liao in (1998) constructs the socalled zero-order deformation equation: where ℏ is a non-zero auxiliary parameter, H is an auxiliary function, is an auxiliary linear operator, u 0 (⋅) is an initial guess of u(⋅), Φ is an 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 Φ to a power series with respect to the embedding parameter p, one has where (2.18) 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. In this way, it is easy to obtain u m for m ≥ 1, at mth-order and finally get the solution as follows:

mth-order deformation
when m → ∞, we get an accurate approximation of the original Equation (3.1). If Equation (3.1) admits unique solution, then this method will produce the unique solution. If Equation (3.1) does not possess unique solution, the HAM will give a solution among many other (possible) solutions.
The constants are as follows: The parameters are as follows:

ℏ-Curve
According to the theory of HAM, the convergence and the rate of solution series are dependent on the convergent control parameter ℏ. This means that this parameter gives one a convenient way to adjust and to control the convergent region of the solutions. This subsection briefly describes how to choose this parameter based on (Liao, 2003).
(3.11) Solving the mth-order deformation one obtains 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 following curves: The curves Υ i (i = 1, 2, 3) are a function of ℏ and thus can be plotted by a curve Υ ≈ ℏ. According to (Liao, 2003), 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 convergence of the solutions. Thus, if we choose any value in the valid region of ℏ, we are sure that the corresponding solution series are convergent. For given initial approximations u 0 (⃗ r, t), v 0 (⃗ r, t) the auxiliary linear operator , and the auxiliary function H(⃗ r, t), the valid region of ℏ for different special quantities is often nearly the same for a given problem. Hence, the so-called ℏ-curve provides us with a convenient way to show the influence of ℏ on the convergence region of the solutions series.
As proved by (Shuijun, 2004) in general, if ℏ is properly chosen so that the series of solutions are convergent based on the following theorem: Theorem 1 Suppose that A ⊂ R be a Banach space donated with a suitable norm ‖ ⋅ ‖, over which the sequence u k (t) of u(r, p) = ∑ ∞ k=0 u k (r)p k is defined for a prescribed value of ℏ. Assume also that the initial approximation u 0 (t) remains inside the ball of the solution u(t). Taking r ∈ R be a constant, the following statements hold true: (1) If ‖u k+1 (t)‖ ≤ r‖u k (t)‖ for all k, given some 0 < r < 1, then the series solution u(t, p) = ∑ ∞ k=0 u k (t)p k converges absolutely at p = 1 to u(r) = ∑ ∞ k=0 u k (r) over the domain of definition of t, (2) If ‖u k+1 (⃗ r, t)‖ ≥ r‖u k (⃗ r, t)‖ for all k, given some r > 1, then the series solution u(t, p) = ∑ ∞ k=0 u k (t)p k diverges at p = 1 over the domain of definition of t.

Proof
(1) Let S n (t) denote the sequence of partial sum of the series u(r) = ∑ ∞ k=0 u k (r). We need to show that S n (t) is a Cauchy sequence in A, hence, consider, It should be remarked that owing to these inequalities, all the approximations produced by the homotopy method will lie within the ball of u(t). In order to show the Cauchy sequence test let n, m be two natural numbers such that n ≥ m. Using the above inequalities and the triangle inequality successively, we have, This means that S n (t) is a Cauchy sequence in the Banach space A, and this implies that the series solution u(r) = ∑ ∞ k=0 u k (r) is convergent. ✷ (2) Under the hypothesis supplied in (2), there exist a number d, d > r > 1, so that the interval of convergence of the power series u(r, p) = ∑ ∞ k=0 u k (r)p k is |p| < 1 d < 1 which excludes the case of p = 1. ✷ Theorem 2 If the series solution defined by u(r) = ∑ ∞ k=0 u k (r) is convergent, then it converges to an exact solution of the nonlinear problem N(u(r, t)) = 0. (3.13) Proof The proof presented in (Liao, 2003). ✷ Theorem 3 Assume that the series solution ∑ ∞ n=0 u n (t) is convergent to the solution u(t) for a prescribed value of ℏ. If the truncated series ∑ M n=0 u n (t) is used as an approximation to the solution u(t) of the problem N(u(r, t)) = 0, then an upper bound for the error, E M (t), is estimated as Proof Using the inequality in the proof of Theorem 1, we receive Since 0 < r < 1 then (1 − r n−M ) < 1 which implies the desired upper bound for the error, E M (t). ✷ Remark 1 An optimal value of the convergence control parameter ℏ can be found by means of the exact square residual error integrated in the whole region of interest Γ, at the order of approxima- Hence, the more quickly Res(ℏ) decreases to zero, the faster the corresponding homotopy series solution u(r) = ∑ ∞ k=0 u k (r) converges to the exact solution of the problem N(u(r, t)) = 0. So, at the given order of approximation M, the corresponding optimal value of the convergence control parameter ℏ is given by the minimum of Res(ℏ), corresponding to a nonlinear algebraic equation of the form dRes(ℏ)∕dt = 0 (see Table 4).

Analysis and results
In order to implement the HAM to the model (2.1)-(2.10), we write the expressions R m for each dynamical variables. The dynamical variables of the physical model are Hence, we defined the following series for the dynamical variables as follows: According to our model the linear operator [⋅] will be in the form of: and the nonlinear operator N[⋅] as follows: (4.1) According to (Liao, 2009) let: then the expressions for R m (⋅) are as follows: The results depend on the combination of the parameters k and according to the form Γ ≡ √ k 3 ∕ .
The effect of each of these parameters was weakly on the results. Another important parameter influencing on the results was the time which taken for the fireball to be visible. We denoted this time as t FireBall . The connections between Γ and t FireBall is when Γ increase the t FireBall decrease (numerical values of these parameters: Γ = 0.7m and t FireBall = 4s). The results of our calculations are presented in Figures 1 and 2. The injection period was about 0.2s. The fireball presented at the time instants t = 0.05, 0.15, 0.25, 0.95, 0.5, 1.5, 2.5, and 3.5. The dashed line in Figures 1 and 2 presented the visible fireball boundary in the experiment. Timely development of the fireball presented in these figures counterclockwise. The shape and location of the fireball are depends on the injection time. At the beginning of the process, the fireball is located close to the edge and it is small Figure 1. And as time progresses the location of the fireball approaches to the center and becomes more symmetric (Figure 1 counterclockwise). An interesting phenomenon occurs in Figure 2. The first three Figures 2(a)-(c) the fireball left at the center and is symmetric, but its length is shortened. In Figure  2(d), the fireball splits into two fireballs, this is because the high temperature and the breach diameter which effects on the shape of the fireball as long as the release time is relatively long. When applying the HAM we took into account different values of ℏ. The flat portion of the ℏ-curve called corresponds to a region of convergence of the solutions. Our simulations corresponds to 30th-order approximation ( Figure 3). As shown in Figures 1 and 2 our calculations agree very well with experimental data.
We compute the residual error according to the expression Res(ℏ) that corresponds to each nonlinear operator N[⋅] given in Equations (4.3)-(4.12). The results for the error are summarized in Tables 1-3. As shown in these tables when the m-th order deformation is increase the error decrease. (4.21)

Conclusions
A mathematical model was extended based on (Utyuzhnikov, 2002) by describing the droplets size using continuous model to simulate the combustion of polydisperse fuel spray vapor-droplets releases in the open atmosphere with hydrocarbon fuels at high pressures. The model takes into account two-phase clouds (fireball) combustion. The k-model is used to described the turbulent processes. The two phase model was described using Eulerian-Lagrangian approach for the transient flows of fuel vapor-droplets mixture. We also take into account soot formation, radiative heat transfer mass transfer, and momentum exchange between the gaseous and the dispersed phase. The main process of the ignition of the fireball is the evolution from the reaction initiation until the total fuel is burnout.
In order to investigate the proposed model, we applied the well-known analytic approximate method the HAM and we compared the results to experimental data. An optimal value approach for the convergence control parameter has also been given. Via the theorems provided here, not only the question of the convergence of the homotopy series is answered, but also the region of validity of the space variable ensuring the convergence is determined. The calculations of the dynamic of the fireball as well as the radiative fraction of the total combustion energy are shown a good agreement with the experimental fireball data.