MATHEMATICAL MODELING OF ETHANE PYROLYSIS USING ENO SCHEMES

Abstract

In this paper an algorithm for constructing ENO schemes on unstructured grids for Navier-Stokes equations supplemented by the equations for the concentrations of chemicals is considered. A second order accuracy scheme, based on the linear reconstruction of a function in cells, is presented. The mathematical model was verified on the task of analyzing the gas-dynamic parameters of the gross reaction of ethane pyrolysis.

Keywords

• ENO scheme
• unstructured grid
• high order accuracy
• pyrolysis of ethane

Introduction

Mathematical modeling of gas-dynamic processes taking into account the chemical reactions is an urgent problem, the solution of which will allow one not only to describe the observed full-scale chemical experiments, but also to explore the behavior of the process under the conditions that cannot be realized experimentally [ 1 ]. Such processes are of complex nonlinear nature [ 2 ] and described by the system of NavierStokes equations, supplemented by the Arrhenius equations. To solve such systems one need to use numerical methods of high order accuracy [ 3 ] to get results, closest to the actual results of physical and chemical processes. In most cases studied processes occur in areas of complex shapes, which necessitate the use of unstructured grids. This paper describes a method of constructing the ENO scheme on a triangular grid. The idea of the algorithm is to restore the data in the cells by using a piecewise linear function with preservation of properties of TVD-schemes. This scheme allows one to calculate the flow of gas dynamics more accurately.

Mathematical model

The mathematical model of the processes under consideration is a set of conservation laws (conservation of mass, momentum and energy), supplemented by the equations for the concentrations of chemicals [ 4 ]:

   (v )  0,  t (v )

t (E) t    (v  v )  p  0,    (v (E  p))    (j h j J j )  Sh , (Yi )    (v Yi )    Ji  Ri .

t The equation of state of an ideal gas mixture is as follows: p  RT  (Yj / Mwj )

j moment:  T  x,t0   T0 ,  Y  x, t0   Y0.

 x, t0   0 , The initial conditions define the state of the object under consideration at the initial v  x,t0   (u  x,t0 , v  x,t0 , w x,t0 )T ,  p  x, t0   p0 , Sh  

j M wj The boundary conditions define the state of the gas on the border of the area under consideration. They can be the conditions of free inflow or outflow or the impermeability condition.

Here  is the gas density (kg / m3 ) , v is the velocity of the gas (m / s) , p is the pressure (Pa) , E is the total energy per unit volume (J / m3 ) , hj  TTref c pj dT is the enthalpy (J / kg) , cpj is the specific heat of the substance j (J / (kg  K )) , h j R j is the heating (cooling) as a result of a chemical reaction and external influences (J / (s  m3 )) , Yi is the local mass fraction of the substance i , Ri  M wi N Rˆ

r1 ir is the volume rate of the substance i (kg / (s  m3 )) , Rˆir is the molar Arrhenius coefficient of the formation (the flow rate) of the substance i in the course of the reaction r (1/ (s  m3 )) , M wi is the molecular mass of the substance i (kg), R is the universal gas constant (J / (mol  K )) , T is the temperature (K ) . In the energy equation and the equation for concentrations of chemicals the value of J j is the diffusion flux (kg / (s  m2 )) , which can be written as [ 4 ] T J j  DimYj  DiT . Here

DiT

is the thermal diffusivity coefficient

T (kg / (s  m)) , Dim is the coefficient of mass diffusion of the substance i into the substance j (m2 / s).

The simulation of gas flows on the example of ethane pyrolysis using ANSYS Fluent package was carried out in [ 5 ]. In this article the difference scheme on the basis of essentially non-oscillatory high-accuracy algorithm for the simulation of gas flows is proposed.

Difference scheme

We introduce a triangular grid   {Pi  (xi , yi ), i  1, 2,..., N} , containing all internal and boundary points of the computational domain. On the grid  we build a triangulation ()  {m  (Pim , Pjm , Pkm ), Pim , Pjm , Pkm  , m  1,..., M .

Then we construct a difference scheme that approximates the system of equations (1) - (2) using the integro-interpolation method: d dt

Ui (t)  1 i i F ndl  1 i i H ndl   Qdxdy 0

i

T Here Ui (t)  ,u,v,E,Yi  , F  (F1, F2 )T are the convective terms in the system (1), H  (H1, H2 )T are the diffusive terms in the system (1), Q is the source, n is the outer normal to the edge of the cell i .

 u2u p   uvv   00   00   00  F1   (Euvp)u  , F2   (Ev2pp)v  , H1    h0j J jx ; H2    h0j J jy , Q   S0h .

 uYi   vYi   j Jix   j Jiy   Ri  The integral in (3) is calculated by the Gauss quadrature formula:

3 q i F  n ds   i j1  j F (u  (G(jk) ,t), u  (G(jk) ,t))nk ,

k 1 (3) (4) where F(u (G(jk) ,t),u (G(jk) ,t))  nk are the discrete flows, which are calculated according to the scheme of rupture decay [ 6 ], u (G(jk) ,t), u (G(jk) ,t) are "left" and "right" vector values on the boundary of the cell, the algorithm for finding which is described in the next section.

Given a triangulation of the area {1, 2 ,...,  N } and given the average values of a function u(x, y) for each cell i (i  1, 2,..., N ) of the grid: (5) (6) (7) ENO algorithm ui  1 i  u(x, y)dxdy

i 1 i i p(x, y)dxdy ui following scheme [ 7 ]: u(0, j,k) (x, y)  u0  ux(j,k) (x  x0 )  u(yj,k) ( y  y0 ), where ux , uy are found from the solution to the system: x j x0 y j  y0 ux( j,k) xk  x0 yk  y0 u(yj,k)  u j  u0 uk  u0

. For each cell i (i  1, 2,..., N ) we construct a polynomial p(x, y) of degree at most k , which interpolates the function u(x, y) with the order of accuracy k 1 such that In order to construct the polynomials pi (x, y) of the k –th order for the cell i we (k 1)(k  2) 2 consider patterns Si  {1,..., K } , where K  is the number of cells in each of the patterns.

To construct a second order accuracy scheme for finding the gas-dynamic parameters in a cell 0 , we use the following pattern: S  {0 , i , ia , ib ,  j ,  ja ,  jb , k , ka , kb} (see the figure 1).

Choose the following patterns-candidates: S1  {0 ,  j , k }, S2  {0 , i , k }, S3  {0 , i ,  j },

S4  {0 , i , ia },

S5  {0 , i , ib},

S6  {0 ,  j ,  ja}, S7  {0 ,  j ,  jb}, S8  {0 , k , ka }, S9  {0 , k , kb}. Construct planes on these patterns. For example, for the pattern S1  {0 ,  j , k } we construct a plane by the Having multiple patterns-candidates, one can make a choice by the method described below. Let an operator F be defined: F[u(0, j,k) ,u(0,i,k) ,u(0,i, j) ,u(0,i,ia) ,u(0,i,ib) ,u(0, j, ja) ,u(0, j, jb) ,u(0,k,ka) ,u(0,k,kb) ]  u(0, j*,k*) . Here u(0, j*,k*) is the function, describing the plane having the smallest angle with the plane OXY. This is to ensure that the scheme had non-increasing property of the total variation (TVD-scheme) [ 8, 9 ].

Model of a chemical reaction of ethane pyrolysis

To check the adequacy of the constructed scheme the mathematical model based on the gross reaction of ethane pyrolysis was built. The table below shows the reaction scheme, the pre-exponential factor Ai and the activation energy Ei for each stage: Here A1  C2H6 , A2  C2H4 , A3  H2 , A4  CH4 , X i  Yi / M wi is the local molar fraction of the substance i , k j is the velocity constant of the j -th direct step (note that the dimension values are as follows: k1  [1 / s], k2  [l / (mol  s)] ). Dependence of the step rate constant on temperature is expressed from the Arrhenius equation [ 10 ]: Sh  (

RC2H6  The volumetric rate of formation of a substance i is determined from the relation:

N Ri  M wi r1ir wr . Here ir are the stoichiometric coefficients of the component i in the reaction r .

For such a mechanism of ethane pyrolysis heating or cooling as a result of chemical reactions is calculated as follows [ 4 ]:

hC02H6 hC02H4 hH02 hC0H4

M w,C2H6 M w,C2H4 M w,H2 M w,CH4 For each component of the mixture the heat capacity is approximated with respect to temperature by the following dependencies: cpC2H6  0, 0013T 2  4,8436T  483, 69,

RC2H4 

RH2 

RCH4 ). cpC2H4  0, 0015T 2  3, 9641T  505, 22, cpH2  0, 0016T 2  0, 366T 14274, 97, cpCH4  3, 3505T 1230, 9.

The thermal diffusivity coefficient and the coefficient of mass diffusion of the substance i into the substance j are calculated according to formulae [ 4 ]: DiT  2, 59 107 T 0,659  MMw0,wi05,i1511X1Xi i Yi   MMw0w,i05,i1418X9Xi i  , Dim 

1 X i  X i / Dij    3  1 T  , Dij  0, 00188   M wi

0,5  1 

 M wj 

, pabs 0,5(i  j ) D j, j i where pabs is the absolute pressure, D is the integral of diffusion collisions [ 4 ]: 0,1930 1,03587

1,76474 

  exp(0,47635TD* ) exp(1,52996TD* ) exp(3,89411TD* ) ,    ,  j are the Lennard-Jones parameters [ 4 ].    kb  j The technologies of hydrocarbon processing by controlling the chemical reactions of pyrolysis through the additional generation of radicals in the reacting gas using laser radiation [ 11 ], which leads to lower threshold temperature of the reaction and lower outlet temperature of the target products are being developed in the Boreskov Institute of Catalysis (the Siberian branch of the Russian Academy of Sciences). The scheme of experimental setup of the metal reactor with the externally heated reaction zone and the generation of radicals in the gas volume using laser radiation was previously presented in [ 4, 11, 12 ].

In contrast to the above-mentioned work, for the verification of numerical ENOmodels the reactor is considered as closed. This allows to demonstrate at the initial stage the performance of the numerical scheme for the calculation of gas dynamics taking into account the chemical transformations and heat and mass transfer. The calculations for the reaction zone are carried out. On the walls of the reaction zone the temperature T  900К. In the calculations we assume that the zone is closed and at the initial time it is filled with ethane of the temperature of T  300К. Figures shows the change in time of concentrations of the reaction products.

Conclusion

In this paper the numerical model for solving problems of gas dynamics, taking into account the chemical processes based on ENO schemes on unstructured grids is constructed. The essentially non-oscillatory scheme of higher order of accuracy is developed. The model is intended to describe the two-dimensional dynamics of the reacting gas taking into account heat exchange and diffusion processes as well as the chemical reactions of ethane pyrolysis. Computational experiments for comparison of the calculated and experimental data for ethane pyrolysis in a closed reactor with the heated reaction zone without laser radiation is carried out [ 11 ]. The calculated changes in concentrations of the substances involved in the chemical reaction of ethane pyrolysis, meet the basic physical and chemical laws. The degree of conversion of ethane at the time of the setting at these temperatures is about 40% (see the figure 2), which corresponds to the experimental data [ 4, 11 ]. For more accurate accounting of chemical transformations in the unsteady flow of the reacting gas it is planned to perform the modification of the model with the introduction of equations that are based on the radical mechanism of ethane pyrolysis [ 5 ]. It is assumed that such an accounting will allow one to calculate with higher accuracy not only the conversion of ethane with the change of the target product (ethylene) concentration, but also the relations between the reaction byproducts such as methane, hydrogen and others.

Acknowledgements

The work of Stadnichenko O.A. was supported by base budget project of Boreskov Institute of Catalysis SBRAS (Novosibirsk) V.44.2.11.

The work of Gubaidullin I.M. was supported by grant RFBR № 15-07-01764 A "Optimal control of chemical reactions metal complex catalysis."

References

1. Khusnutdinov RI , Baiguzina AR , Mukminov RR , Akhmetov IV , Gubaidullin IM , Spivak SI , Dzhemilev UM . New synthesis of pyrrole-2-carboxylic and pyrrole-2.5-dicarboxylic acid esters in the presence of iron-containing catalysts . Russian Journal of Organic Chemistry , 2010 ; 46 ( 7 ): 1053 - 1059 .

2. Balaev AV , Parfenova LV , Gubaidullin IM , Rusakov SV , Spivak SI , Khalilov LM , Dzhemilev UM . The mechanism of cp2zrcl2-catalyzed alkene cycloalumination with triethylaluminum to give alumacyclopentanes . Doklady Physical Chemistry , 2001 ; 381 ( 1 ): 279 - 282 .

3. Zhalnin RV . On the construction of a parallel computational algorithm for the direct numerical simulation of complex gas dynamic flows . Journal of Middle Volga Math. Soc. , 2008 ; 10 ( 1 ): 137 - 146 . [in Russian]

4. Stadnichenko OA , Snytnikov VN , Snytnikov VlN . Mathematical modeling of multicomponent gas flows with energy intensive chemical processes by the example of ethane pyrolysis . Numerical methods and Programming , 2014 ; 15 ( 4 ): 658 - 668 . [in Russian]

5. Nurislamova LF , Gubaidullin IM , Novichkova AV , Stoyanovskaya OP , Stadnichenko OA , Snytnikov VN . Few-step kinetic model of gaseous autocatalytic ethane pyrolysis and its evaluation by means of uncertainty and sensitivity analysis . Chemical Product and Process Modeling , 2014 ; 9 ( 2 ): 143 - 154 .

6. Godunov SK , Zabrodin AV , Ivanov MYa , Krayko AN , Prokopov GP . Numerical solution of multi-dimensional problems of gas dynamics . Moscow: Nauka, 1976 ; 400 p. [in Russian]

7. Ivanov IE , Kryukov IA . Quasimonotonous method of high order accuracy for the calculation of internal and jet flows of inviscid gas . Mathematical modeling , 1996 ; 8 ( 6 ): 47 - 55 . [in Russian]

8. Vyaznikov KV , Tishkin VF , Favorsky AP . Construction of monotone difference schemes of higher order of approximation for systems of hyperbolic equations . Mathematical modeling , 1987 ; 1 ( 5 ): 95 - 120 . [in Russian]

9. Harten A , Engquist B , Osher S , Chakravarthy S. Uniformly high order essentially nonoscillatory schemes . Journal of Computational Physics , 1987 ; 71 : 231 - 303 .

10. Stiller V . Arrhenius equation and nonequilibrium kinetics . Moscow: Mir, 2000 ; 176 p. [in Russian]

11. Snytnikov VN , Mishchenko TI , Snytnikov VlN , Malykhin SE , Avdeev VI , Parmon VN . Autocatalytic gas - phase dehydrogenation of ethane . Research on Chemical Intermediates , 2012 ; 38 ( 3 ): 1133 - 1147 .

12. Stadnichenko OA , Snytnikov VN , Snytnikov VlN , Masyuk NS. Mathematical modeling of ethane pyrolysis in a flow reactor with allowance for laser radiation effects . Chemical Engineering Research and Design , 2016 ; 109 : 405 - 413 . DOI: 10 .1016/j.cherd. 2016 . 02 .008.