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. 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.
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]:
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 (
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 with initial condotions = ( − ), = (1 − )exp(
), 1+ = (1 − )exp( 1+
) − ( − ) − ,
(0) = 0/(1 + 0) = ̅0, (0) = (0) = 0.
(
Mathematical Modeling / O. Vidilina, E. Shchepakina
(
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 (
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(
The intersection of the slow surface with the surface of irregular points (see Figure 2), given by the expression (1 − )exp( 1+ )
(1+1 )2 − − = 0
determines a breakdown curve. The breakdown curve separates the stable ( ) and unstable ( )subsets of the slow surface S,
see Figure 3. System (
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 (
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).
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.
To calculate the critical value of the parameter = ∗ and the asymptotic expressions for the corresponding canard
∗ = 0 + 1 + ( ),
( , )= 0( )+ 1( )+ ( ),
( , )= 0( )+ 1( )+ ( ),
we use the usual method of eliminating an independent variable. In this case, the system (
(1 − )exp(1+ ) = (1 − )exp(1+ )− ( − )− , (1 − )exp(1+ ) = ( − ).
( 0′ + 1′) (1 − )exp(1+ 0 0)[1 + (1+ 10)2] = exp(1+ 0 0)[1 + (1+ 1 0)2]
−( 0 + 1)( 0 − 0 + ( 1 − 1))− ( 0 + 1)+ ( ),
( 0′ + 1′) (1 − )exp(1+ 0 0)[1 + (1+ 1 0)2]
Setting = 0 in (
From the equation of the breakdown curve, taking into account (
1 1+ 0∗ (1+ 0∗)2 − ( 0 + )= 0, where ( ∗, ∗, ∗) is the gluing point of the slow integral manifolds.
(1 − 2 ∗)exp(1+ 0∗ 0∗)+ 0 0′∗ = 0, 0′∗ = 0′( ∗).
= ( 0 + 1)( 0 − 0 + ( 1 − 1))+ ( ).
After double differentiation (
Thus, the expressions (
The expressions (
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 (
In the case =0 system (
− = 0.
). (1 − )exp( 0 ) [ 0′ + 1′ +
(11(+ 0′0−)21)] = − 1.
case have the form
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
(
∗ = 0, ) + = 0, 1 1+ ∗) (1+ ∗)2 − = 0.
From (
∗ = ∗ − (1 + ∗)2 .
From here it follows that (1 − 2 ∗)( ∗ − с−1 ∗)+ с−1 ∗(1 − ∗) = 0. (1 + ∗)4 = ∗2 − ∗. 0 = (2 ∗ − 1)exp( ∗
1+ ∗). 0∗0 = 1 2
( с−1 + √4 + с−2), 0∗0 = ( 0∗0 − 1), 00 = exp( 0∗0) .
2+√4+ с−2
Equations (
Using the smallness of the parameter , from (
that we obtained from (
0 0 (1− )exp(1+ 0 0) + (1−2 ) (1− ) exp( 0 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 (
1 1+ 0∗) (1+ 0∗)2 − ( 0 + ) = 0 , (1 − 2 ∗)exp( 0∗ 1+ 0∗ ) + 0 ∗(1− ∗) exp(1+ 0∗0 ∗)− 0∗
∗ ∗(1− ∗) exp(1+ 0 0∗)
= 0. 0∗∗, we can find 0 with arbitrary precision.
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 (
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.
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.
The contribution of O. Vidilina was supported by the Russian Foundation for Basic Research and Samara region (grant 1641-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).