Mathematical Modeling Critical phenomena in a model of fuel's heating in a porous medium Shchepakina E.A. Samara State Airspace University Abstract. The autoignition of flammable liquid in an inert porous medium are studied. In this paper we concentrated on the critical case which is concerned with the phenomenon of delayed loss of stability in the dynamical model. The realizability conditions for the critical regime are obtained. It is shown that critical regime is modelled by a canard − a trajectory of slow–fast system, which first move near the stable part of the slow invariant manifold, then move near the unstable part of it. Keywords: ignition, critical phenomenon, canard, invariant manifold, delayed loss of stability. Citation: Shchepakina E.A. Critical phenomena in a model of fuel's heating in a porous medium. Proceedings of Information Technology and Nanotechnology (ITNT-2015), CEUR Workshop Proceedings, 2015; 1490: 179-189. DOI: 10.18287/1613-0073-2015-1490-179-189 Introduction This paper deals with the investigation of the critical conditions for autoignition of combustible fluids in porous insulation materials [1, 2]. This phenomenon is usually caused by a leaking of a combustible liquid into lagging material surrounding a hot pipework. Due to highly insulation environment heat losses are remarkable low and autoignition may occur as a result of exothermic oxidation reaction. The investigation of the autoignition process has been carried out by many authors, see for instance [3–12] and references therein. We shall study a process which may be defined as the autoignition in two-phase medium (combustible liquid and inert porous matrix). The possible depletion of oxygen or its diffusion into porous structure and transport of the liquid or its vapour within the insulation are all ignored, in order to focus attention on the competitive effects of the reactive term of the dispersed liquid and evaporative heat loss. The dimensionless model in this case has the form [1] du  QK1 xe1/ u  (u  ua )  Qc K 2 xe βe / u , dt (1) dx   K 2 xeβe / u  K1 xe 1/ u . dt 179 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… Here, u is a dimensionless temperature of the reactant phase; the dimensionless concentration x represents the mass fraction of combustible liquid present in the porous material; the dimensionless parameters Q and K1 characterize the heat of reaction and the reaction frequency, respectively, for the exothermic oxidation reaction, while Qc and K 2 are the similar terms for the endothermic evaporation reaction; β e is the ratio of the enthalpy of vaporization to the activation energy of the oxidation reaction; ua is ambient temperature. Introducing the new variables θ and τ by u  θ  β2 , t  τ exp(1/ β), β  u(0), and taking into account  1   1   1 exp     exp   exp    ,  β(θβ  1)   1+βθ   β leads (1) to the form QK1 x  θ   ua  β  1 θ exp     θ  2  exp   β2  1+βθ   β  β Qc K 2 x  1-β   βθ   exp  e  exp  e  , (2) β2  β   1+βθ   1-β   βθ   θ  x   K 2 x exp  e  exp  e   K1 x exp  .  β   1+βθ   1+βθ  We introduce the new parameters  1  1-β  u β ε= exp    , a  K 2 exp  e  , θ a  a 2 ,  β  β  β QK  1 QK  -β  μ  2 1 exp    , ν= c 2 2 exp  e  . β  β β  β  Due to the smallness of the parameter ε for typical combustible liquids, system (2) can be rewritten in the singularly perturbed form (see, for instance, [13−17]):  θ   βe θ  εθ  μx exp     θ  θ a   νx exp  , (3)  1+βθ   1+βθ   βθ   θ  x  ax exp  e   K1 x exp  . (4)  1+βθ   1+βθ  The chemically relevant phase space  of system (3), (4) is defined by   x  0, θ  1/  . In [2] system (1) was investigated numerically under quasi-steady-state assumption that corresponds to the assumption ε  0 for system (3), (4). This approach allows determining the main types of chemical regimes of the investigated process. The quasi- steady-state assumption is widely used in the theory of combustion and gives good results 180 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… to draw conclusions about the qualitative behavior of the full system for sufficiently small ε . However, the critical phenomena are highly sensitive with respect to the parameters. Hence, this fact implies the considerable difficulties under the numerical calculations. Thus, detailed study of this mathematical object is possible with taking into account the small perturbations and using of asymptotic methods, for example, methods of the integral manifolds theory for singularly perturbed systems. Criteria for the critical regimes The trivial solution is the final steady state of the system. The degenerate equation  θ   βe θ  0  μx exp     θ  θ a   νx exp    F ( x,θ) describes the slow curve S of  1+βθ   1+βθ  (3), (4) (see, for example, [17, 18]). The subset S s ( S u ) of S with F ( x,θ) 0 ( 0) θ is called the stable or attractive (unstable or repulsive) part of S. A point A on S in which F / θ=0 is called the jump or turning point. Stable and unstable parts of the slow curve are zeroth order approximations of corresponding stable and unstable slow invariant manifolds. The invariant manifolds lie in an ε-neighborhood of the slow curve, except near jump or turning points (see [17] and references therein). The slow curve has an asymptote θ  θ , where 1 θ  ln  μ / ν  , βe  1  β ln  μ / ν  and intersects the axis Oθ in the point with θ  θ0 . The shape of the curve S varies with the relation between values of the parameters, which leads to a change in qualitative behavior of the system. So, if we change the value of one parameter, with fixed values of the other parameters, we can change the type of chemical reaction. Following [2], we consider β e as a control parameter. For θ  θ0 we have βe  b0 where 1  βθ b0  1  ln  μ / ν  . θa 181 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… a) b) c) d) Fig. 1. – The slow curve of (3), (4) for b0  1 , β>ua , and a) βe  bT  b0 ; b) bT  βe  b0 ; c) βe  b0 ; d) βe  b0 182 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… a) b) Fig. 2. – a) The trajectory (the solid line) and the slow curve (the dashed line) of (3), (4), and b) the x− and θ −components of the solution in the case of the slow regime: b=0.3685 Consider the case β>ua and b0  1 . For βe  bT the lower branch of S can consist of two stable parts ( S1s ) and one unstable part ( S1u ), see Fig. 1a. These parts are divided by two turning points, which merge with one another and disappear at a value βe  bT [2, 11], see Fig. 1b. In both cases the trajectories of system (3), (4) move along the stable part of the slow curve to the final steady state. These trajectories correspond to the slow regimes, which are safe, see Fig. 2. For βe  b0 the upper and lower branches of S consist of stable ( S1s and S 2s ) and unstable ( S1u and S 2u ) parts, which are divided by the turning points A1 and A2 , see Fig. 1d. And the system's trajectories starting at any point of the basin of attraction of S 2s correspond to the slow regimes, see Fig. 3. In other case, when the initial point is out of the basin of attraction of S 2s , we can observe the thermal explosion (Fig. 4) or thermal explosion with delay [8, 18]. The thermal explosion with delay occurs when the initial point belongs to the basin of attraction of S1s and the system's trajectories having reached the jump point A1 along S1s at the tempo of the slow variable jump into the explosive regime. 183 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… a) b) Fig. 3. – The case of the slow regime: b=1; the initial point belongs to the basin of attraction of S 2s For βe  b0 the point A1 merges with A2 to give one self–intersection point A of the slow curve, see Fig. 1c. As it was noted above, in ε-neighborhood of the subset s u S1s ( S 2u ) there exists a stable (unstable) slow invariant manifold S1,ε ( S2,ε ). For some value βe  b*  b0  O(ε), (ε  0) [11], the stable and unstable slow invariant s manifolds S1,ε u and S2,ε are glued at the point A . As a result for βe  b* system (3), (4) has a canard trajectory [17-22] which, at first, follows an attractive invariant manifold, and then a repulsive one. In both cases the distances travelled are O(1) as ε  0 , see Fig. 5. This canard simulates the critical regime, separating slow chemical regimes from regimes with a self–acceleration in the case b0  1 . If β>ua , one can observe the similar transformation of the slow curve (and the qualitative behavior of the system) as shown in Fig. 1 but with a decreasing value of the parameter β e ( bT  b0 in this case). For b0  1 plots of the slow curve are mirror images with respect to the vertical axis of the graphs shown in Fig. 1, and for nonsignificant values of the initial concentration of a combustible liquid the thermal behavior of the chemical system is safe. Otherwise, the value βe  bT determines the boundary of the safe region [2]. 184 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… a) b) Fig. 4. – The case of thermal explosion: b=1; the initial point is out of the basin of attraction of S 2s Our goal is to reveal the sufficient conditions for realization of the critical regime for the case b0  1 . As it has been noted above, the main feature here consists in fact that during the critical regime the temperature attains a high value but without explosion. The interest in critical phenomena is occasioned by not only for reasons of safety, but in many cases the critical regime is the most effective in technological processes [6-9, 18, 20-22]. Realizability conditions for the critical regime Using the method of integral manifolds and the canard techniques in [17,18] it is possible to find the critical value of the parameter βe  b* and corresponding trajectory in the form of the asymptotic representation θ  φ( x,ε)  φ0 ( x)  εφ1 ( x)  o(ε), (5) βe  b*  b0  εb1  o(ε). (6) We write (3), (4) as 185 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model…   βθ   θ  εθ  ax exp  e   K1 x exp     1+βθ   1+βθ    βθ   θ   θ  θ a +νx exp  e   μx exp  ,  1+βθ   1+βθ  or, taking into account (5), (6),   φ0   φ1  x  εφ0  ε 2 φ1   K1 exp   1  ε   1+βφ0   2  1+βφ 0    b φ  1+βφ0  b1φ0  b0 φ1   a exp  0 0  1  ε  1+βφ0   1+βφ0   2   b φ  1+βφ0  b1φ0  b0 φ1   φ0  εφ1  θ a +νx exp  0 0  1  ε  1+βφ0   1+βφ0   2  (7)  φ0   φ1  μx exp   1  ε   o(ε). 1+βφ0   2  1+βφ0   a) b) Fig. 5. – a) The trajectory (the solid line) and the slow curve (the dashed line) of (3), (4), and (b) the x− and θ −components of the solution in the case of critical regime: b  b*  0.56827 186 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… Setting ε=0 in (7) we obtain the slow curve equation:  φ0  F ( x, φ0 )  μx exp    φ0  1+βφ0  (8)  bφ   θ a  νx exp  0 0   0.  1+βφ0  Conditions for self-intersection of the slow curve in point A  xs ,φ0 ( xs )  F F  0 (9) x  xs ,φ0 ( xs ) φ0  x ,φ ( x ) s 0 s give us the coordinates of the self-intersection point and the zeroth-order approximations for critical value b* . Indeed, from (8), (9) we obtain 1  βθa  θa xs  , φ0 ( xs )  θ a (10) μ exp θ a / 1  βθ a  ln  ν / μ  and 1  βθ a b0  1  ln  μ / ν  . (11) θa Equating the coefficients with ε1 in (7) we get   bφ   φ0   xφ0  a exp  0 0   K1 exp     1+βφ 0   1+βφ 0    x   bφ   φ0     φ1 1+ b ν exp  0 0   μ exp  2  0   (12)  1+βφ0    1+βφ 0   1+βφ 0     b φ  bφ  νx exp  0 0  1 0 .  1+βφ 0  1+βφ 0  From (9) we note that the expression in brackets in r.h.s. of (12) is equal to zero at point A . To avoid a discontinuity in function φ1 ( x) at xs we put   b φ (x )   φ0 ( xs )   φ0 ( xs )  a exp  0 0 s   K1 exp     1+βφ0 ( xs )   1+βφ 0 ( xs )    b φ (x )  b φ (x )  ν exp  0 0 s  1 0 s .  1+βφ0 ( xs )  1+βφ 0 ( xs )  From this and (10) we get φ0 ( xs ) 1  βθ a    (1- b0 )θ a   b1 =  K1 exp    a . (13) νθ a   1  βθ a   187 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… Note, that the value φ0 ( xs ) can be found from equation (8) after double differentiation with respect to x, taking into account (9):  μ 1  βθ a  ν 3  θ  φ0 ( xs )  2  ln exp  a    θa μ  1+βθ a   (14)  1  βθ a  μ    2β 1  βθ a   2  ln  .  θa ν Expressions (13) and (14) give us 2μ 3 1  βθ a  4 ν  3θ a  K1 a  b1 = ln 3 exp     θa 4 μ  1+βθ a   μ ν  (15)  1  βθ a  ν    2β 1  βθ a   2  ln  .  θa μ Thus, the expressions (8), (11), (12), and (15) determine the first−order approximation for canard and corresponding critical value βe  b* . It should be noted that it is not possible to explicitly solve equation (8) with respect to φ 0 , while the critical value b* has been found in the explicit form. However, one can use the implicit or parametric representation for slow invariant manifold [17] to obtain an approximation of the canard. Conclusion In this paper the model of autoignition of combustible fluids in an inert porous medium has been studied. The realizability conditions for the critical regime have been obtained as the explicit asymptotic expression for the control parameter. It was shown that the critical regime is modelled by the canard. This regime plays the role of a watershed between the safe processes and regimes with self-acceleration that leads to the explosion. It should be noted that the critical regime is not a slow regime, since the temperature may attain a high value, and is not explosive, as the temperature increases at the tempo of the slow variable. Thus, for the examined model the new type of the safe regime has been revealed. Acknowledgements This work is supported in part by the Russian Foundation for Basic Research (grants 13-01-97002-p, 14-01-97018-p) and the Ministry of education and science of the Russian Federation in the framework of the implementation of Program of increasing the competitiveness of SSAU for 2013−2020 years. 188 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Shchepakina E.A. Critical phenomena in a model… References 1. McIntosh AC, Bains M, Crocombe W, Griffiths JF. Autoignition of combustible fluids in porous insulation materials. Combust Flame, 1994; 99: 541-550. 2. McIntosh AC, Griffiths JF. On the thermal runaway of combustible fluids in lagging material. IMA Journal of Applied Mathematics, 1996; 55: 83-96. 3. Spalding Von DB. Combustion and mass transfer. Oxford, New York: Pergamon Press, 1979. 4. Zeldovich YaB, Barenblatt GI, Librovich VB, Makhviladze GM. The mathematical theory of combustion and explosions. New York: Consultants Bureau, 1985. 5. Babushok VI, Goldshtein VM, Sobolev VA. Critical condition for the thermal explosion with reactant consumption. Combustion Science and Technology, 1990; 70: 81-89. 6. Gorelov GN, Sobolev VA. Mathematical modeling of critical phenomena in thermal explosion theory. Combust Flame, 1991; 87: 203-210. 7. Gorelov GN, Sobolev VA. Duck−trajectories in a thermal explosion problem. Applied Mathematics Letters, 1992; 5(6): 3-6. 8. Gol'dshtein V, Zinoviev A, Sobolev V, Shchepakina E. Criterion for thermal explosion with reactant consumption in a dusty gas. Proceedings of the Royal Society of London A, 1996; 452: 2103-2119. 9. Shchepakina E. Black swans and canards in self−ignition problem. Nonlinear Analysis: Real Word Applications, 2003; 4: 45-50. 10. Kuo K−K. Principles of Combustion. New York: Wiley, 2005. 11. Sobolev VA, Shchepakina EA. Model reduction and critical phenomena in macrokinetics. Moscow: “Fizmatlit” Publisher, 2010. [in Russian] 12. Sazhin S. Droplets and sprays. London, New York, Heidelberg: Springer, 2014. 13. Mishchenko EF, Rozov NKh. Differential equations with small parameters and relaxation oscillations. New York: Plenum Press, 1980. 14. O'Malley RE. Singular perturbation methods for ordinary differential equations. Applied Mathematical Sciences, 1991; 89. 15. Mishchenko EF, Kolesov YuS, Kolesov AYu, Rozov NKh. Asymptotic methods in singularly perturbed systems. New York: Plenum Press, 1995. 16. Vasilieva AB, Butuzov VF, Kalachev LV. The boundary function method for singular perturbation problems. In: Series in Applied Mathematics. Philadelphia: SIAM, 1995; 14 p. 17. Shchepakina E, Sobolev V, Mortell MP. Singular perturbations. Introduction to system order reduction methods with applications. In: Lecture Notes in Mathematics. Cham − Heidelber − NewYork − Dordrecht − London: Springer, 2014; 2114 p. 18. Shchepakina E, Sobolev V. Black swans and canards in laser and combustion models. In: Mortell, M.P., O'Malley, R.E., Pokrovskii, A., Sobolev, V. (eds) Singular perturbation and hysteresis. Philadelphia: SIAM, 2005; 207-255. 19. Benoit E, Callot JL, Diener F, Diener M. Chasse au canard. Collectanea Mathematica, 1981-1982; 31-32(1-3): 37-119. [in French] 20. Sobolev VA, Shchepakina EA. Duck trajectories in a problem of combustion theory. Differential Equations, 1996; 32: 1177-1186. 21. Shchepakina EA, Sobolev VA. Integral manifolds, canards and black swans. Nonlinear Analysis: Theory, Methods & Applications - Series A, 2001; 44(7): 897-908. 22. Gorelov GN, Shchepakina EA, Sobolev VA. Canards and critical behavior in autocatalytic combustion models. Journal of Engineering Mathematics, 2006; 56: 143-160. 189 Information Technology and Nanotechnology (ITNT-2015)