- Email: [email protected]
An equation-oriented approach to modeling heterogeneous catalytic reactors
Accepted Manuscript An equation-oriented approach to modeling heterogeneous catalytic reactors N.V. Vernikovskaya PII: DOI: Reference:
S1385-8947(17)30844-6 http://dx.doi.org/10.1016/j.cej.2017.05.095 CEJ 16993
To appear in:
Chemical Engineering Journal
Please cite this article as: N.V. Vernikovskaya, An equation-oriented approach to modeling heterogeneous catalytic reactors, Chemical Engineering Journal (2017), doi: http://dx.doi.org/10.1016/j.cej.2017.05.095
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
AN EQUATION-ORIENTED APPROACH TO MODELING HETEROGENEOUS CATALYTIC REACTORS
N. V. Vernikovskaya
Boreskov Institute of Catalysis SB RAS, Novosibirsk 630090, Pr. Ak. Lavrentieva 5, Russia Novosibirsk State University, Novosibirsk, Russia Novosibirsk State Technical University, Novosibirsk, Russia E-mail [email protected]
Abstract At multiscale modeling of heterogeneous catalytic reactors there appear systems comprising a large number of equations. The solution of such systems is an arduous task and available algorithms require voluminous computations. A modeling strategy has been suggested for such systems. An effective solution algorithm has been developed for a large class of models based on the application of an identical set of numerical tools such as integro-interpolation method, method of straight lines, a special case of a second-order Rosenbrock method, tridiagonal matrix algorithm or Thomas algorithm on each scale of a multiscale reactor model. Step size control is implemented with account for the rate of change of the variables on each scale. The efficiency and robustness of this algorithm were demonstrated at multiscale modeling of heterogeneous catalytic reactors such as tubular reactors, monolith catalytic reactor, and fluidized bed reactor. When switching between reactor types at multiscale modeling it is not necessary to modify the algorithm. The algorithm can also be used for multiscale modeling of heterogeneous catalytic reactors of other types.
Keywords: Catalytic reactors, multiscale modeling, numerical methods.
1. Introduction Mathematical modeling of diverse chemical processes running in various reactors can help to reduce experimental efforts during chemical process development. In chemical reactor modeling chemical reactors are represented as space-time hierarchical structures. The space scale ranges from 10-11 to 10 3 m, the associated time scale ranges from 10 -15 to 108 s [1,2]. The hierarchical multiscale approach is a backbone of the mathematical modeling of chemical reactors [1,3]. This approach consists in decomposing a complex chemical-technological process into chemical and physical components, studying these components independently and subsequently carrying out the synthesis of a general mathematical model from the models of separate parts of a complex process. Atomic-molecular processes and heat- and mass transfer in porous catalyst pellets or in the catalytic layer represent the lower two scales of the mathematical models of heterogeneous catalytic reactors. The processes of heat- and masstransfer in the reactor, i.e. in the catalyst bed or monolith channel, represent the third scale, and modeling of a switch device with regard to the processes of mixing, heat exchange, etc – the forth scale. A multiscale mathematical model can describe the preceding levels and given level simultaneously. The history of development and application of the hierarchical approach to the modeling of catalytic processes is described in detail in [1]. Recently there have been many studies on the multiscale modeling of catalytic processes in reactors of various types. Generally, to treat a developed mathematical model there are two groups of mathematical modeling tools based on two different approaches: block-oriented (BO) and equation-oriented (EO) [4]. The BO approaches address modeling on the flowsheet or reactor levels [5,6]. In the first case, every process is abstracted by a flowsheet consisting of the coupled units, e.g. heatexchanger, compressor, reactor, absorber, etc. However for the reactor most simulators
2
provide yield based reactor models or a set of idealized models like continuous stirred tank reactor (CSTR), plug flow, etc. with limited capabilities for defining reactor kinetics [6]. In the second case, every reactor model is decomposed into simpler objects such as reaction kinetics, thermodynamics, reactor geometry, reactor configuration, etc. Data are exchanged between the individual objects. BO approaches have many advantages. However reactor configurations are limited to a fixed set of reactors with given models: adiabatic bed reactor, tubular reactor, fluidized bed reactor, and every unit or object is solved one by one in a predefined sequence. This approach is basically used for modeling processes on the same scale. Besides, BO approaches do not support the implementation of more detailed models of standard reactors or models of non-standard reactors: monolith reactor, radial flow reactor, moving bed, multiphase reactors, membrane processes etc. EO approaches have some advantages over BO approaches, for they support the implementation of more complex reactor models and the simulators of EO process solve the system of equations simultaneously. The last-mentioned advantage indicates that EO approaches are more flexible as compared to BO approaches and, therefore, are more suitable for multiscale modeling of heterogeneous catalytic reactors. Among the techniques within equation-oriented approaches to reactor modeling one can distinguish computational fluid dynamics (CFD). CFD (commercial or in-house software) provides an excessive data set with a limited physical meaning, so CFD flow models can be validated only with a great deal of experimental data. This technique is too time-consuming. Another technique within EO approaches is using the models of necessary and sufficient degree of sophistication. One cannot always improve a model by adding parameters. These models satisfy Occam’s razor principle, which means the models “contain the elements that are needed, but without extra embellishments” [7]. Therefore, compared to CFD-based models EO approach-based models use less parameters and can be measured or estimated
3
independently. Validation of such models requires less experimental data. Transfer coefficients can be obtained from the correlations which capture also the effect of reactor scale. The models and simulations that accurately treat both the physics and chemistry of real systems and the parameters independently obtained from measured or estimated data sets are currently considered to be the most appropriate approach to modeling various chemical reactors [8,9]. In the case of multiscale modeling of heterogeneous catalytic reactors the EO approach-based technique allows constructing models of dimensionality ≤ 2D with the firstorder derivative with respect to a coordinate on each scale for a large class of reactors. Simultaneous solution of a system of equations as a whole simultaneously on all scales under consideration is an advantage of EO approaches, but the solution cost increases significantly as the size of the model increases [10]. In this study a novel strategy is proposed for multiscale EO process simulation. This strategy is based on the use of a single set of numerical methods on each scale. This strategy is applicable to models of dimensionality ≤ 2D with the first-order derivative with respect to a coordinate on each scale under consideration.
2. Field of applications of the proposed strategy In the case of multiscale modeling of heterogeneous catalytic reactors the use of the models with the elements that are needed, but without extra embellishments for a large class of reactors permits the model to be presented on each of the scales under consideration as a system of differential equations of dimensionality ≤ 2D with the first-order derivative with respect to a coordinate. This means that the system of partial differential equations can contain both first- and second-order derivatives with respect to one of the coordinates and only first-order derivatives with respect to the other coordinate (Fig. 1 a) or time (Fig. 1b).
4
All types of heterogeneous catalytic gas-solid reactors can be subdivided into three groups: packed bed, fluidized bed, and barrier wall reactors [11].
2.1. Processes on the third (reactor) scale The class of packed bed reactors includes two specific subgroups: tubular (1) and adiabatic (2) reactors. For steady-state processes in the case (1) it is necessary to take account of radial convection and diffusion heat- and mass transfer processes, whereas along the reactor length, which is normally several times larger than the radius, it is sufficient to take account of only the convection transfer, i.e. first-order derivative with respect to length. In the case (2) for steady-state processes it is not necessary to take account of radial heat- and masstransfer processes, while it is possible to consider axial heat- and mass transfer with any degree of sophistication, i.e. by taking account of both convection and diffusion components. In adiabatic reactors radial heat- and mass transfer can be neglected because there is usually no heat supply or removal through the reactor wall and, as a result, no gradient of temperature or substances concentrations across the reactor radius. As for the effects of reactor wall on layer porosity, the layer porosity can fluctuate within the distance of several particle diameters from the reactor wall for regularly packed beds and only one diameter for randomly packed beds, which is typical of all reactor types [12,13]. In contrast to tubular reactors, in the adiabatic ones the reactor to particle diameter ratios are normally very large so the effect of porosity fluctuations on the interstitial gas velocity and, consequently, the radial heat- and mass transfer can be neglected. To obtain the first-order derivative, in the model it is possible to use the method of transition to the non-steady state problem. Furthermore, in this case it is also possible to solve the non-steady state problem. This can be used in the description of the process in the monolith reactor.
5
(a)
(b)
Figure 1. Schematic representation of the structure of partial differential equations: (a) – with respect to space coordinate “y” contains only first-order derivatives, (b) – the first-order derivative with respect to time.
For steady-state processes in fluidized bed reactors it is not necessary to take account of the radial heat- and mass transfer, whereas it is possible to consider axial heat- and mass transfer with any degree of sophistication, i.e. by taking account of both convection and diffusion components. In fluidized bed reactor radial heat- and mass transfer can be neglected because of intensive catalyst mixing. To obtain the first-order derivative in the model it is possible to use the method of transition to the non-steady state problem and also to solve the non-steady state problem. Barrier wall reactors can be characterised by a barrier that incorporates a catalyst placed between two gas spaces or chambers with different gas composition [11]. Within this barrier the chemical reaction takes place and product flows may be manipulated to go to one or both sides using selective permeability or pressure gradient. In the reactors of this type for steadystate processes it is necessary to take account of the processes of radial convective and diffusive heat- and mass transfer, including the barrier, whereas along the reactor’s length it sufficient to consider the convective transfer, i.e. the first-order derivative with respect to length.
6
2.2. Description of processes on the second (catalyst) scale In the foregoing descriptions of reactor scale processes the mathematical models for catalyst pellets or catalytic layer, if necessary, can be simplified as one-dimensional. In this case to determine the first-order derivative in the model it is possible to use the method of transition to the non-steady state for a steady-state reactor model or to solve the non-steady state problem for a non-steady state reactor model.
2.3. Description of processes on the first (atomic-molecular) scale On this scale one considers either a detailed molecular nonsteady-state kinetic model for surface reactions and gets the first-order derivative with respect to time or a detailed molecular steady-state kinetic model for surface reactions, and in the latter case it is possible to use the method of transition to the non-steady state problem. A multiscale reactor model may incorporate the models of all three scales or the models of only two scales. In any case the equations of a lower-scale model are to be solved in each higher scale grid point. The combination of numerical methods suggested to solve the equations is identical for the models on every scale. This combination can be simplified only for the model of the first-scale processes, in which heat- and mass transfer processes are not considered. The suggested combination of numerical methods can be used in modeling the process for a single-scale model.
3. Solution method Without the loss of generality, put down the equations for a third-scale model in cylindrical coordinates. For clearness and simplification of the algorithm description, assume that radial and axial transfer velocities are constant and can be taken outside the derivative:
7
∂ui ∂y
+ ai ⋅
∂u 1 ∂(x ⋅ ui ) 1 ∂ bi ⋅ x ⋅ i = ci , i = 1, M ⋅ − ⋅ x ∂x x ∂x ∂x
(1)
Initial conditions:
y = 0 : ui = ui0
(2)
Boundary conditions:
∂ui
x =0: x =R:
∂x ∂ui ∂x
bM ⋅
= 0; i = 1, M ; = 0; i = 1, M − 1; ∂uM
∂x
(3)
= α ⋅ (Tw − uM )
where u = (u1,..., uM ) is the vector of the concentrations of substances and temperature,
ai , bi , ci are the coefficients depending nonlinearly on the solution, i.e. on the concentrations of substances and temperature, α is the heat-transfer coefficient from the tube wall to the bed, J/(m2.s.K), Tw is tube wall temperature, K.
3.1. Integro-interpolation method Discrete analogues of the model equations were constructed by using the integrointerpolation method. Transition to discrete models cannot be made in a purely formal way. Such a transition must leave in force as many basic properties of the source object models as possible [14]. In respect to the equations of system (1) with initial (2) and boundary (3) conditions such a transition means that for the corresponding difference scheme it is necessary that the discrete analog of the mass and energy conservation laws be fulfilled. A common approach to the construction of such difference scheme is based on writing this law in an integral form for cells x j −1 ≤ x < x j , x j = j ⋅ h j , j = 1, N of the difference grid 8
followed by substitution of the obtained integrals and derivatives with approximate difference equations, step h j can be obtained from the condition of equal cylinder section areas. A complete description of integro-interpolation method can be found online in Inline Supplementary Algorithm Description S1.
3.2. Method of straight lines We do not introduce discrete formulas for derivatives with respect to “y”, thus transforming partial differential equation (PDE) for each ui into a set of ordinary differential equations (ODEs) for nodal values along the reactor radius. The process of transformation of equation (1) into system of ODEs for each variable ui is an example of application of the method of straight lines or the method of semi-discretization [15]. By successively applying this procedure to all inner nodes and taking account of the boundary conditions (3) in the relation at j=1 and j=N, one gets a system of equations for each variable ui = (u1,..., uM ) written in the form: dui = A ⋅ u + C , ∀i = 1, M dy i i i
(4)
where three diagonal matrix Ai and vector Сi depend on the solution in a non-linear manner. The method of straight lines is advantageous because to solve the semidiscrete form of initial differential equation in partial derivatives it is possible to use different methods for ordinary differential equations. Let us put down system (4) for one variable ui , for simplicity, skipping index “i”: du = f (u ), f (u ) = A ⋅ u + C , u(0) = u 0 dy
(5)
9
where A is a tridiagonal matrix of dimensionality NxN, С is a vector of dimensionality N which are nonlinearly dependable on the solution, u, f(u) are vectors of dimensionality N.
3.3. A special case of a second-order Rosenbrock method When solving the problems with initial data described by ordinary differential equations the best results are normally obtained with either linear multistep or Runge-Kutta methods. In 1963, Rosenbrock [16] proposed the methods that differ from the explicit Runge-Kutta schemes by the regularization of the right part of a differential problem. Such methods can be derived from a class of semi-implicit methods if each kni is computed with the use of not more than one Newton’s iteration. Later on the methods of the specified type were called Rosenbrock-type methods [17]. They can be specified as follows: m
un +1 = un + ∑ pi ⋅ kni , I − ai ⋅ h ⋅ fu
i =1
i −1 ⋅ un + ∑ γij ⋅ knj ⋅ kni = h ⋅ f j =1
i −1 ⋅ un + ∑ βij ⋅ knj , j =1
(6)
where I is a unity matrix, fu is Jacobi matrix of the right hand size of system (5), ai , pi , γij , βij are the parameters. For Rosenbrock-type methods to be implemented most effectively all ai must be mutually equal, and γij = 0. In this case at each step it is necessary to compute and invert only one matrix of dimensionality N. To solve the system of equations (5) let us make use of a second-order L-stable formula as follows [17]: un +1 = un + (0.5 + 0.25 2) ⋅ kn1 + (0.5 − 0.25 2) ⋅ kn 2
(7)
Dn = I − (1 − 0.5 2) ⋅ h ⋅ fu (un )
(8)
Dn ⋅ kn 1 = h ⋅ f (un )
(9)
Dn ⋅ kn 2 = h ⋅ f (un + 2 ⋅ kn1 )
(10)
10
In this case the Jacobi matrix is of dimensionality NxN and systems of linear equations (9) and (10) can be solved using some laborious algorithms such as LU decomposition with the number of operations proportional to N3. As applicable to our system of equations (5), in which matrix A is tridiagonal, it may be preferable to use the tridiagonal matrix algorithm concurrently with the iteration process. A complete description of this procedure can be found online in Inline Supplementary Algorithm Description S2.
3.4. Step size control Step h is chosen so that for the global error the following relation is fulfilled [17-19]: (11)
kn 2 − kn1 / un < 15 ⋅ ε
The stepsize adjustment algorithm is based on the computation of the values of the right part and difference between kn1 and kn2 for each equation of system (7): R = un + (0.5 + 0.25 2) ⋅ kn1 + (0.5 − 0.25 2) ⋅ kn 2
(12)
RE = kn 1 − kn 2
(13)
RQ = kn 1 − kn 2 / un
(14)
Then maximum values of vectors R, RE, RQ, are determined for each variable of system (4): Ri = max j (Rj ),
(15)
REi = max j (RE j ), RQi = max j (RQ j )
The stepsize adjustment algorithm is presented in Inline Supplementary Fig. S1.
3.5. Representation of the reaction term The reaction term in the mass balance equations of system (4) was represented as a sum of two terms to extract a linear part with respect to the i-th component. This allows the
11
assigned linear part to be transferred to diagonal of the equations system, thus increasing the monotony of the system and securing non-negativity of the concentrations [20].
3.6. Distinctive features of the method of solution for a two-phase model In the case of a two-phase model the equations for the second phase are added to system (1). The combination of numerical methods remains the same. In some cases, when it is a lack of the terms describing convective and diffusive transfer in the second phase, it is possible to implement a special run algorithm (tridiagonal matrix algorithm or Thomas algorithm) securing the solution at step n+1 in both phases.
3.7. Distinctive features of the method of solution on the first scale The method of solution was described using a third scale model in cylindrical coordinates as an example. The transformation from the cylindrical to the spherical or Cartesian coordinates does not lead to a more sophisticated solution algorithm. The combination of numerical methods used for the solution of the equations of a second-scale model will be identical to that for the solution of equations of a third-scale model. The solution algorithm can be simplified only at modeling the processes on the first atomic-molecular scale when heat- and mass-transfer processes are not considered. For this scale the equations of the dynamics of adsorbed and intermediate species are put down as follows: d θ = f (θ), θ(0) = θ 0 dt
(16)
The system of equations (16) can be solved with an L–stable formula ensuring the second-order accuracy (7)-(10). The dimensionality of the system is equal to the number of adsorbed and intermediate species, Nθ. As a rule, for heterogeneous catalytic reactors this number is sufficiently small, at least in comparison with the system dimensionality when using the method of straight lines. Therefore, the dimensionality of a Jacobi matrix and matrix 12
Dn, which is equal to NθxNθ, is also small so it is possible to apply LU decomposition method to linear equations (9) and (10).
3.8. Distinctive features of the method of solution for a multiscale reactor model The solution of equations of a lower-scale is implemented in each higher scale grid point. The step size is chosen with account for changes in the solution on each scale.
4. Numerical examples The suggested strategy for multiscale EO process simulation was practically applied at the Boreskov Institute of Catalysis SB RAS to multiscale modeling of heterogeneous catalytic reactors. It was used to perform computations of the models of dimensionality ≤ 2D with a first-order derivative with respect to one of the coordinates on each studied scale of a multiscale reactor model. The introduced algorithm for the solution of a system of equations on each scale and its application to a multiscale model significantly saves the computation time and provides adequate accuracy and stability. Evidence to support these advantages can be found in the following examples of application of the proposed strategy to multiscale modeling of tubular reactors, fluidized bed reactor, and monolith catalytic reactor.
4.1. Multiscale mathematical modeling of tubular reactor The main features of the multiscale reactor model are as follows: considering mass transfer processes inside the catalyst taking into account pressure variation in the pellet (second scale); correlations for radial heat- and mass transfer parameters and wall heat transfer coefficients taking into account gas flow through holes in the grains and between them (third scale). The mathematical model is described in [21]. 4.1.1. Algorithm realization
13
The details of the algorithm realization in the case of tubular reactor can be found online in Inline Supplementary Algorithm Description S3. The programs are written in Intel(R) Visual Fortran Compiler 11.1. The time of computation of one variant ranges from 20 to 25 seconds on a Core i5-4570 class computer depending on catalytic reactions under consideration. The flowchart of the algorithm showing primarily the order of computational modules in the algorithm is given in Inline Supplementary Fig. S2. This realization of the algorithm makes it possible to obtain solutions at widely variable technological parameters of the process. During the realization of the algorithm it is not difficult to choose the initial approximation. With complex dependencies of reaction rates on the concentrations of substances it is possible to avoid shifting of the linear part to the diagonal of the system.
4.1.2. Modeling results of tubular reactors The conditions of modeling, the results of modeling using the above realizations of the algorithm and the results of modeling versus experimental data for some processes in tubular reactors are given in Table 1. From Table 1 it follows that the process parameters varied were inlet temperature, pressure, concentrations, linear velocity, cooling agent temperature as well as reactor dimensions, catalyst activity, shape and size of catalyst pellets. Modeling of methane steam reforming in a tubular reactor makes it possible to determine the optimal shape and size of catalyst pellets for typical technological conditions and feedstock compositions. Based on the results of modeling of oxidation of formaldehyde into formic acid in a tubular reactor recommendations on pilot plant design and optimum technological conditions for the process have been derived. For catalytic dehydration of ethanol to ethylene the technological
14
parameters ensuring maximum product yield at reduced catalyst loading have been determined. Validation of multiscale mathematical model and numerical algorithm was done by comparing the modeling results of these processes with pilot scale experimental data (Table 1).
15
Table 1. The results of multiscale mathematical modeling of tubular reactors. Conditions
Modeling results
Results of modeling versus experimental data
Methane steam reforming [21,22] Т in = 520 °C, Pin = 25 atm, Ltube = 12 m,
Optimum pellet sizes have been determined: Lcat =20 mm, Experimental data were obtained at a pilot plant
Dtube = 0.1 m, Т w = 880 °C, u = 9.9 m/s,
Lcat / Dcat = 1.5. Optimum pellet shapes:
in in = 24, CHin2O = 73, CCO = 1.5, CHin2 = CCH 4
in CO2
1.2, C
= 0.02, C
in N2
= 0.28 (vol.%).
Catalyst pellets: Raschig rings, cylinders, trilobe; Lcat / Dcat = 0.66–1.5
Dtube / max( Lcat , Dcat ) ≥ 5, N hole = 3–7. Catalyst pellet efficiency criteria:
1) Cylinder, N hole = 7 due to high η and λr ; 2) Trilobe.
designed by Haldor Topsoe, USA [23,24]. Relative error of temperatures along the axis and radius of the reactor and pressure drop along the bed were 3,4; 0,14 and 4,08 %, respectively [25].
The use of size-optimized catalyst pellets at a constant
In computations, the composition of products at
reactor performance allows decreasing Т w , which
the outlet of the tubular reactor is equal to
significantly reduces wall replacement costs, or to boost
experimental values.
reactor performance by decreasing residual methane and/or increasing natural gas consumption.
X CH 4 ( at l = 2.4 m) → max , ∆P → min
Oxidation of formaldehyde into formic acid [26]
u = 0.935–1.55 m/s, Dtube = 14–22 mm,
Resident time can be reduced more than twice at the same Experimental data were obtained at a pilot plant.
Т cool = 79–115 °C, Relative catalyst
X CH 2O of 98,5 % using two cascade reactors which differ
Relative errors of temperatures along the reactor
either by catalyst activity or Т cool . Based on modeling
axis, total X CH 2O and total SCH 2O2 were 2,8; 1,43
results, a pilot plant was designed and constructed.
and 1,27 %, especially in the case of using two
in activity, 0.25–1, CCH = 5.2–6.1 vol.%. 2O
reactors with different Т cool . Catalytic dehydration of ethanol to ethylene [22,27] 16
Dtube = 0.03 m, Ltube = 1.6–3.5 m, Pin = 1–
With an increase in Т w from 400 to 460 oC the
Experimental data were obtained at a pilot plant
3 atm, Т w = 400–460 oC, Residence time
residence time and, consequently, Ltube , required for
with a single-tube reactor. Relative errors of
= 2.1–3.9 s, u = 0.6–1.0 m/s, CCin2 H 5OH =
X C2 H 5OH = 98% are decreased. However, if Т w > 420
90.2, C
in H 2O
= 9.8 (vol.%).
conversion and ethylene selectivity were 0.4; 0.4
o
C YC 2 H 4 decreases because of the increase in
Catalyst pellets: Cylinder Lcat = 5 mm,
butylenes and other by-products formation. The
Dcat = 4 mm, Technological restrictions:
highest X C2 H 5OH = 98.7% and SC2 H 4 = 98.2% were
o X C2 H 5OH > 98%, SC2 H 4 > 97%, ∆P ≤ 0.2 observed at Т w = 420 C. The increase in u leads to the
atm, Optimization criterion: N tube → min.
temperatures along the reactor axis, ethanol
decrease in residence time which is explained by intensification of λr .
17
and 0.2 %, respectively.
As an example the process of methane steam reforming presented in Fig. 2 demonstrates three-dimensional dependencies of methane concentration (Fig. 2a) and temperature (Fig. 2b) on the reactor length and radius and also three-dimensional dependencies of methane concentration on reactor length (near the reactor wall) and catalyst radius in each point lengthwise (Fig. 2c). Kinetic equations and constants were taken from [28].
Figure 2. Three-dimensional dependencies of methane concentration (a) and temperature (b) on the reactor length and radius and three-dimensional dependencies of methane concentration on reactor length (near the reactor wall) and catalyst radius in each point lengthwise (c). Т in = in 793 K, Pin = 25 atm, Ltube = 3 m, Dtube = 0.1 m, u = 9.9 m/s, Т w = 1153 K, CCH = 24.0, 4 in in = 1.6, CHin2 = 1.3, CCO = 0.01, CNin2 = 3.49 (vol.%), Catalyst: cylinder Lcat = CHin2O = 69.6, CCO 2
0.019 m, Dcat = 0.014 m, N hole =4, Dhole = 0.003 m, Deq = 0.002 m.
The methane concentration and temperature profiles are shown in Figure 2. The program computes the concentrations of all substances (СН4, H2O, CO, H2, CO2, N2) across the reactor radius and length and also across the pellet radius. In addition to the above-listed processes the multiscale mathematical model and numerical algorithm were used to model a number of other processes in tubular reactors [22,29].
4.2. Multiscale mathematical modeling of monolith catalytic reactor
A monolith catalytic reactor with short contact time for partial oxidation of methane to syngas has been analysed by means of a dynamic one-dimensional two-phase reactor model with account for both transport limitations in the boundary layer of a fluid near the catalyst surface (third scale) and detailed kinetics via elementary reactions (first scale). The reaction mechanism contains 32 elementary steps describing methane oxidation on platinum, 14 gaseous compounds and 13 intermediates [22,30]. A detailed mathematical model is described in [30]. 4.2.1. Algorithm realization
The details of the algorithm realization in the case of monolith reactor can be found online in Inline Supplementary Algorithm Description S4. The programs are written in Intel(R) Visual Fortran Compiler 11.1. Time of computation of one variant is about a few hours on a Core i5-4570 class computer. Such a long time is appropriate for the non-steady state problem under consideration for which it takes a sufficiently long time to determine a solid-state thermal profile. A flowchart of the algorithm representing mainly the order of computational modules in the algorithm is shown in Inline Supplementary Fig. S3. The algorithm realization involves no difficulties. It is notable that during the warm-up process, when the values for intermediates undergo rapid changes, the size of the time step is
19
automatically set small. When the surface reaches a steady state while the monolith solid phase is slowly reaching its steady state the size of the time step is automatically increased.
4.2.2. Modeling results of monolith reactor
The program provides for varying the shape of the channels thereby modifying the dependences for calculating the coefficients of heat and mass exchange between the gaseous phase in the monolith channels and catalytic surface as well as the geometric sizes of the monolith channels and technological parameters of the process. The effects of such parameters as linear velocity (0.3 – 0.5 Nm/s), equivalent channel diameter (0.38 – 0.67 mm), inlet temperature (400 – 420 oC), inlet methane concentration (0.24 – 0.27 mole fraction in air) and effective thermal conductivity of the monolith (nearly 10 Jm-1s-1K-1) on the process dynamics, specifically on the ignition dynamics have been examined. It was shown that increasing linear velocity and equivalent channel diameter as well as decreasing axial conductivity of solid phase favors decreasing a time delay in syngas production in the Pt/Ce-Zr-La/α-Al2O3 honeycomb monolith with a triangular shape channels [22,30]. The examples of three-dimensional dependencies on time are presented in Fig. 3 for the fraction of surface hydrogen H* (Fig. 3a), concentration of hydrogen in gaseous phase of monolith channel (Fig. 3b) as well as temperature in gaseous phase (Fig. 3c) along the channel length.
20
Figure 3. Three-dimensional dependencies on time for the fraction of intermediate hydrogen H* (a), concentration of hydrogen in gaseous phase of monolith channel (b) and temperature in gaseous phase (c) along the channel length. Monolith with triangular channel cross-section. in = 25 mol.% in air, Т in = 400 oC, Pin = 1 atm, u = 0.3 m/s, deq = 0.67 mm, ε = 0.59. CCH 4
The concentration profiles for one of the thirteen intermediates and one of fourteen gaseous compounds and temperature profile in the monolith channel are presented in Figure 3. The program also performs computations for the remaining intermediates and gaseous compounds both in the monolith channel and on the monolith catalytic surface (in subsurface layer). The monolith solid phase temperature is also calculated. 21
4.2.3. Comparison of modeling results with experimental data
Simulation results are compared with transient experimental data on the reaction ignition in a full-size monolith 0.05 m in length. The time-dependent temperature profiles along the monolith channel and concentrations of substances at the monolith outlet show reasonable agreement between calculated and experimental values. However for the first 20 s the predicted temperature in the front part of the monolith increases slightly faster than the experimentally measured temperature. Both in the experiments and simulations during light-off, carbon dioxide, as a product of complete oxidation of methane, appears initially. Then, syngas selectivity slowly increases with rising temperature [22,30].
4.3. Multiscale mathematical modeling of fluidized bed reactor
Mathematical modelling of a fluidized bed reactor for propane-isobutane dehydrogenation has been performed on the base of a two-phase non-isothermal reactor model taking into account kinetics, convective transfer and axial gas-heat dispersion (third scale). Unique feature of the reactor is feeding of the regenerated hot catalyst into the top of the reactor, thereby increasing the reactor top temperature. The mathematical model is described in [31]. 4.3.1. Algorithm realization
The mass transfer equations are stationary one-dimensional second order differential equations in bubble and dense phases. Assuming the equal temperature in bubble and dense phases the heat transfer equation is a single-phase stationary one-dimensional second order differential equation. To apply the algorithm of Paragraph 3 let us make use of the method of transition to the non-steady state problem. Using the integro-interpolation method and the method of straight lines we get three systems of ODEs for concentrations in dense and bubble phases and for the temperature.
22
At each integration step the systems of ODEs are sequentially solved by method (7)-(10) using tridiagonal matrices and the stepsize adjustment algorithm. The programs are written in Intel(R) Visual Fortran Compiler 11.1. Time of computation of one variant is less than 30 seconds on a Core i5-4570 class computer. 4.3.2. Modeling results of fluidized bed reactor
Process simulation parameters are given in Table 2. An increase in the rate of heated catalyst introduction after regeneration resulted in a growth of average bed temperature ( Т av ). The effect of Т av and inlet propane mass fraction on the isobutylene and propylene yield is shown in Figure 4. Isobutylene yield was calculated as weight fraction in isobutane–isobutylene fraction.
Table 2. The parameters used in the simulation of industrial fluidized bed reactor Parameter Value Reactor diameter, m
4.5
Reactor height, m
4.9
Inlet gas velocity, Nm/s
0.15
Propane fraction in inlet mixture, %
0 – 60
Velocity of downward catalyst moving, m/s
0.01 – 0.1
Inlet catalyst temperature, oC
630
Fraction of bubble phase cross section
0.35
The isobutylene yield goes up along with increasing propane mass fraction and average layer temperature (Fig. 4). The rise of the inlet propane concentration leads to reduction
of
hydrogen
concentration
in
the
reaction
mixture
shifting
isobutane
dehydrogenation equilibrium to the right. The dependences of propylene yield (Fig. 4.) have more complicated behavior.
23
Isobutylene and propylene yield, % weight
70 60 50 40 30 20 0
15
30
45
60
Inlet propane fraction, % weight Figure 4. The effect of inlet propane mass fraction and Т av on the isobutylene (solid lines) and propylene (dashed lines) yield. Black lines (square symbols) correspond to Т av = 560 oC, red lines (circle symbols) correspond to Т av = 585 oC.
4.3.3. Comparison of modeling results with experimental data
The model and numerical algorithm are validated against experimental data obtained for commercial fluid bed isobutane dehydrogenation reactor [22,31]. From Table 3, it can be seen that the model is able to reproduce satisfactorily the dependence of isobutylene yield and selectivity on the average bed temperature.
Table 3. An influence of the average bed temperature on isobutylene selectivity and yield. Modelling results and experimental data. u = 0.15 m/s, Ciin−C4 H10 = 95.5 wt%. Average bed
Isobutylene selectivity, % weight
temperature ,
Experimental
Numerical
Experimental
Numerical
data
results
data
results
o
C
Isobutylene yield, % weight
585
83.2
83.5
38.5
38.0
575
84.2
85.0
34.0
35.0
24
554
85.2
86.0
31.0
30.0
4.4. A comparison of modeling approaches
To compare computation time for the presented approach and conventional approaches let us consider, for instance, the application of the straight line method in their realization. Then in case of two-scale process modeling with model dimensionality ≤ 2D on each scale, for instance, on reactor and catalyst pellet scales, the presented approach ensures a significant decrease in system dimensionality, thereby a decrease in computation time (Table 4).
Table 4. Comparison of computation time for different approaches. Approaches
Сonstraint
Conventional Second-order derivatives approach.
Order of solution
Computation time
Simultaneous solution of
f(N·NL·NP)
with respect to all coordinates a system of equations on on each scale are normal.
all scales.
Presented
First-order derivative with
Separate solution of a
approach.
respect to a coordinate on
system of equations on
each scale.
each scale.
f(NP) ·N+ f(N)
Besides, to solve the resultant systems of linear equations of type (9) traditionally either the LU decomposition or Gaussian elimination method is used, otherwise, if the system is represented in form (5), the tridiagonal matrix algorithm is used concurrently with the iteration process. In these cases the computation time is proportional either to N3eqs or 8Neqs·Niter, respectively. In our case, following the representation of the system in form (5) the tridiagonal matrix algorithm is used and there is no need to proceed with iterations if one uses the stepsize adjustment algorithm. The computation time is then proportional to 8Neqs.
25
5. Conclusions
EO approach-based multiscale models of heterogeneous catalytic reactors, which are noteworthy as the degree of complexity of such models is dictated only by the necessity, which is sufficient too, may have dimensionality ≤ 2D with the first-order derivative with respect to a coordinate on each scale for a large class of reactors, including nonstandard ones. Simultaneous solution of a system of equations, including all the scales under consideration, which is characterized by high dimensionality, is extremely time-consuming because bulky data storage is required. A strategy of modeling such systems has been suggested based on the use of an identical set of numerical methods, independently on each scale, with step size control implemented with account for the rate of change of the variables on each scale. The introduced algorithm of numerical solution significantly saves computation time and provides adequate accuracy and stability. The developed strategy and the combination of numerical methods suggested for the solution of equations for the model on each scale proved effective for multiscale modeling of heterogeneous catalytic reactors such as tubular reactors, monolith catalytic reactor, and fluidized bed reactor. When switching between the reactor types the parameters of computation algorithms do not need to be corrected. The developed strategy can also be applied to multiscale modeling of heterogeneous catalytic reactors of other types, including nonstandard ones. The models used for modeling of such reactors agree with the modern trends in chemical reaction engineering, such as [8,9,32]: an estimation of transfer coefficients with correlations capturing also the effect on reactor scale; simultaneous modeling of transport-kinetic interactions on a catalytic particle and reactor scales; validation of the models against more than one output parameters and the data obtained under representative current process operating conditions.
26
Acknowledgments.
The author highly acknowledges her colleagues Ludmilla Bobrova, Victor Chumachenko, Nataliya Chumakova, Aigana Kagyrmanova, Vitalii Kashkin, Elena Ovchinnikova, Ilya Zolotarskii (Boreskov Institute of Catalysis SB RAS) for collaboration as well as for fruitful discussions of the heterogeneous catalytic reactors modeling strategy and results presented in this paper. This work was conducted within the framework of budget project No. 0303-2016-0017 for Boreskov Institute of Catalysis.
Nomenclature ai , pi , γij , βij – parameters (Eq. 6);
A, B, C , f (u ) – coefficients in equations (Fig. 1); Ai – tridiagonal matrix depend on the solution in a non-linear manner (Eq. 4); A – tridiagonal matrix of dimensionality NxN (Eq. 5); Ci – vector depend on the solution in a non-linear manner (Eq. 4) or concentration of i-th
compound in gas phase, (mol.%, vol.%, wt.%); C – vector of dimensionality N (Eq. 5); Dn – matrix (Eqs 8-10); D – diameter, (m); Deq – equivalent catalyst pellet diameter as a plate, (m);
deq – equivalent channel diameter, (mm); fu – Jacobi matrix; f – function in Table 4;
27
I – unity matrix; h – stepsize; l – the axial coordinate of the reactor, (m); L – length of the reactor (m);
M – number of variables in Eq. 1; N – number of grid points along the reactor radius; NL – number of grid points along the reactor length; NP – number of grid points along the pellet radius; N text – number (see Subscripts);
P – pressure, (atm); ∆P – pressure drop along the reactor, (atm); R – radius of reactor, (m); S – selectivity, (%); Т – gas temperature, (oC, K); t – time, (s);
u – linear velocity under normal temperature, (m s-1);
ui – variable in Eq. 1; X – conversion, (%);
x – the radial coordinate of the reactor, (m). Y – yield of product, (%); y – the axial coordinate of the reactor, (m).
Greek letters:
α – the heat-transfer coefficient from the tube wall to the bed, (J m-2s-1K-1); ε – porosity or accuracy of the algorithm in Eq. 11;
η – effectiveness factor of catalyst pellet; 28
λr – radial heat conductivity in the catalyst bed, (J m-1s-1K-1); θ – intermediate species, (undimensional); Subscripts: 0
– initial value;
av
– average;
cat
– catalyst;
cool
– cooling agent;
eqs
– equations;
hole
– catalyst hole;
– variable;
i
in
– inlet; – iteration;
iter
j
– point along reactor radius;
tube
w
– tubular reactor;
– tube wall;
References
[1] M. G. Slin’ko, History of the development of mathematical modeling of catalytic processes and reactors, Theoretic. Found. Chem. Eng. 41 (1) (2007) 13–29. [2] J. J. Lerou, K. M. Ng, Chemical reaction engineering: A multiscale approach to a multiobjective task, Chem. Eng. Sci. 51 (10) (1996) 1595–1614. [3] M. G. Slin’ko, Modelirovanie khimicheskikh reaktorov (Modeling of Chemical Reactors), Nauka: Novosibirsk, 1968 (in Russian). [4] W. Marquardt, Trends in computer-aided process modeling, Computers Chem. Engng 20 (1996) 591-609.
29
[5] R. Bogusch, B. Lohmann, W. Marquardt, Computer-aided process modeling with MODKIT, Computers Chem. Engng 25 (2001) 963-995. [6] A. S. Moharir, S. S. Shah, R. D. Gudi, et al, Generalized reactor model: an object oriented approach to reactor modeling, Europ.Symp. on Comp. Aided Process Engng - 11 (2001) 243248. [7] R. Aris, The optimal design of chemical reactors; a study in dynamic programming, New York, Academic Press, 1961. [8] H. Stitt, M. Marigo, S. Wilkinson and T. Dixon, How good is your model?, Johnson Matthey Technol. Rev. 59 (2) (2015) 74-89. [9] M. P. Dudukovic, Reaction engineering: Status and future challenges, Chem. Eng. Sci. 65 (2010) 3-11. [10] F. Zhao, X. Chen, A hybrid numerical-symbolic solving strategy for equation-oriented process simulation and optimization, AIChE J. DOI 10.1002/aic.15622. [11] W. P. M. van Swaaij, A. G. J. van der Ham, A. E. Kronberg, Evolution patterns and family relations in G–S reactors, Chem. Eng. J. 90 (2002) 25–45. [12] L. H. Roblee, R. M. Baird, J. W. Tierney, Radial porosity variations in packed beds, AIChE J. 4 (1958) 460-464. [13] I. I. Ioffe, L. M. Pismen, Engineering Chemistry of Heterogeneous Catalysis, Chimiya, Leningrad, 2nd ed., 1972 (in Russian). [14] A. A. Samarskii, A. P. Mikhailov, Matematicheskoe modelirovanie. Idei. Metody. Primery (Mathematical Modeling: Ideas, Methods, and Examples), Moscow: Fizmatlit, 2001. [15] C. A. J. Fletcher, Computational techniques for fluid dynamics, vol. 1, Springer-Verlag, 1991. [16] H. H. Rosenbrock, General implicit processes for the numerical solution of differential equations, Computer J. 5 (1965) 329-330.
30
[17] E.A. Novikov, Numerical methods for solution of differential equations in chemical kinetics, in: V. I. Bykov (Ed.), Mathematical methods in chemical kinetics, Nauka, Novosibirsk, 1990, 53-68 (in Russian). [18] V. N. Bibin, In Proceedings of XII International Conference on Chemical Reactors, Boreskov Institute of Catalysis, Novosibirsk, 1 (1996) 127–132 (in Russian). [19] N. V. Vernikovskaya, T. L. Pavlova, N. A. Chumakova, A. S. Noskov, Mathematical modelling of filtration and catalytic oxidation of diesel particulates in filter porous media, in: J.C. Cortes, L. Jodar, R. Villanueva (Eds.), Mathematical Modeling in Social Sciences and engineering, Nova Sciences Publishers Inc., Hauppauge, NY, 2014 27–40. [20] L. A. Rapatskii, L. V. Yausheva, Mathematical Modeling of Catalytic Processes in Porous Media. Preprint of Computation Center, Siberian Div., Russ. Acad. Sci.: Novosibirsk 1019 1994 (in Russian). [21] A. P. Kagyrmanova, I.A. Zolotarskii, E. I. Smirnov, et al, Optimum dimensions of shaped steam reforming catalysts, Chem. Eng. J. 134 (2007) 228-234. [22] L. N. Bobrova, N. V. Vernikovskaya, A. P. Kagyrmanova, V. N. Kashkin, I. A. Zolotarskii, Application of Mathematical Modeling and Simulation to Study the Complexity of the Basic Chemical and Physical Phenomena in Catalytic Reactors, in Book: Mathematical Modelling. Editors: Christopher R. Brennan. Nova Science Publishers 2012, 475-577. [23] J. R. Rostrup-Nielsen, Catalytic Steam Reforming, Berlin: Springer, 1984. [24] J. R. Rostrup-Nielsen, L. J. Christiansen, J. H. B. Hansen, Activity of Steam Reforming Catalysts: Role and Assessment, Applied Catalysis 43 (1988) 287-303. [25] A. P. Kagyrmanova, Optimization of catalyst shape and size in tubular packed bed reactor, PhD Thesis, Novosibirsk, 2009.
31
[26] I. A. Zolotarskii, T. V. Andrushkevich, G. Ya. Popova, et al, Modeling, design and operation of pilot plant for two-stage oxidation of methanol into formic acid, Chem. Eng. J. 238 (2014) 111-119. [27] A. P.
Kagyrmanova,
V. A.
Chumachenko,
V. N.
Korotkikh, et al, Catalytic
dehydration of bioethanol to ethylene: Pilot-scale studies and process simulation, Chem. Eng. J. 176 –177 (2011) 188 –194. [28] J. Xu, G. F. Froment, Methane steam reforming: II. Diffusional limitations and reactor simulation, AIChE J. 35 (1989) 97. [29] E. V. Ovchinnikova, N. V. Vernikovskaya, T. V. Andrushkevich, V. A. Chumachenko, Mathematical modeling of β-picoline oxidation to nicotinic acid in multitubular reactor: Effect of the gas recycle, Chem. Eng. J. 176-177 (2011) 114-123. [30] N. V. Vernikovskaya, L. N. Bobrova, L. G. Pinaeva, et al, Transient behavior of the methane partial oxidation in a short contact time reactor: Modeling on the base of catalyst detailed chemistry, Chem. Eng. J. 134 (2007) 180-189. [31] N. V. Vernikovskaya, I. G. Savin, V. N. Kashkin, et al, Dehydrogenation of propane– isobutane mixture in a fluidized bed reactor over Cr2O3/Al2O3 catalyst: Experimental studies and mathematical modelling, Chem. Eng. J. 176-177 (2011) 158-164. [32] M. P. Dudukovic, Trends in catalytic reaction engineering, Catal. Today 48 (1999) 5-15.
Captions:
Figure 1. Schematic representation of the structure of partial differential equations: (a) – with respect to space coordinate “y” contains only first-order derivatives, (b) – the first-order derivative with respect to time.
32
Figure 2. Three-dimensional dependencies of methane concentration (a) and temperature (b) on the reactor length and radius and three-dimensional dependencies of methane concentration on reactor length (near the reactor wall) and catalyst radius in each point lengthwise (c). Т in = in 793 K, Pin = 25 atm, Ltube = 3 m, Dtube = 0.1 m, u = 9.9 m/s, Т w = 1153 K, CCH = 24.0, 4
in in = 1.6, CHin2 = 1.3, CCO = 0.01, CNin2 = 3.49 (vol.%), Catalyst: cylinder Lcat = CHin2O = 69.6, CCO 2
0.019 m, Dcat = 0.014 m, N hole =4, Dhole = 0.003 m, Deq = 0.002 m.
Figure 3. Three-dimensional dependencies on time for the fraction of intermediate hydrogen H* (a), concentration of hydrogen in gaseous phase of monolith channel (b) and temperature in gaseous phase (c) along the channel length. Monolith with triangular channel cross-section. in CCH = 25 mol.% in air, Т in = 400 oC, Pin = 1 atm, u = 0.3 m/s, deq = 0.67 mm, ε = 0.59. 4
Figure 4. The effect of inlet propane mass fraction and Т av on the isobutylene (solid lines) and propylene (dashed lines) yield. Black lines (square symbols) correspond to Т av = 560 oC, red lines (circle symbols) correspond to Т av = 585 oC.
33
A novel strategy is proposed for multiscale modeling of catalytic reactors. This strategy was applied for modeling tubular, monolith, fluidized bed reactors. Robustness of the algorithm is shown for multiscale modeling of these reactors.
34
S1385-8947(17)30844-6 http://dx.doi.org/10.1016/j.cej.2017.05.095 CEJ 16993
To appear in:
Chemical Engineering Journal
Please cite this article as: N.V. Vernikovskaya, An equation-oriented approach to modeling heterogeneous catalytic reactors, Chemical Engineering Journal (2017), doi: http://dx.doi.org/10.1016/j.cej.2017.05.095
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
AN EQUATION-ORIENTED APPROACH TO MODELING HETEROGENEOUS CATALYTIC REACTORS
N. V. Vernikovskaya
Boreskov Institute of Catalysis SB RAS, Novosibirsk 630090, Pr. Ak. Lavrentieva 5, Russia Novosibirsk State University, Novosibirsk, Russia Novosibirsk State Technical University, Novosibirsk, Russia E-mail [email protected]
Abstract At multiscale modeling of heterogeneous catalytic reactors there appear systems comprising a large number of equations. The solution of such systems is an arduous task and available algorithms require voluminous computations. A modeling strategy has been suggested for such systems. An effective solution algorithm has been developed for a large class of models based on the application of an identical set of numerical tools such as integro-interpolation method, method of straight lines, a special case of a second-order Rosenbrock method, tridiagonal matrix algorithm or Thomas algorithm on each scale of a multiscale reactor model. Step size control is implemented with account for the rate of change of the variables on each scale. The efficiency and robustness of this algorithm were demonstrated at multiscale modeling of heterogeneous catalytic reactors such as tubular reactors, monolith catalytic reactor, and fluidized bed reactor. When switching between reactor types at multiscale modeling it is not necessary to modify the algorithm. The algorithm can also be used for multiscale modeling of heterogeneous catalytic reactors of other types.
Keywords: Catalytic reactors, multiscale modeling, numerical methods.
1. Introduction Mathematical modeling of diverse chemical processes running in various reactors can help to reduce experimental efforts during chemical process development. In chemical reactor modeling chemical reactors are represented as space-time hierarchical structures. The space scale ranges from 10-11 to 10 3 m, the associated time scale ranges from 10 -15 to 108 s [1,2]. The hierarchical multiscale approach is a backbone of the mathematical modeling of chemical reactors [1,3]. This approach consists in decomposing a complex chemical-technological process into chemical and physical components, studying these components independently and subsequently carrying out the synthesis of a general mathematical model from the models of separate parts of a complex process. Atomic-molecular processes and heat- and mass transfer in porous catalyst pellets or in the catalytic layer represent the lower two scales of the mathematical models of heterogeneous catalytic reactors. The processes of heat- and masstransfer in the reactor, i.e. in the catalyst bed or monolith channel, represent the third scale, and modeling of a switch device with regard to the processes of mixing, heat exchange, etc – the forth scale. A multiscale mathematical model can describe the preceding levels and given level simultaneously. The history of development and application of the hierarchical approach to the modeling of catalytic processes is described in detail in [1]. Recently there have been many studies on the multiscale modeling of catalytic processes in reactors of various types. Generally, to treat a developed mathematical model there are two groups of mathematical modeling tools based on two different approaches: block-oriented (BO) and equation-oriented (EO) [4]. The BO approaches address modeling on the flowsheet or reactor levels [5,6]. In the first case, every process is abstracted by a flowsheet consisting of the coupled units, e.g. heatexchanger, compressor, reactor, absorber, etc. However for the reactor most simulators
2
provide yield based reactor models or a set of idealized models like continuous stirred tank reactor (CSTR), plug flow, etc. with limited capabilities for defining reactor kinetics [6]. In the second case, every reactor model is decomposed into simpler objects such as reaction kinetics, thermodynamics, reactor geometry, reactor configuration, etc. Data are exchanged between the individual objects. BO approaches have many advantages. However reactor configurations are limited to a fixed set of reactors with given models: adiabatic bed reactor, tubular reactor, fluidized bed reactor, and every unit or object is solved one by one in a predefined sequence. This approach is basically used for modeling processes on the same scale. Besides, BO approaches do not support the implementation of more detailed models of standard reactors or models of non-standard reactors: monolith reactor, radial flow reactor, moving bed, multiphase reactors, membrane processes etc. EO approaches have some advantages over BO approaches, for they support the implementation of more complex reactor models and the simulators of EO process solve the system of equations simultaneously. The last-mentioned advantage indicates that EO approaches are more flexible as compared to BO approaches and, therefore, are more suitable for multiscale modeling of heterogeneous catalytic reactors. Among the techniques within equation-oriented approaches to reactor modeling one can distinguish computational fluid dynamics (CFD). CFD (commercial or in-house software) provides an excessive data set with a limited physical meaning, so CFD flow models can be validated only with a great deal of experimental data. This technique is too time-consuming. Another technique within EO approaches is using the models of necessary and sufficient degree of sophistication. One cannot always improve a model by adding parameters. These models satisfy Occam’s razor principle, which means the models “contain the elements that are needed, but without extra embellishments” [7]. Therefore, compared to CFD-based models EO approach-based models use less parameters and can be measured or estimated
3
independently. Validation of such models requires less experimental data. Transfer coefficients can be obtained from the correlations which capture also the effect of reactor scale. The models and simulations that accurately treat both the physics and chemistry of real systems and the parameters independently obtained from measured or estimated data sets are currently considered to be the most appropriate approach to modeling various chemical reactors [8,9]. In the case of multiscale modeling of heterogeneous catalytic reactors the EO approach-based technique allows constructing models of dimensionality ≤ 2D with the firstorder derivative with respect to a coordinate on each scale for a large class of reactors. Simultaneous solution of a system of equations as a whole simultaneously on all scales under consideration is an advantage of EO approaches, but the solution cost increases significantly as the size of the model increases [10]. In this study a novel strategy is proposed for multiscale EO process simulation. This strategy is based on the use of a single set of numerical methods on each scale. This strategy is applicable to models of dimensionality ≤ 2D with the first-order derivative with respect to a coordinate on each scale under consideration.
2. Field of applications of the proposed strategy In the case of multiscale modeling of heterogeneous catalytic reactors the use of the models with the elements that are needed, but without extra embellishments for a large class of reactors permits the model to be presented on each of the scales under consideration as a system of differential equations of dimensionality ≤ 2D with the first-order derivative with respect to a coordinate. This means that the system of partial differential equations can contain both first- and second-order derivatives with respect to one of the coordinates and only first-order derivatives with respect to the other coordinate (Fig. 1 a) or time (Fig. 1b).
4
All types of heterogeneous catalytic gas-solid reactors can be subdivided into three groups: packed bed, fluidized bed, and barrier wall reactors [11].
2.1. Processes on the third (reactor) scale The class of packed bed reactors includes two specific subgroups: tubular (1) and adiabatic (2) reactors. For steady-state processes in the case (1) it is necessary to take account of radial convection and diffusion heat- and mass transfer processes, whereas along the reactor length, which is normally several times larger than the radius, it is sufficient to take account of only the convection transfer, i.e. first-order derivative with respect to length. In the case (2) for steady-state processes it is not necessary to take account of radial heat- and masstransfer processes, while it is possible to consider axial heat- and mass transfer with any degree of sophistication, i.e. by taking account of both convection and diffusion components. In adiabatic reactors radial heat- and mass transfer can be neglected because there is usually no heat supply or removal through the reactor wall and, as a result, no gradient of temperature or substances concentrations across the reactor radius. As for the effects of reactor wall on layer porosity, the layer porosity can fluctuate within the distance of several particle diameters from the reactor wall for regularly packed beds and only one diameter for randomly packed beds, which is typical of all reactor types [12,13]. In contrast to tubular reactors, in the adiabatic ones the reactor to particle diameter ratios are normally very large so the effect of porosity fluctuations on the interstitial gas velocity and, consequently, the radial heat- and mass transfer can be neglected. To obtain the first-order derivative, in the model it is possible to use the method of transition to the non-steady state problem. Furthermore, in this case it is also possible to solve the non-steady state problem. This can be used in the description of the process in the monolith reactor.
5
(a)
(b)
Figure 1. Schematic representation of the structure of partial differential equations: (a) – with respect to space coordinate “y” contains only first-order derivatives, (b) – the first-order derivative with respect to time.
For steady-state processes in fluidized bed reactors it is not necessary to take account of the radial heat- and mass transfer, whereas it is possible to consider axial heat- and mass transfer with any degree of sophistication, i.e. by taking account of both convection and diffusion components. In fluidized bed reactor radial heat- and mass transfer can be neglected because of intensive catalyst mixing. To obtain the first-order derivative in the model it is possible to use the method of transition to the non-steady state problem and also to solve the non-steady state problem. Barrier wall reactors can be characterised by a barrier that incorporates a catalyst placed between two gas spaces or chambers with different gas composition [11]. Within this barrier the chemical reaction takes place and product flows may be manipulated to go to one or both sides using selective permeability or pressure gradient. In the reactors of this type for steadystate processes it is necessary to take account of the processes of radial convective and diffusive heat- and mass transfer, including the barrier, whereas along the reactor’s length it sufficient to consider the convective transfer, i.e. the first-order derivative with respect to length.
6
2.2. Description of processes on the second (catalyst) scale In the foregoing descriptions of reactor scale processes the mathematical models for catalyst pellets or catalytic layer, if necessary, can be simplified as one-dimensional. In this case to determine the first-order derivative in the model it is possible to use the method of transition to the non-steady state for a steady-state reactor model or to solve the non-steady state problem for a non-steady state reactor model.
2.3. Description of processes on the first (atomic-molecular) scale On this scale one considers either a detailed molecular nonsteady-state kinetic model for surface reactions and gets the first-order derivative with respect to time or a detailed molecular steady-state kinetic model for surface reactions, and in the latter case it is possible to use the method of transition to the non-steady state problem. A multiscale reactor model may incorporate the models of all three scales or the models of only two scales. In any case the equations of a lower-scale model are to be solved in each higher scale grid point. The combination of numerical methods suggested to solve the equations is identical for the models on every scale. This combination can be simplified only for the model of the first-scale processes, in which heat- and mass transfer processes are not considered. The suggested combination of numerical methods can be used in modeling the process for a single-scale model.
3. Solution method Without the loss of generality, put down the equations for a third-scale model in cylindrical coordinates. For clearness and simplification of the algorithm description, assume that radial and axial transfer velocities are constant and can be taken outside the derivative:
7
∂ui ∂y
+ ai ⋅
∂u 1 ∂(x ⋅ ui ) 1 ∂ bi ⋅ x ⋅ i = ci , i = 1, M ⋅ − ⋅ x ∂x x ∂x ∂x
(1)
Initial conditions:
y = 0 : ui = ui0
(2)
Boundary conditions:
∂ui
x =0: x =R:
∂x ∂ui ∂x
bM ⋅
= 0; i = 1, M ; = 0; i = 1, M − 1; ∂uM
∂x
(3)
= α ⋅ (Tw − uM )
where u = (u1,..., uM ) is the vector of the concentrations of substances and temperature,
ai , bi , ci are the coefficients depending nonlinearly on the solution, i.e. on the concentrations of substances and temperature, α is the heat-transfer coefficient from the tube wall to the bed, J/(m2.s.K), Tw is tube wall temperature, K.
3.1. Integro-interpolation method Discrete analogues of the model equations were constructed by using the integrointerpolation method. Transition to discrete models cannot be made in a purely formal way. Such a transition must leave in force as many basic properties of the source object models as possible [14]. In respect to the equations of system (1) with initial (2) and boundary (3) conditions such a transition means that for the corresponding difference scheme it is necessary that the discrete analog of the mass and energy conservation laws be fulfilled. A common approach to the construction of such difference scheme is based on writing this law in an integral form for cells x j −1 ≤ x < x j , x j = j ⋅ h j , j = 1, N of the difference grid 8
followed by substitution of the obtained integrals and derivatives with approximate difference equations, step h j can be obtained from the condition of equal cylinder section areas. A complete description of integro-interpolation method can be found online in Inline Supplementary Algorithm Description S1.
3.2. Method of straight lines We do not introduce discrete formulas for derivatives with respect to “y”, thus transforming partial differential equation (PDE) for each ui into a set of ordinary differential equations (ODEs) for nodal values along the reactor radius. The process of transformation of equation (1) into system of ODEs for each variable ui is an example of application of the method of straight lines or the method of semi-discretization [15]. By successively applying this procedure to all inner nodes and taking account of the boundary conditions (3) in the relation at j=1 and j=N, one gets a system of equations for each variable ui = (u1,..., uM ) written in the form: dui = A ⋅ u + C , ∀i = 1, M dy i i i
(4)
where three diagonal matrix Ai and vector Сi depend on the solution in a non-linear manner. The method of straight lines is advantageous because to solve the semidiscrete form of initial differential equation in partial derivatives it is possible to use different methods for ordinary differential equations. Let us put down system (4) for one variable ui , for simplicity, skipping index “i”: du = f (u ), f (u ) = A ⋅ u + C , u(0) = u 0 dy
(5)
9
where A is a tridiagonal matrix of dimensionality NxN, С is a vector of dimensionality N which are nonlinearly dependable on the solution, u, f(u) are vectors of dimensionality N.
3.3. A special case of a second-order Rosenbrock method When solving the problems with initial data described by ordinary differential equations the best results are normally obtained with either linear multistep or Runge-Kutta methods. In 1963, Rosenbrock [16] proposed the methods that differ from the explicit Runge-Kutta schemes by the regularization of the right part of a differential problem. Such methods can be derived from a class of semi-implicit methods if each kni is computed with the use of not more than one Newton’s iteration. Later on the methods of the specified type were called Rosenbrock-type methods [17]. They can be specified as follows: m
un +1 = un + ∑ pi ⋅ kni , I − ai ⋅ h ⋅ fu
i =1
i −1 ⋅ un + ∑ γij ⋅ knj ⋅ kni = h ⋅ f j =1
i −1 ⋅ un + ∑ βij ⋅ knj , j =1
(6)
where I is a unity matrix, fu is Jacobi matrix of the right hand size of system (5), ai , pi , γij , βij are the parameters. For Rosenbrock-type methods to be implemented most effectively all ai must be mutually equal, and γij = 0. In this case at each step it is necessary to compute and invert only one matrix of dimensionality N. To solve the system of equations (5) let us make use of a second-order L-stable formula as follows [17]: un +1 = un + (0.5 + 0.25 2) ⋅ kn1 + (0.5 − 0.25 2) ⋅ kn 2
(7)
Dn = I − (1 − 0.5 2) ⋅ h ⋅ fu (un )
(8)
Dn ⋅ kn 1 = h ⋅ f (un )
(9)
Dn ⋅ kn 2 = h ⋅ f (un + 2 ⋅ kn1 )
(10)
10
In this case the Jacobi matrix is of dimensionality NxN and systems of linear equations (9) and (10) can be solved using some laborious algorithms such as LU decomposition with the number of operations proportional to N3. As applicable to our system of equations (5), in which matrix A is tridiagonal, it may be preferable to use the tridiagonal matrix algorithm concurrently with the iteration process. A complete description of this procedure can be found online in Inline Supplementary Algorithm Description S2.
3.4. Step size control Step h is chosen so that for the global error the following relation is fulfilled [17-19]: (11)
kn 2 − kn1 / un < 15 ⋅ ε
The stepsize adjustment algorithm is based on the computation of the values of the right part and difference between kn1 and kn2 for each equation of system (7): R = un + (0.5 + 0.25 2) ⋅ kn1 + (0.5 − 0.25 2) ⋅ kn 2
(12)
RE = kn 1 − kn 2
(13)
RQ = kn 1 − kn 2 / un
(14)
Then maximum values of vectors R, RE, RQ, are determined for each variable of system (4): Ri = max j (Rj ),
(15)
REi = max j (RE j ), RQi = max j (RQ j )
The stepsize adjustment algorithm is presented in Inline Supplementary Fig. S1.
3.5. Representation of the reaction term The reaction term in the mass balance equations of system (4) was represented as a sum of two terms to extract a linear part with respect to the i-th component. This allows the
11
assigned linear part to be transferred to diagonal of the equations system, thus increasing the monotony of the system and securing non-negativity of the concentrations [20].
3.6. Distinctive features of the method of solution for a two-phase model In the case of a two-phase model the equations for the second phase are added to system (1). The combination of numerical methods remains the same. In some cases, when it is a lack of the terms describing convective and diffusive transfer in the second phase, it is possible to implement a special run algorithm (tridiagonal matrix algorithm or Thomas algorithm) securing the solution at step n+1 in both phases.
3.7. Distinctive features of the method of solution on the first scale The method of solution was described using a third scale model in cylindrical coordinates as an example. The transformation from the cylindrical to the spherical or Cartesian coordinates does not lead to a more sophisticated solution algorithm. The combination of numerical methods used for the solution of the equations of a second-scale model will be identical to that for the solution of equations of a third-scale model. The solution algorithm can be simplified only at modeling the processes on the first atomic-molecular scale when heat- and mass-transfer processes are not considered. For this scale the equations of the dynamics of adsorbed and intermediate species are put down as follows: d θ = f (θ), θ(0) = θ 0 dt
(16)
The system of equations (16) can be solved with an L–stable formula ensuring the second-order accuracy (7)-(10). The dimensionality of the system is equal to the number of adsorbed and intermediate species, Nθ. As a rule, for heterogeneous catalytic reactors this number is sufficiently small, at least in comparison with the system dimensionality when using the method of straight lines. Therefore, the dimensionality of a Jacobi matrix and matrix 12
Dn, which is equal to NθxNθ, is also small so it is possible to apply LU decomposition method to linear equations (9) and (10).
3.8. Distinctive features of the method of solution for a multiscale reactor model The solution of equations of a lower-scale is implemented in each higher scale grid point. The step size is chosen with account for changes in the solution on each scale.
4. Numerical examples The suggested strategy for multiscale EO process simulation was practically applied at the Boreskov Institute of Catalysis SB RAS to multiscale modeling of heterogeneous catalytic reactors. It was used to perform computations of the models of dimensionality ≤ 2D with a first-order derivative with respect to one of the coordinates on each studied scale of a multiscale reactor model. The introduced algorithm for the solution of a system of equations on each scale and its application to a multiscale model significantly saves the computation time and provides adequate accuracy and stability. Evidence to support these advantages can be found in the following examples of application of the proposed strategy to multiscale modeling of tubular reactors, fluidized bed reactor, and monolith catalytic reactor.
4.1. Multiscale mathematical modeling of tubular reactor The main features of the multiscale reactor model are as follows: considering mass transfer processes inside the catalyst taking into account pressure variation in the pellet (second scale); correlations for radial heat- and mass transfer parameters and wall heat transfer coefficients taking into account gas flow through holes in the grains and between them (third scale). The mathematical model is described in [21]. 4.1.1. Algorithm realization
13
The details of the algorithm realization in the case of tubular reactor can be found online in Inline Supplementary Algorithm Description S3. The programs are written in Intel(R) Visual Fortran Compiler 11.1. The time of computation of one variant ranges from 20 to 25 seconds on a Core i5-4570 class computer depending on catalytic reactions under consideration. The flowchart of the algorithm showing primarily the order of computational modules in the algorithm is given in Inline Supplementary Fig. S2. This realization of the algorithm makes it possible to obtain solutions at widely variable technological parameters of the process. During the realization of the algorithm it is not difficult to choose the initial approximation. With complex dependencies of reaction rates on the concentrations of substances it is possible to avoid shifting of the linear part to the diagonal of the system.
4.1.2. Modeling results of tubular reactors The conditions of modeling, the results of modeling using the above realizations of the algorithm and the results of modeling versus experimental data for some processes in tubular reactors are given in Table 1. From Table 1 it follows that the process parameters varied were inlet temperature, pressure, concentrations, linear velocity, cooling agent temperature as well as reactor dimensions, catalyst activity, shape and size of catalyst pellets. Modeling of methane steam reforming in a tubular reactor makes it possible to determine the optimal shape and size of catalyst pellets for typical technological conditions and feedstock compositions. Based on the results of modeling of oxidation of formaldehyde into formic acid in a tubular reactor recommendations on pilot plant design and optimum technological conditions for the process have been derived. For catalytic dehydration of ethanol to ethylene the technological
14
parameters ensuring maximum product yield at reduced catalyst loading have been determined. Validation of multiscale mathematical model and numerical algorithm was done by comparing the modeling results of these processes with pilot scale experimental data (Table 1).
15
Table 1. The results of multiscale mathematical modeling of tubular reactors. Conditions
Modeling results
Results of modeling versus experimental data
Methane steam reforming [21,22] Т in = 520 °C, Pin = 25 atm, Ltube = 12 m,
Optimum pellet sizes have been determined: Lcat =20 mm, Experimental data were obtained at a pilot plant
Dtube = 0.1 m, Т w = 880 °C, u = 9.9 m/s,
Lcat / Dcat = 1.5. Optimum pellet shapes:
in in = 24, CHin2O = 73, CCO = 1.5, CHin2 = CCH 4
in CO2
1.2, C
= 0.02, C
in N2
= 0.28 (vol.%).
Catalyst pellets: Raschig rings, cylinders, trilobe; Lcat / Dcat = 0.66–1.5
Dtube / max( Lcat , Dcat ) ≥ 5, N hole = 3–7. Catalyst pellet efficiency criteria:
1) Cylinder, N hole = 7 due to high η and λr ; 2) Trilobe.
designed by Haldor Topsoe, USA [23,24]. Relative error of temperatures along the axis and radius of the reactor and pressure drop along the bed were 3,4; 0,14 and 4,08 %, respectively [25].
The use of size-optimized catalyst pellets at a constant
In computations, the composition of products at
reactor performance allows decreasing Т w , which
the outlet of the tubular reactor is equal to
significantly reduces wall replacement costs, or to boost
experimental values.
reactor performance by decreasing residual methane and/or increasing natural gas consumption.
X CH 4 ( at l = 2.4 m) → max , ∆P → min
Oxidation of formaldehyde into formic acid [26]
u = 0.935–1.55 m/s, Dtube = 14–22 mm,
Resident time can be reduced more than twice at the same Experimental data were obtained at a pilot plant.
Т cool = 79–115 °C, Relative catalyst
X CH 2O of 98,5 % using two cascade reactors which differ
Relative errors of temperatures along the reactor
either by catalyst activity or Т cool . Based on modeling
axis, total X CH 2O and total SCH 2O2 were 2,8; 1,43
results, a pilot plant was designed and constructed.
and 1,27 %, especially in the case of using two
in activity, 0.25–1, CCH = 5.2–6.1 vol.%. 2O
reactors with different Т cool . Catalytic dehydration of ethanol to ethylene [22,27] 16
Dtube = 0.03 m, Ltube = 1.6–3.5 m, Pin = 1–
With an increase in Т w from 400 to 460 oC the
Experimental data were obtained at a pilot plant
3 atm, Т w = 400–460 oC, Residence time
residence time and, consequently, Ltube , required for
with a single-tube reactor. Relative errors of
= 2.1–3.9 s, u = 0.6–1.0 m/s, CCin2 H 5OH =
X C2 H 5OH = 98% are decreased. However, if Т w > 420
90.2, C
in H 2O
= 9.8 (vol.%).
conversion and ethylene selectivity were 0.4; 0.4
o
C YC 2 H 4 decreases because of the increase in
Catalyst pellets: Cylinder Lcat = 5 mm,
butylenes and other by-products formation. The
Dcat = 4 mm, Technological restrictions:
highest X C2 H 5OH = 98.7% and SC2 H 4 = 98.2% were
o X C2 H 5OH > 98%, SC2 H 4 > 97%, ∆P ≤ 0.2 observed at Т w = 420 C. The increase in u leads to the
atm, Optimization criterion: N tube → min.
temperatures along the reactor axis, ethanol
decrease in residence time which is explained by intensification of λr .
17
and 0.2 %, respectively.
As an example the process of methane steam reforming presented in Fig. 2 demonstrates three-dimensional dependencies of methane concentration (Fig. 2a) and temperature (Fig. 2b) on the reactor length and radius and also three-dimensional dependencies of methane concentration on reactor length (near the reactor wall) and catalyst radius in each point lengthwise (Fig. 2c). Kinetic equations and constants were taken from [28].
Figure 2. Three-dimensional dependencies of methane concentration (a) and temperature (b) on the reactor length and radius and three-dimensional dependencies of methane concentration on reactor length (near the reactor wall) and catalyst radius in each point lengthwise (c). Т in = in 793 K, Pin = 25 atm, Ltube = 3 m, Dtube = 0.1 m, u = 9.9 m/s, Т w = 1153 K, CCH = 24.0, 4 in in = 1.6, CHin2 = 1.3, CCO = 0.01, CNin2 = 3.49 (vol.%), Catalyst: cylinder Lcat = CHin2O = 69.6, CCO 2
0.019 m, Dcat = 0.014 m, N hole =4, Dhole = 0.003 m, Deq = 0.002 m.
The methane concentration and temperature profiles are shown in Figure 2. The program computes the concentrations of all substances (СН4, H2O, CO, H2, CO2, N2) across the reactor radius and length and also across the pellet radius. In addition to the above-listed processes the multiscale mathematical model and numerical algorithm were used to model a number of other processes in tubular reactors [22,29].
4.2. Multiscale mathematical modeling of monolith catalytic reactor
A monolith catalytic reactor with short contact time for partial oxidation of methane to syngas has been analysed by means of a dynamic one-dimensional two-phase reactor model with account for both transport limitations in the boundary layer of a fluid near the catalyst surface (third scale) and detailed kinetics via elementary reactions (first scale). The reaction mechanism contains 32 elementary steps describing methane oxidation on platinum, 14 gaseous compounds and 13 intermediates [22,30]. A detailed mathematical model is described in [30]. 4.2.1. Algorithm realization
The details of the algorithm realization in the case of monolith reactor can be found online in Inline Supplementary Algorithm Description S4. The programs are written in Intel(R) Visual Fortran Compiler 11.1. Time of computation of one variant is about a few hours on a Core i5-4570 class computer. Such a long time is appropriate for the non-steady state problem under consideration for which it takes a sufficiently long time to determine a solid-state thermal profile. A flowchart of the algorithm representing mainly the order of computational modules in the algorithm is shown in Inline Supplementary Fig. S3. The algorithm realization involves no difficulties. It is notable that during the warm-up process, when the values for intermediates undergo rapid changes, the size of the time step is
19
automatically set small. When the surface reaches a steady state while the monolith solid phase is slowly reaching its steady state the size of the time step is automatically increased.
4.2.2. Modeling results of monolith reactor
The program provides for varying the shape of the channels thereby modifying the dependences for calculating the coefficients of heat and mass exchange between the gaseous phase in the monolith channels and catalytic surface as well as the geometric sizes of the monolith channels and technological parameters of the process. The effects of such parameters as linear velocity (0.3 – 0.5 Nm/s), equivalent channel diameter (0.38 – 0.67 mm), inlet temperature (400 – 420 oC), inlet methane concentration (0.24 – 0.27 mole fraction in air) and effective thermal conductivity of the monolith (nearly 10 Jm-1s-1K-1) on the process dynamics, specifically on the ignition dynamics have been examined. It was shown that increasing linear velocity and equivalent channel diameter as well as decreasing axial conductivity of solid phase favors decreasing a time delay in syngas production in the Pt/Ce-Zr-La/α-Al2O3 honeycomb monolith with a triangular shape channels [22,30]. The examples of three-dimensional dependencies on time are presented in Fig. 3 for the fraction of surface hydrogen H* (Fig. 3a), concentration of hydrogen in gaseous phase of monolith channel (Fig. 3b) as well as temperature in gaseous phase (Fig. 3c) along the channel length.
20
Figure 3. Three-dimensional dependencies on time for the fraction of intermediate hydrogen H* (a), concentration of hydrogen in gaseous phase of monolith channel (b) and temperature in gaseous phase (c) along the channel length. Monolith with triangular channel cross-section. in = 25 mol.% in air, Т in = 400 oC, Pin = 1 atm, u = 0.3 m/s, deq = 0.67 mm, ε = 0.59. CCH 4
The concentration profiles for one of the thirteen intermediates and one of fourteen gaseous compounds and temperature profile in the monolith channel are presented in Figure 3. The program also performs computations for the remaining intermediates and gaseous compounds both in the monolith channel and on the monolith catalytic surface (in subsurface layer). The monolith solid phase temperature is also calculated. 21
4.2.3. Comparison of modeling results with experimental data
Simulation results are compared with transient experimental data on the reaction ignition in a full-size monolith 0.05 m in length. The time-dependent temperature profiles along the monolith channel and concentrations of substances at the monolith outlet show reasonable agreement between calculated and experimental values. However for the first 20 s the predicted temperature in the front part of the monolith increases slightly faster than the experimentally measured temperature. Both in the experiments and simulations during light-off, carbon dioxide, as a product of complete oxidation of methane, appears initially. Then, syngas selectivity slowly increases with rising temperature [22,30].
4.3. Multiscale mathematical modeling of fluidized bed reactor
Mathematical modelling of a fluidized bed reactor for propane-isobutane dehydrogenation has been performed on the base of a two-phase non-isothermal reactor model taking into account kinetics, convective transfer and axial gas-heat dispersion (third scale). Unique feature of the reactor is feeding of the regenerated hot catalyst into the top of the reactor, thereby increasing the reactor top temperature. The mathematical model is described in [31]. 4.3.1. Algorithm realization
The mass transfer equations are stationary one-dimensional second order differential equations in bubble and dense phases. Assuming the equal temperature in bubble and dense phases the heat transfer equation is a single-phase stationary one-dimensional second order differential equation. To apply the algorithm of Paragraph 3 let us make use of the method of transition to the non-steady state problem. Using the integro-interpolation method and the method of straight lines we get three systems of ODEs for concentrations in dense and bubble phases and for the temperature.
22
At each integration step the systems of ODEs are sequentially solved by method (7)-(10) using tridiagonal matrices and the stepsize adjustment algorithm. The programs are written in Intel(R) Visual Fortran Compiler 11.1. Time of computation of one variant is less than 30 seconds on a Core i5-4570 class computer. 4.3.2. Modeling results of fluidized bed reactor
Process simulation parameters are given in Table 2. An increase in the rate of heated catalyst introduction after regeneration resulted in a growth of average bed temperature ( Т av ). The effect of Т av and inlet propane mass fraction on the isobutylene and propylene yield is shown in Figure 4. Isobutylene yield was calculated as weight fraction in isobutane–isobutylene fraction.
Table 2. The parameters used in the simulation of industrial fluidized bed reactor Parameter Value Reactor diameter, m
4.5
Reactor height, m
4.9
Inlet gas velocity, Nm/s
0.15
Propane fraction in inlet mixture, %
0 – 60
Velocity of downward catalyst moving, m/s
0.01 – 0.1
Inlet catalyst temperature, oC
630
Fraction of bubble phase cross section
0.35
The isobutylene yield goes up along with increasing propane mass fraction and average layer temperature (Fig. 4). The rise of the inlet propane concentration leads to reduction
of
hydrogen
concentration
in
the
reaction
mixture
shifting
isobutane
dehydrogenation equilibrium to the right. The dependences of propylene yield (Fig. 4.) have more complicated behavior.
23
Isobutylene and propylene yield, % weight
70 60 50 40 30 20 0
15
30
45
60
Inlet propane fraction, % weight Figure 4. The effect of inlet propane mass fraction and Т av on the isobutylene (solid lines) and propylene (dashed lines) yield. Black lines (square symbols) correspond to Т av = 560 oC, red lines (circle symbols) correspond to Т av = 585 oC.
4.3.3. Comparison of modeling results with experimental data
The model and numerical algorithm are validated against experimental data obtained for commercial fluid bed isobutane dehydrogenation reactor [22,31]. From Table 3, it can be seen that the model is able to reproduce satisfactorily the dependence of isobutylene yield and selectivity on the average bed temperature.
Table 3. An influence of the average bed temperature on isobutylene selectivity and yield. Modelling results and experimental data. u = 0.15 m/s, Ciin−C4 H10 = 95.5 wt%. Average bed
Isobutylene selectivity, % weight
temperature ,
Experimental
Numerical
Experimental
Numerical
data
results
data
results
o
C
Isobutylene yield, % weight
585
83.2
83.5
38.5
38.0
575
84.2
85.0
34.0
35.0
24
554
85.2
86.0
31.0
30.0
4.4. A comparison of modeling approaches
To compare computation time for the presented approach and conventional approaches let us consider, for instance, the application of the straight line method in their realization. Then in case of two-scale process modeling with model dimensionality ≤ 2D on each scale, for instance, on reactor and catalyst pellet scales, the presented approach ensures a significant decrease in system dimensionality, thereby a decrease in computation time (Table 4).
Table 4. Comparison of computation time for different approaches. Approaches
Сonstraint
Conventional Second-order derivatives approach.
Order of solution
Computation time
Simultaneous solution of
f(N·NL·NP)
with respect to all coordinates a system of equations on on each scale are normal.
all scales.
Presented
First-order derivative with
Separate solution of a
approach.
respect to a coordinate on
system of equations on
each scale.
each scale.
f(NP) ·N+ f(N)
Besides, to solve the resultant systems of linear equations of type (9) traditionally either the LU decomposition or Gaussian elimination method is used, otherwise, if the system is represented in form (5), the tridiagonal matrix algorithm is used concurrently with the iteration process. In these cases the computation time is proportional either to N3eqs or 8Neqs·Niter, respectively. In our case, following the representation of the system in form (5) the tridiagonal matrix algorithm is used and there is no need to proceed with iterations if one uses the stepsize adjustment algorithm. The computation time is then proportional to 8Neqs.
25
5. Conclusions
EO approach-based multiscale models of heterogeneous catalytic reactors, which are noteworthy as the degree of complexity of such models is dictated only by the necessity, which is sufficient too, may have dimensionality ≤ 2D with the first-order derivative with respect to a coordinate on each scale for a large class of reactors, including nonstandard ones. Simultaneous solution of a system of equations, including all the scales under consideration, which is characterized by high dimensionality, is extremely time-consuming because bulky data storage is required. A strategy of modeling such systems has been suggested based on the use of an identical set of numerical methods, independently on each scale, with step size control implemented with account for the rate of change of the variables on each scale. The introduced algorithm of numerical solution significantly saves computation time and provides adequate accuracy and stability. The developed strategy and the combination of numerical methods suggested for the solution of equations for the model on each scale proved effective for multiscale modeling of heterogeneous catalytic reactors such as tubular reactors, monolith catalytic reactor, and fluidized bed reactor. When switching between the reactor types the parameters of computation algorithms do not need to be corrected. The developed strategy can also be applied to multiscale modeling of heterogeneous catalytic reactors of other types, including nonstandard ones. The models used for modeling of such reactors agree with the modern trends in chemical reaction engineering, such as [8,9,32]: an estimation of transfer coefficients with correlations capturing also the effect on reactor scale; simultaneous modeling of transport-kinetic interactions on a catalytic particle and reactor scales; validation of the models against more than one output parameters and the data obtained under representative current process operating conditions.
26
Acknowledgments.
The author highly acknowledges her colleagues Ludmilla Bobrova, Victor Chumachenko, Nataliya Chumakova, Aigana Kagyrmanova, Vitalii Kashkin, Elena Ovchinnikova, Ilya Zolotarskii (Boreskov Institute of Catalysis SB RAS) for collaboration as well as for fruitful discussions of the heterogeneous catalytic reactors modeling strategy and results presented in this paper. This work was conducted within the framework of budget project No. 0303-2016-0017 for Boreskov Institute of Catalysis.
Nomenclature ai , pi , γij , βij – parameters (Eq. 6);
A, B, C , f (u ) – coefficients in equations (Fig. 1); Ai – tridiagonal matrix depend on the solution in a non-linear manner (Eq. 4); A – tridiagonal matrix of dimensionality NxN (Eq. 5); Ci – vector depend on the solution in a non-linear manner (Eq. 4) or concentration of i-th
compound in gas phase, (mol.%, vol.%, wt.%); C – vector of dimensionality N (Eq. 5); Dn – matrix (Eqs 8-10); D – diameter, (m); Deq – equivalent catalyst pellet diameter as a plate, (m);
deq – equivalent channel diameter, (mm); fu – Jacobi matrix; f – function in Table 4;
27
I – unity matrix; h – stepsize; l – the axial coordinate of the reactor, (m); L – length of the reactor (m);
M – number of variables in Eq. 1; N – number of grid points along the reactor radius; NL – number of grid points along the reactor length; NP – number of grid points along the pellet radius; N text – number (see Subscripts);
P – pressure, (atm); ∆P – pressure drop along the reactor, (atm); R – radius of reactor, (m); S – selectivity, (%); Т – gas temperature, (oC, K); t – time, (s);
u – linear velocity under normal temperature, (m s-1);
ui – variable in Eq. 1; X – conversion, (%);
x – the radial coordinate of the reactor, (m). Y – yield of product, (%); y – the axial coordinate of the reactor, (m).
Greek letters:
α – the heat-transfer coefficient from the tube wall to the bed, (J m-2s-1K-1); ε – porosity or accuracy of the algorithm in Eq. 11;
η – effectiveness factor of catalyst pellet; 28
λr – radial heat conductivity in the catalyst bed, (J m-1s-1K-1); θ – intermediate species, (undimensional); Subscripts: 0
– initial value;
av
– average;
cat
– catalyst;
cool
– cooling agent;
eqs
– equations;
hole
– catalyst hole;
– variable;
i
in
– inlet; – iteration;
iter
j
– point along reactor radius;
tube
w
– tubular reactor;
– tube wall;
References
[1] M. G. Slin’ko, History of the development of mathematical modeling of catalytic processes and reactors, Theoretic. Found. Chem. Eng. 41 (1) (2007) 13–29. [2] J. J. Lerou, K. M. Ng, Chemical reaction engineering: A multiscale approach to a multiobjective task, Chem. Eng. Sci. 51 (10) (1996) 1595–1614. [3] M. G. Slin’ko, Modelirovanie khimicheskikh reaktorov (Modeling of Chemical Reactors), Nauka: Novosibirsk, 1968 (in Russian). [4] W. Marquardt, Trends in computer-aided process modeling, Computers Chem. Engng 20 (1996) 591-609.
29
[5] R. Bogusch, B. Lohmann, W. Marquardt, Computer-aided process modeling with MODKIT, Computers Chem. Engng 25 (2001) 963-995. [6] A. S. Moharir, S. S. Shah, R. D. Gudi, et al, Generalized reactor model: an object oriented approach to reactor modeling, Europ.Symp. on Comp. Aided Process Engng - 11 (2001) 243248. [7] R. Aris, The optimal design of chemical reactors; a study in dynamic programming, New York, Academic Press, 1961. [8] H. Stitt, M. Marigo, S. Wilkinson and T. Dixon, How good is your model?, Johnson Matthey Technol. Rev. 59 (2) (2015) 74-89. [9] M. P. Dudukovic, Reaction engineering: Status and future challenges, Chem. Eng. Sci. 65 (2010) 3-11. [10] F. Zhao, X. Chen, A hybrid numerical-symbolic solving strategy for equation-oriented process simulation and optimization, AIChE J. DOI 10.1002/aic.15622. [11] W. P. M. van Swaaij, A. G. J. van der Ham, A. E. Kronberg, Evolution patterns and family relations in G–S reactors, Chem. Eng. J. 90 (2002) 25–45. [12] L. H. Roblee, R. M. Baird, J. W. Tierney, Radial porosity variations in packed beds, AIChE J. 4 (1958) 460-464. [13] I. I. Ioffe, L. M. Pismen, Engineering Chemistry of Heterogeneous Catalysis, Chimiya, Leningrad, 2nd ed., 1972 (in Russian). [14] A. A. Samarskii, A. P. Mikhailov, Matematicheskoe modelirovanie. Idei. Metody. Primery (Mathematical Modeling: Ideas, Methods, and Examples), Moscow: Fizmatlit, 2001. [15] C. A. J. Fletcher, Computational techniques for fluid dynamics, vol. 1, Springer-Verlag, 1991. [16] H. H. Rosenbrock, General implicit processes for the numerical solution of differential equations, Computer J. 5 (1965) 329-330.
30
[17] E.A. Novikov, Numerical methods for solution of differential equations in chemical kinetics, in: V. I. Bykov (Ed.), Mathematical methods in chemical kinetics, Nauka, Novosibirsk, 1990, 53-68 (in Russian). [18] V. N. Bibin, In Proceedings of XII International Conference on Chemical Reactors, Boreskov Institute of Catalysis, Novosibirsk, 1 (1996) 127–132 (in Russian). [19] N. V. Vernikovskaya, T. L. Pavlova, N. A. Chumakova, A. S. Noskov, Mathematical modelling of filtration and catalytic oxidation of diesel particulates in filter porous media, in: J.C. Cortes, L. Jodar, R. Villanueva (Eds.), Mathematical Modeling in Social Sciences and engineering, Nova Sciences Publishers Inc., Hauppauge, NY, 2014 27–40. [20] L. A. Rapatskii, L. V. Yausheva, Mathematical Modeling of Catalytic Processes in Porous Media. Preprint of Computation Center, Siberian Div., Russ. Acad. Sci.: Novosibirsk 1019 1994 (in Russian). [21] A. P. Kagyrmanova, I.A. Zolotarskii, E. I. Smirnov, et al, Optimum dimensions of shaped steam reforming catalysts, Chem. Eng. J. 134 (2007) 228-234. [22] L. N. Bobrova, N. V. Vernikovskaya, A. P. Kagyrmanova, V. N. Kashkin, I. A. Zolotarskii, Application of Mathematical Modeling and Simulation to Study the Complexity of the Basic Chemical and Physical Phenomena in Catalytic Reactors, in Book: Mathematical Modelling. Editors: Christopher R. Brennan. Nova Science Publishers 2012, 475-577. [23] J. R. Rostrup-Nielsen, Catalytic Steam Reforming, Berlin: Springer, 1984. [24] J. R. Rostrup-Nielsen, L. J. Christiansen, J. H. B. Hansen, Activity of Steam Reforming Catalysts: Role and Assessment, Applied Catalysis 43 (1988) 287-303. [25] A. P. Kagyrmanova, Optimization of catalyst shape and size in tubular packed bed reactor, PhD Thesis, Novosibirsk, 2009.
31
[26] I. A. Zolotarskii, T. V. Andrushkevich, G. Ya. Popova, et al, Modeling, design and operation of pilot plant for two-stage oxidation of methanol into formic acid, Chem. Eng. J. 238 (2014) 111-119. [27] A. P.
Kagyrmanova,
V. A.
Chumachenko,
V. N.
Korotkikh, et al, Catalytic
dehydration of bioethanol to ethylene: Pilot-scale studies and process simulation, Chem. Eng. J. 176 –177 (2011) 188 –194. [28] J. Xu, G. F. Froment, Methane steam reforming: II. Diffusional limitations and reactor simulation, AIChE J. 35 (1989) 97. [29] E. V. Ovchinnikova, N. V. Vernikovskaya, T. V. Andrushkevich, V. A. Chumachenko, Mathematical modeling of β-picoline oxidation to nicotinic acid in multitubular reactor: Effect of the gas recycle, Chem. Eng. J. 176-177 (2011) 114-123. [30] N. V. Vernikovskaya, L. N. Bobrova, L. G. Pinaeva, et al, Transient behavior of the methane partial oxidation in a short contact time reactor: Modeling on the base of catalyst detailed chemistry, Chem. Eng. J. 134 (2007) 180-189. [31] N. V. Vernikovskaya, I. G. Savin, V. N. Kashkin, et al, Dehydrogenation of propane– isobutane mixture in a fluidized bed reactor over Cr2O3/Al2O3 catalyst: Experimental studies and mathematical modelling, Chem. Eng. J. 176-177 (2011) 158-164. [32] M. P. Dudukovic, Trends in catalytic reaction engineering, Catal. Today 48 (1999) 5-15.
Captions:
Figure 1. Schematic representation of the structure of partial differential equations: (a) – with respect to space coordinate “y” contains only first-order derivatives, (b) – the first-order derivative with respect to time.
32
Figure 2. Three-dimensional dependencies of methane concentration (a) and temperature (b) on the reactor length and radius and three-dimensional dependencies of methane concentration on reactor length (near the reactor wall) and catalyst radius in each point lengthwise (c). Т in = in 793 K, Pin = 25 atm, Ltube = 3 m, Dtube = 0.1 m, u = 9.9 m/s, Т w = 1153 K, CCH = 24.0, 4
in in = 1.6, CHin2 = 1.3, CCO = 0.01, CNin2 = 3.49 (vol.%), Catalyst: cylinder Lcat = CHin2O = 69.6, CCO 2
0.019 m, Dcat = 0.014 m, N hole =4, Dhole = 0.003 m, Deq = 0.002 m.
Figure 3. Three-dimensional dependencies on time for the fraction of intermediate hydrogen H* (a), concentration of hydrogen in gaseous phase of monolith channel (b) and temperature in gaseous phase (c) along the channel length. Monolith with triangular channel cross-section. in CCH = 25 mol.% in air, Т in = 400 oC, Pin = 1 atm, u = 0.3 m/s, deq = 0.67 mm, ε = 0.59. 4
Figure 4. The effect of inlet propane mass fraction and Т av on the isobutylene (solid lines) and propylene (dashed lines) yield. Black lines (square symbols) correspond to Т av = 560 oC, red lines (circle symbols) correspond to Т av = 585 oC.
33
A novel strategy is proposed for multiscale modeling of catalytic reactors. This strategy was applied for modeling tubular, monolith, fluidized bed reactors. Robustness of the algorithm is shown for multiscale modeling of these reactors.
34