Calculation of critical conditions for the filtration combustion model
                                                          O. Vidilina1, E. Shchepakina1
                                  1
                                      Samara National Research University, 34 Moskovskoe Shosse, 443086, Samara, Russia


Abstract

The paper is devoted to the study of the dynamic model of the autocatalytic combustion reaction in an inert medium with partial heat removal
from the reaction phase to the environment. We pay particular attention to modelling of the critical regime, which is a kind of a watershed
between the slow burning regimes and explosion modes. New algorithm for computing a critical value of the control parameter is presented.

Keywords: filtration combustion; thermal explosion; critical phenomena; singular perturbations; integral manifolds; canards


1. Introduction

   In last few years there was an increase in researches concerning multiphase combustions systems. The results of the studies
are widely used in the problems of safety of gas emissions, explosive dust clouds, mixture detonations, transportation and use of
combustible and explosive substances.
   In the present paper we consider a mathematical model of autocatalytic combustion reaction in a multiphase medium. The
multiphase nature of the process arises from the existence of inert phase alongside with the reactant phase. The inert medium
could correspond, for example, to a dusty medium or a porous matrix. We paid particular attention to the modelling of the
critical regime that is kind of a watershed between the slow burning regimes and thermal explosions.
   The main goal of the mathematical theory of thermal explosion [1-4] is to study the dynamics of the combustion process for a
given dimensions of the reactor, thermophysical and kinetic characteristics, heat transfer coefficient. These characteristics
correspond to parameters of the differential system, which is a mathematical model of the process. Under certain conditions of
values of these parameters, the reaction proceeds for as long as possible without transition into the explosion or the slow burning
mode. We call such regime a critical one.
   The goal of the present work is to determine the values of the parameters that correspond to the critical regime. In order to
find the values we consistently applied analytical and numerical methods. The main result of the paper is the derivation of the
algorithm used in calculation of a value of the control parameter of the system that corresponds to the critical regime. The
critical value of the control parameter is a solution of an algebra-differential system.

2. Model

   We consider combustion model of a rarefied gas mixture in an inert porous, or in a dusty, medium. We assume that the
temperature distribution and phase-to-phase heat exchange are uniform. The chemical conversion kinetics are represented by a
one-stage, irreversible reaction. The dimensionless model in this case has the form [5-7]:

               𝑑𝜃                           𝜃
           𝛾        = 𝜂 (1 − 𝜂)exp (              ) − 𝛼(𝜃 − 𝜃𝑐 ) − 𝛿𝜃,                                                               (1)
               𝑑𝜏                       1+𝛽𝜃

                𝑑𝜃𝑐
           𝛾𝑐         = 𝛼(𝜃 − 𝜃𝑐 ),                                                                                                  (2)
                𝑑𝜏

           𝑑𝜂                           𝜃
                = 𝜂 (1 − 𝜂)exp (             ),                                                                                      (3)
           𝑑𝜏                         1+𝛽𝜃


with initial condotions

           𝜂(0) = 𝜂0 /(1 + 𝜂0 ) = 𝜂̅0 , 𝜃(0) = 𝜃𝑐 (0) = 0.                                                                           (4)

   Here 𝜃 and 𝜃с are the dimensionless temperatures of the reactant phase and of the inert phase, respectively; 𝜂 is the depth of
conversion; 𝜏 denotes the dimensionless time; 𝜂0 is the parameter for autocatalyticity (this kinetic parameter characterizes the
degree of self-acceleration of the reaction: the lower the value, the more marked the autocatalytic reaction will be). The terms
−𝛿𝜃 and −𝛼(𝜃 − 𝜃𝑐 ) reflect the external heat dissipation and phase-to-phase heat exchange. The parameters 𝛾 and 𝛾𝑐
characterize the physical features of the reactant phase and of the inert phase, respectively. System (1)-(3) is singularly perturbed
since 𝛽 and 𝛾 are the small for typical combustible gas mixture [1-3].
   Depending on the relation between values of the parameters, the chemical reaction either moves to a slow regime with decay
of the reaction, or into a regime of self-acceleration which leads to an explosion. So, if we change the value of one parameter,
with fixed values of the other parameters, we can change the type of chemical reaction. Let us consider 𝛼 as a control parameter.
For some value of 𝛼 (we call it critical) the reaction is maintained and gives rise to a rather sharp transition from slow motions to
explosive ones. The transition region from slow regimes to explosive ones exists due to the continuous dependence of the system


3rd International conference “Information Technology and Nanotechnology 2017”                                                           151
                                                  Mathematical Modeling / O. Vidilina, E. Shchepakina
(1)-(3) on the parameter 𝛼. To find the critical value of the parameter 𝛼, it is possible to use special asymptotic formulae [4, 7,
8]. That approach was used in [5-7, 9] for system (1)-(3), in [10-19] for other laser and chemical systems, and in [20-24] for
some biological problems. In the next section the main results concerning this approach obtained for system (1)-(3) are given.
The realizability conditions for the critical regime were obtained in the form of a system of non-linear algebra-differential
equations, but the problem of calculating the critical value of the control parameter with the help of this system had not been
solved. The paper is devoted to develop an algorithm for calculating the critical parameter value. The readers can find details of
this algorithm in Sections 4 and 5.

3. Modelling of the critical regime

    We will consider the case 𝜂0 = 0 for the sake of simplicity, taking into account that for case 𝜂0 ≠ 0 the correction to the
initial conditions can be found with help of fast integral manifolds [4].
    The slow surface 𝑆 of system (1)-(3) is described by the equation (see Figure 1):
                              𝜃
             𝜂 (1 − 𝜂)exp (     ) − 𝛼(𝜃 − 𝜃𝑐 ) − 𝛿𝜃 = 0.
                               1+𝛽𝜃




                                                       Fig. 1. The slow surface of system (1)-(3).

This surface is a zero-order (𝛾 =0) approximation of a slow integral manifold of the system [4, 7, 8]. Recall, that the slow
integral manifold of a singularly perturbed system is defined as an invariant surface of slow motions, i.e., the flow on it has the
order O(1) as 𝛾 →0. Far from the slow surface, the fast variables of the system vary very rapidly, with a speed of order O(1/𝛾) as
γ→0.
   The intersection of the slow surface with the surface of irregular points (see Figure 2), given by the expression

                                 𝜃        1
              𝜂 (1 − 𝜂)exp (          )         −𝛼−𝛿 =0
                               1+𝛽𝜃 (1+𝛽𝜃)2


determines a breakdown curve. The breakdown curve separates the stable (𝑆 𝑠 ) and unstable (𝑆 𝑢 ) subsets of the slow surface S,
see Figure 3. System (1)-(3) has an stable integral manifold (𝑆𝜀𝑠 ) and an unstable integral manifold (𝑆𝜀𝑢 ) near 𝑆 𝑠 and 𝑆 𝑢 ,
respectively.
    When 𝛼 > 𝛼 ∗ the trajectories of the system starting at the initial point move along the stable branch 𝑆 𝑠 and the temperature
𝜃 does not reach relatively large values (see Figure 4). These trajectories correspond to the slow burning regimes.
    When 𝛼 < 𝛼 ∗ the system's trajectories, having reached the breakdown curve along 𝑆 𝑠 at the tempo of the slow variable, jump
into the explosive regime (see Figure 5).
    Due to the continuous dependence of the right-hand side of (1)-(3) on the parameter 𝛼 there are some intermediate trajectories
in the region between those shown above. For some value 𝛼 = 𝛼 ∗ we can glue the stable and unstable slow integral manifolds at
a point of the breakdown curve to get a canard [7, 8, 25, 26], i.e., the system’s trajectory which at first move along the stable
slow integral manifold and then continue for a while along the unstable slow integral manifold, see Figure 6.
    The canard describes the critical regime that separates the domain of slow burning modes and the domain of thermal
explosion. A deviation from the value of 𝛼 ∗ leads to the destruction of the gluing of stable and unstable slow integral manifolds
with a subsequent reaction’s transition either to the slow regime (when the trajectory of the system unfolds along a stable
manifold from the breakdown curve) or into the thermal explosion mode (when the trajectory, reaching the breakdown curve,
jumps from the slow manifold and rapidly runs away from it).




3rd International conference “Information Technology and Nanotechnology 2017”                                                 152
                                                      Mathematical Modeling / O. Vidilina, E. Shchepakina




                                                               Fig. 2. The surface of irregular points.




  Fig. 3. The intersection of the slow surface of system (1)-(3) with the surface of irregular points. On the intersection (the breakdown line) the stability of the
     slow manifold changes. The upper and lower sheets of the slow surface are stable, the part enclosed between the surface of irregular points is unstable.




                                  Fig. 4. The trajectory (left) and the 𝜃-component (right) in the case of a slow birning regime:
                                                   𝛼 = 3, 𝛽 = 0.1, 𝛾 = 0.001, 𝛾𝑐 = 0.7, 𝜂̅0 = 0.02, 𝛿 = 0.02.




                                       Fig. 5. trajectory (left) and the 𝜃-component (right) in the case of thermal explosion:
                                                    𝛼 = 0.7, 𝛽 = 0.1, 𝛾 = 0.001, 𝛾𝑐 = 0.7, 𝜂̅0 = 0.02, 𝛿 = 0.02.




3rd International conference “Information Technology and Nanotechnology 2017”                                                                                     153
                                                           Mathematical Modeling / O. Vidilina, E. Shchepakina




        Fig. 6. The canard and the 𝜃-component (right) in the case of critical regime: 𝛼 = 0.949, 𝛽 = 0.1, 𝛾 = 0.001, 𝛾𝑐 = 0.7, 𝜂̅0 = 0.02, 𝛿 = 0.02.

   To calculate the critical value of the parameter 𝛼 = 𝛼 ∗ and the asymptotic expressions for the corresponding canard

              𝛼 ∗ = 𝛼0 + 𝛾𝛼1 + 𝑜(𝛾),

              𝜃(𝜂, 𝛾) = 𝜑0 (𝜂) + 𝛾𝜑1 (𝜂) + 𝑜(𝛾),                                                                                                   (5)

              𝜃𝑐 (𝜂, 𝛾) = 𝜓0 (𝜂) + 𝛾𝜓1 (𝜂) + 𝑜(𝛾),

we use the usual method of eliminating an independent variable. In this case, the system (1)-(3) takes the form

                  𝑑𝜃                       𝜃                                 𝜃
              𝛾        𝜂 (1 − 𝜂)exp (              ) = 𝜂 (1 − 𝜂)exp (                ) − 𝛼(𝜃 − 𝜃𝑐 ) − 𝛿𝜃,
                  𝑑𝜂                    1+𝛽𝜃                               1+𝛽𝜃

                   𝑑𝜃𝑐                         𝜃
              𝛾𝑐         𝜂 (1 − 𝜂)exp (             ) = 𝛼(𝜃 − 𝜃𝑐 ).
                   𝑑𝜂                      1+𝛽𝜃


   We substitute (5) into these equations to get

                                                          𝜑0                 𝜑                    𝜑0                𝜑
              𝛾(𝜑0′ + 𝛾𝜑1′ )𝜂 (1 − 𝜂)exp (                                   1
                                                                ) [1 + 𝛾 (1+𝛽𝜑 )2
                                                                                  ] = exp (                          1
                                                                                                        ) [1 + 𝛾 (1+𝛽𝜑 )2
                                                                                                                          ]
                                                      1+𝛽𝜑0                          0          1+𝛽𝜑0                   0


                                          −(𝛼0 + 𝛾𝛼1 )(𝜑0 − 𝜓0 + 𝛾(𝜑1 − 𝜓1 )) − 𝛿(𝜑0 + 𝛾𝜓1 ) + 𝑜(𝛾),                                               (6)

                                                           𝜑0                    𝜑
              𝛾𝑐 (𝜑0′ + 𝛾𝜑1′ )𝜂 (1 − 𝜂)exp (                                   1
                                                                  ) [1 + 𝛾 (1+𝛽𝜑 )2
                                                                                    ]
                                                          1+𝛽𝜑0                      0


                                                      = (𝛼0 + 𝛾𝛼1 )(𝜑0 − 𝜓0 + 𝛾(𝜑1 − 𝜓1 )) + 𝑜(𝛾).                                                 (7)

   Setting 𝛾 = 0 in (6) and (7), we obtain

                                   𝜑0
              𝜂 (1 − 𝜂)exp (              ) − 𝛼0 (𝜑0 − 𝜓0 ) − 𝛿𝜑0 = 0,                                                                             (8)
                                 1+𝛽𝜑0

                                           𝜑0
              𝛾𝑐 𝜑0′ 𝜂 (1 − 𝜂)exp (                 ) = 𝛼0 (𝜑0 − 𝜓0 ).                                                                             (9)
                                          1+𝛽𝜑0


From the equation of the breakdown curve, taking into account (5), we have

                                     𝜑0∗              1
              𝜂 ∗ (1 − 𝜂 ∗ )exp (              )                − (𝛼0 + 𝛿) = 0,                                                                  (10)
                                    1+𝛽𝜑0∗ (1+𝛽𝜑0∗ )2


where (𝜂 ∗ , 𝜑 ∗ , 𝜓 ∗ ) is the gluing point of the slow integral manifolds.
  After double differentiation (8) with respect to 𝜂, with taking into account (10), we get one more condition at the gluing point:

                                    𝜑0∗
              (1 − 2𝜂 ∗ ) exp (            ) + 𝛼0 𝜓0′∗ = 0, 𝜓0′∗ = 𝜓0′ (𝜂 ∗ ).                                                                   (11)
                                  1+𝛽𝜑0∗


   Thus, the expressions (8)-(11) give us the zeroth order approximations of the critical value of the control parameter and the
canard.
   Further, in order to find the first order approximations, we equate the coefficients of 𝛾 in the first degree in the system (6),
(7). As a result, we get:
                                𝜑                         𝜑       𝜑1
             𝜑0 𝜂 (1 − 𝜂)exp ( 0 ) = [𝜂 (1 − 𝜂)exp ( 0 )               2 −
                                                                           (𝛼0 + 𝛿)] 𝜑1 + 𝛼0 𝜓1 − 𝛼1 (𝜑0 − 𝜓0 ),         (12)
                                    1+𝛽𝜑0                                  1+𝛽𝜑0 (1+𝛽𝜑0 )


3rd International conference “Information Technology and Nanotechnology 2017”                                                                            154
                                                                    Mathematical Modeling / O. Vidilina, E. Shchepakina
                                           𝜑0      𝜑 (𝛾 𝜓 −1)                      ′
              𝜂 (1 − 𝜂)exp (      ) [𝜑0′ + 𝛾𝑐 𝜓1′ + 1 𝑐 0 2 ] = −𝛿𝜑1 .                                                     (13)
                             1+𝛽𝜑               0   (1+𝛽𝜑 )                            0
                            1                                                          𝜑0∗
              𝛼1 =                [𝛼0 𝜓1′∗ − 𝜑0′∗ 𝜂 ∗ (1 − 𝜂 ∗ )exp (                        )].                           (14)
                       𝜑0∗ −𝜓0∗                                                   1+𝛽𝜑0∗


Here 𝜑0′∗ = 𝜑0′ (𝜂 ∗ ).
   The expressions (12)-(14) determine the first order approximations of the critical value of the control parameter and the
canard. It should be noted that to calculate the values of 𝛼0 and 𝛼1 from (8)-(14) it is necessary to apply numerical methods. The
development of the algorithm for finding the critical value of the control parameter is our next goal.

4. The gluing point

   In order to verify the correctness of the algorithm, developed in the present paper, we can use some specific case when an
analytical solution of (8)-(14) is available and compare it to the one yielded by our method. For this goal we now consider the
case 𝛿 =0 which corresponds the absence of external heat dissipation.
   In the case 𝛿 =0 system (1)-(3) possesses a first integral

              𝜂 −𝛾𝜃 − 𝛾𝑐 𝜃𝑐 = 0.

With the help of this first integral, we can reduce the order of system (1)-(3) by eliminating the variable 𝜃𝑐 . As a result we obtain
a plane system:

                  𝑑𝜃                                     𝜃
              𝛾        = 𝜂 (1 − 𝜂)exp (                         ) − 𝛼(1 + 𝛾/𝛾𝑐 )𝜃 + 𝛼/𝛾𝑐 (𝜂 − 𝜂̅0 ),
                𝑑𝜏                                   1+𝛽𝜃
              𝑑𝜂                                     𝜃
                   = 𝜂 (1 − 𝜂)exp (                        ).
              𝑑𝜏                                    1+𝛽𝜃


Here we have to deal with the slow curve rather than then slow surface [5-7, 9]. The coordinates of the gluing point of the
integral manifolds for some value 𝛼 = 𝛼0 can be found from the self-intersection conditions of the slow curve, which in our
case have the form

                                             𝜃∗                         𝛼
              𝜂 ∗ (1 − 𝜂 ∗ ) exp (                   ) − 𝛼𝜃 ∗ +             𝜂 ∗ = 0,                                       (15)
                                           1+𝛽𝜃 ∗                      𝛾𝑐
                                           𝜃∗              𝛼
              (1 − 2𝜂 ∗ )exp (                      )+          = 0,                                                       (16)
                                       1+𝛽𝜃 ∗             𝛾𝑐
                                          𝜃∗                    1
              𝜂 ∗ (1 − 𝜂 ∗ ) exp (                   )                 − 𝛼 = 0.                                            (17)
                                           1+𝛽𝜃 ∗ (1+𝛽𝜃 ∗ )2


   From (15) and (17) we get

              𝜂 ∗ = 𝛾𝑐 𝜃 ∗ − (1 + 𝛽𝜃 ∗ )2 𝛾𝑐 .                                                                             (18)

   Using (15) and (16), we obtain

              (1 − 2𝜂 ∗ )(𝜃 ∗ − 𝛾с−1 𝜂 ∗ ) + 𝛾с−1 𝜂 ∗ (1 − 𝜂 ∗ ) = 0.

   From here it follows that

              𝛾𝑐 (1 + 𝛽𝜃 ∗ )4 = 𝛾𝑐 𝜃 ∗ 2 − 𝜃 ∗ .                                                                           (19)

   Equations (18) and (19) give us the values 𝜃 ∗ and 𝜂 ∗ . Further, from equation (16) we obtain
                                      𝜃∗
           𝛼0 = 𝛾𝑐 (2𝜂 ∗ − 1) exp (      ∗ ).            1+𝛽𝜃


   Using the smallness of the parameter 𝛽, from (19) we can analytically find the value 𝜃 ∗ in the form of an asymptotic
expansion with respect to the parameter 𝛽: 𝜃 ∗ = 𝜃00
                                                  ∗      ∗
                                                     + 𝛽𝜃01 + o(𝛽). Similar expansions can also be written for 𝛼0 and 𝜂 ∗ . In
the case 𝛽 = 0 we have:

               ∗        1
              𝜃00 = (𝛾с−1 + √4 + 𝛾с−2 ),
                        2

               ∗         ∗
              𝜂00 = 𝛾𝑐 (𝜃00 − 1),
                                 ∗ )
                            exp(𝜃00
              𝛼00 =                    .
                        2+√4+𝛾с−2




3rd International conference “Information Technology and Nanotechnology 2017”                                                     155
                                                               Mathematical Modeling / O. Vidilina, E. Shchepakina
5. Algorithm for calculating the critical value of the control parameter

   Let us describe an algorithm for solving the system of nonlinear algebra-differential equations (8)–(11) with initial conditions
that we obtained from (4)–(5):
                           𝜂0
               𝜂(0) =            = 𝜂̅0 , 𝜑0 (0) = 𝜑0 .
                         1+𝜂0


   We consider numerical solution of the system. Using (8) we get:

                           𝛿           𝜂(1−𝜂)                𝜑0
       𝜓0 = 𝜑0 (1 +            )−                 exp (             ).
                          𝛼0                𝛼0              1+𝛽𝜑0


   Taking the derivative yields:

               1                                                     𝜑0                               𝜑0        𝜑′
       𝜓0′ =        [𝜑0′ (𝛼0 + 𝛿) + (1 − 2η)exp (                           ) − 𝜂(1 − 𝜂) exp (                    0
                                                                                                            ) (1+𝛽𝜑 )2
                                                                                                                       ].
               𝛼0                                                   1+𝛽𝜑0                           1+𝛽𝜑0             0


   Now, let us substitute 𝜓0 and 𝜓0′ into (9):
                                                                                                                      𝜑0
                𝛼                𝛼0 𝛿𝜑0                     (1−2𝜂)            𝜑0                       𝜂(1−𝜂) exp(       )
                                                                                                                   1+𝛽𝜑0
       𝜑0′ = ( 0 −                             𝜑0       +            exp (           )) / (𝛼0 + 𝛿 −                  2       ).
                   𝛾𝑐   𝛾𝑐 𝜂(1−𝜂) exp(             )        𝜂(1−𝜂)           1+𝛽𝜑0                           (1+𝛽𝜑0 )
                                             1+𝛽𝜑0


   We get a first order ordinary differential equation with respect to 𝜑0 = 𝜑0 (𝜂), where 𝜂 is a variable and 𝛼0 is a constant.
   For any given value of 𝛼0 we can solve the equation using Runge-Kutta fourth-order method. This method has a good
precision and, in spite of its laboriousness, is widely used to obtain a numerical solution for ordinary differential equations.
Moreover, we can further benefit from this method by using an adaptive stepsize.
   Next, in order to the coordinates of the gluing point we find 𝜂 ∗ and 𝜑0∗ from (10) and (11):

                                 𝜑0 ∗               1
       𝜂 ∗ (1 − 𝜂∗ ) exp (                   )               − (𝛼0 + 𝛿) = 0 ,
                                1+𝛽𝜑0 ∗ (1+𝛽𝜑0∗ )2

                                                                          𝜑0 ∗
                                𝜑0 ∗                 𝜂 ∗ (1−𝜂 ∗ ) exp(          )−𝛿𝜑0 ∗
                   ∗)                                                   1+𝛽𝜑0 ∗
       (1 − 2𝜂 exp (                     ) + 𝛼0                                𝜑   ∗      = 0.
                            1+𝛽𝜑0 ∗                     𝛾𝑐 𝜂 ∗ (1−𝜂 ∗ ) exp(     0   )
                                                                            1+𝛽𝜑0 ∗


    We can also solve this nonlinear system numerically using any of existing iterative methods. Let us note that sometimes it is
possible to solve such systems by applying the elimination and back substitution method. However, in vast majority of cases
iterative methods are used. Next, we substitute the solution 𝜂 ∗ into (8) to obtain 𝜑0∗∗ . Minimizing the difference between 𝜑0∗ and
𝜑0∗∗ , we can find 𝛼0 with arbitrary precision.
    Analogously, we get 𝛼1 from (12)-(14) as we already computed 𝛼0 .
    The algorithm was implemented using Java. To verify the correctness of the algorithm we considered a special case (𝛽 = 0,
𝛾 = 0, 𝛿 = 0), that allows us to solve the system analytically. Then we compared the analytical solution with the numerical
solution yielded by the algorithm. The results for 𝛼0 are presented in Table 1.

                                   Table 1. The comparison between analytical and numerical solutions for 𝛼0 for various values of 𝛾𝑐 .
                                                                                  𝛼0
                                       𝛾𝑐        Analytical solution                             Numerical solution
                                       0.6       1.837200                                        1.837291
                                       0.7       1.566012                                        1.566038
                                       0.8       1.393921                                        1.393923
                                       0.9       1.275978                                        1.275919
                                       1.1       1.125980                                        1.125923
                                       1.2       1.075605                                        1.075623
                                       1.3       1.035258                                        1.035253


6. Conclusion

   We have studied the mathematical model of filtration combustion of combustible gas in an autocatalytic reaction case. We
have shown that the critical phenomena of the model can be described by the canard. The critical regime is a kind of watershed
between the safe processes and thermal explosion regimes. We have established that it is possible to control the mode and,



3rd International conference “Information Technology and Nanotechnology 2017”                                                             156
                                                    Mathematical Modeling / O. Vidilina, E. Shchepakina
therefore, the combustion process by adjusting the parameter characterizing the heat removal from the reaction phase to the
external environment.
   We have developed a new algorithm that allows to compute the critical value of the control parameter by combining
analytical methods of the geometric theory of singular perturbations and numerical methods. The presented algorithm can be
used in other similar problems for studying critical phenomena in dynamic systems and calculating critical values of control
parameters.

Acknowledgements

   The contribution of O. Vidilina was supported by the Russian Foundation for Basic Research and Samara region (grant 16-
41-630529-p) and the Ministry of Education and Science of the Russian Federation as part of a program of increasing the
competitiveness of SSAU in the period 2013–2020. E. Shchepakina was supported by the Ministry of Education and Science of
the Russian Federation (Project RFMEFI58716X0033).

References

[1] Spalding von DB. Combustion and Mass Transfer. Oxford – New York: Pergamon Press, 1979; 409 p.
[2] Zeldovich YaB, Barenblatt GI, Librovich VB, Makhviladze GM. The Mathematical Theory of Combustion and Explosions. New York: Consultants Bureau,
     1985; 597 p.
[3] Babushok VI, Goldshtein VM, Sobolev VA. Critical Condition for the Thermal Explosion with Reactant Consumption. Combust. Sci. and Tech. 1990; 70:
     81–89.
[4] Sobolev VA, Shchepakina EA. Model Reduction and Critical Phenomena in Macrokinetics. Moscow: “Energoatomizdat” Publisher, 2010; 320 p. (in
     Russian)
[5] Gol'dshtein V, Zinoviev A, Sobolev V, Shchepakina E. Criterion for thermal explosion with reactant consumption in a dusty gas. Proc. of London R. Soc.
     Ser. A 1996; 452: 2103–2119.
[6] Gorelov GN, Shchepakina EA, Sobolev VA. Canards and critical behavior in autocatalytic combustion models. Journal of Engineering Mathematics 2006;
     56(2): 143–160.
[7] Shchepakina E, Sobolev V. Black swans and canards in laser and combustion models. Singular Perturbation and Hysteresis. Philadelphia: SIAM, 2005; 207–
     255.
[8] Shchepakina E, Sobolev V, Mortell MP. Canards and Black Swans. Singular Perturbations. Introduction to System Order Reduction Methods with
     Applications. Lecture Notes in Mathematics 2014; 2114: 141–182.
[9] Golodova ES, Shchepakina EA. Modeling of safe combustion at the maximum temperature. Mathematical Models and Computer Simulations 2009; 1(2):
     322–334.
[10] Gorelov GN, Sobolev VA. Mathematical modeling of critical phenomena in thermal explosion theory. Combust. Flame 1991; 87: 203–210.
[11] Gorelov GN, Sobolev VA. Duck-trajectories in a thermal explosion problem . Appl. Math. Lett. 1992; 5(6): 3–6.
[12] Sobolev VA, Shchepakina EA. Self-ignition of dusty media . Combus Explosion Shock Waves 1993; 29: 378–381.
[13] Sobolev VA, Shchepakina EA. Duck trajectories in a problem of combustion theory. Differential Equations 1996; 32: 1177–1186.
[14] Shchepakina E. Black swans and canards in self-ignition problem. Nonlinear Analysis: Real Word Applications 2003; 4: 45–50.
[15] Shchepakina E. Canards and black swans in model of a 3-D autocatalator. J. Phys.: Conf. Series 2005; 22: 194–207.
[16] Shchepakina E, Korotkova O. Condition for canard explosion in a semiconductor optical amplifier. Journal of the Optical Society of America B: Optical
    Physics 2011; 28(8): 1988–1993.
[17] Shchepakina E, Korotkova O. Canard explosion in chemical and optical systems. Discrete and Continuous Dynamical Systems – Series B 2013; 18(2):
     495–512.
[18] Shchepakina EA. Critical phenomena in a model of fuel's heating in a porous medium. CEUR Workshop Proceedings 2015; 1490: 179–189.
[19] Firstova N, Shchepakina E. Conditions for the critical phenomena in a dynamic model of an electrocatalytic reaction. J. Phys.: Conf. Series 2017; 811:
     012002.
[20] Sobolev V. Canard cascades. Discrete and Continuous Dynamical Systems – Series B 2013; 18(2): 513–521.
[21] Gavin C, Pokrovskii A, Prentice M, Sobolev V. Dynamics of a Lotka-Volterra type model with applications to marine phage population dynamics. J. Phys.:
     Conf. Series 2006; 55(1): 80–93.
[22] Pokrovskii A, Shchepakina E, Sobolev V. Canard doublet in a Lotka-Volterra type model. J. Phys.: Conf. Series 2008; 138: 012019.
[23] Pokrovskii A, Rachinskii D, Sobolev V, Zhezherun A. Topological degree in analysis of canard-type trajectories in 3-D systems. Applicable Analysis 2011;
     90(7): 1123–1139.
[24] Sobolev VA. Canards and the effect of apparent disappearance. CEUR Workshop Proceedings 2015; 1490: 190–197.
[25] Benoit E, Callot JL, Diener F, Diener M. Chasse au canard. Collect. Math. 1981; 31-32: 37–119. (in French)
[26] Shchepakina E, Sobolev V. Integral manifolds, canards and black swans. Nonlinear Analysis. A 2001; 44: 897–908.




3rd International conference “Information Technology and Nanotechnology 2017”                                                                           157