- Email: [email protected]
Large-scale DAE-constrained optimization applied to a modified spouted bed reactor for ethylene production from methane
Accepted Manuscript
Large-Scale DAE-Constrained Optimization Applied to a Modified Spouted Bed Reactor for Ethylene Production from Methane D.M. Yancy-Caballero, L.T. Biegler, R. Guirardello PII: DOI: Reference:
S0098-1354(18)30194-7 10.1016/j.compchemeng.2018.03.017 CACE 6056
To appear in:
Computers and Chemical Engineering
Received date: Revised date: Accepted date:
13 October 2017 16 March 2018 18 March 2018
Please cite this article as: D.M. Yancy-Caballero, L.T. Biegler, R. Guirardello, Large-Scale DAEConstrained Optimization Applied to a Modified Spouted Bed Reactor for Ethylene Production from Methane, Computers and Chemical Engineering (2018), doi: 10.1016/j.compchemeng.2018.03.017
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.
ACCEPTED MANUSCRIPT
Highlights • A spouted bed reactor was proposed to enhance the performance of the oxidative coupling of methane.
CR IP T
• An optimization strategy was developed to demonstrate the theoretical feasibility of the proposed reactor.
• The use a porous membrane as a draft tube in a spouted bed improves the heat and mass
AC
CE
PT
ED
M
AN US
transfer.
1
ACCEPTED MANUSCRIPT
Large-Scale DAE-Constrained Optimization Applied to a Modified Spouted Bed Reactor for Ethylene Production from Methane D. M. Yancy-Caballeroa , L. T. Bieglerc , R. Guirardellob,∗ a Brazilian
CR IP T
Bioethanol Science and Technology Laboratory (CTBE), Brazilian Center for Research in Energy and Materials (CNPEM), R. Giuseppe M´ aximo Scolfaro, 10000, 13083-970, Campinas-SP, Brazil. b School of Chemical Engineering, University of Campinas, Av. Albert Einstein 500, 13083-852, Campinas-SP, Brazil c Department of Chemical Engineering, Carnegie Mellon University, Doherty Hall, 5000 Forbes Avenue Pittsburgh-PA, United States
AN US
Abstract
In this paper, a modified spouted bed reactor is proposed to enhance the yield of the oxidative coupling of methane (OCM). Optimization techniques are used to carry out a theoretical analysis of ethylene production via OCM and define some optimal operating conditions of the reacting system. A model-based DAE-constrained optimization strategy is proposed and applied to the OCM process
M
to illustrate the computational capability of the proposed formulation, and the theoretical feasibility of the proposed reactor. The model developed for the reactor is a one-dimensional model composed
ED
of material, energy, and momentum balances. This model along with the kinetic model constitute a non-linear and differential-algebraic system, which is discretized using orthogonal collocation on
PT
finite elements with continuous profiles approximated by Lagrange polynomials. The resulting
CE
algebraic collocation equations are written as equality constraints in the optimization problem, which is solved with the IPOPT solver within the optimization-modeling platform. An initialization
AC
routine based on simulations was carried out to guarantee convergence in optimizations. Results from simulations and optimizations showed the potential of combining different reactor concepts to ∗ Corresponding
author: Tel: +55 19 3521 3955; fax: +55 19 3521 3965 Email address: [email protected]; P.O. Box 6066. (R. Guirardello)
Preprint submitted to Chemical Engineering Journal
March 19, 2018
ACCEPTED MANUSCRIPT
improve the ethylene production from natural gas via oxidative coupling of methane. Keywords: Dynamic optimization, oxidative coupling of methane, spouted bed reactor, modeling, simulation, ethylene, methane
CR IP T
1. Introduction Natural gas is the third most important energy source in the world nowadays, but only a portion of the extracted natural gas is transported for use since most of it is burned or re-injected into the reservoir [1, 2]. The use of natural gas can be improved if local conversion techniques are developed,
AN US
considering that the natural gas price depends mainly on the transportation. Therefore, conversion processes of methane, which is the principal component of natural gas, are encouraged for natural gas valorization. Several studies have been published on the efficient use of natural gas in chemical and petrochemical industry, taking advantage in this manner of the enormous potential of the world
M
reserves of this natural source of hydrocarbons. Efforts have focused especially in process design problems [3, 4].
ED
The amount of methane used as raw material in the chemical industry represents between 57% of the total consumption [5]. In the long term, with the rapid depletion of crude oil and the
PT
environmental impacts of the oil-derived compounds, the natural gas can become the primary source of energy in the 21st century.
CE
Direct and indirect conversion methods are well known for the conversion of the methane to high value-added chemical products. The indirect use of methane is a more attractive way compared to
AC
the direct one since higher yields are achieved, being these the preferred processes for commercial applications of producing fuels in large-scale [6]. However, these methods are complex processes that require an intermediate step for producing syngas, which consumes enormous amounts of energy causing high production costs [7]. In this sense, direct conversion methods are promising because 3
ACCEPTED MANUSCRIPT
the methane is directly converted without requiring the syngas production step. Partial oxidation and oxidative coupling of methane (OCM) are the principal routes of direct conversion, where methanol, formaldehyde or ethylene can be produced from methane in just one step, but these
CR IP T
processes are challenging because the produced compounds are more active than methane, causing reduced yields [8, 9]. Among the direct conversion methods, the oxidative coupling of methane has received particular attention due to a higher yield than other direct methods. Furthermore, this can be carried out at lower temperatures than the conventional processes for ethylene production [10], and because of it is an exothermic process, the produced energy may be used for purposes of
AN US
steam production.
Two obstacles have hampered the industrial application of the OCM process: the reduced yield of the reaction and the large amounts of energy released [11, 12, 13]. From a practical standpoint, the first limitation can be addressed with the development of more selective catalysts and the
M
advances in the chemical reaction engineering. However, the second limitation, relate to the large heat released, can only be resolved by improving and developing new designs of chemical reactors
ED
since the thermodynamics and fluid dynamics also influence this restriction, and not only the reaction kinetics [14, 15].
PT
Several reactor types and configurations have been proposed for the OCM process, some of the most studied are fixed bed reactors, fluidized bed reactors, membrane reactors, countercurrent
CE
moving bed chromatographic reactors, and electrochemical reactors. All of these types of reactors have its advantages and disadvantages, but only the first three ones have been explored and
AC
developed in detail [14]. Experimental and simulation works have been reported in the literature trying to understand the phenomenological behavior of these reactor types [16, 17, 14, 18, 19, 15]. A distributed oxygen feed along the reactor is just one of the suggested proposals to increase the
4
ACCEPTED MANUSCRIPT
yield of the OCM [20]. In this manner, the application of a membrane reactor to carry out the OCM was proposed to allow flows of small amounts of oxygen into the reaction zone to activate the methane and enable
CR IP T
the coupling of methyl radicals that are formed, limiting at the same time undesired complete combustion reactions. Different membrane types have been used for this purpose (porous membranes, dense membranes, catalytic membranes and perovskite-type membranes). Theoretical studies have reported yields higher than 30% in conventional membrane reactors [21, 22], while experimental ones have reported yields of 28% in porous membrane reactors [23] and around 35% in catalytic
AN US
membrane reactors [24]. Another innovative type of reactor proposed for the OCM is the countercurrent moving bed chromatographic reactor. The highest yields achieved in the oxidative coupling of methane have been reported in this kind of reactor (>50 %) [16, 25]. However, this type of reactor still has to be developed in detail before any large-scale application.
M
Fluidized bed reactors have been considered by some researchers as the best type of reactors for the OCM due to the huge heat released during the reaction, which is more complicated to be
ED
controlled in other types of reactors like the fixed bed reactors [11, 14]. Furthermore, the ability of these reactors to operate isothermally, allow for continuous recirculation and the ease of replacing
PT
the deactivated catalyst, become the fluidized bed reactor a potentially attractive reactor design for the OCM, but most of the selective catalysts have shown agglomerations in fluidized bed reactor
CE
causing problems to get a stable fluidization. Most of the studies about OCM in fluidized bed reactors reported yields less than 20%, which limits the applicability of this type of reactor.
AC
Recently, the Wozny’s research group at Berlin Institute of Technology proposed a reactor
network based on membrane and plug flow/fixed-bed reactors for improving the performance of the OCM process. They investigated feeding policies for an optimal operation of the process combining
5
ACCEPTED MANUSCRIPT
several types of reactors [26, 27, 28, 29, 30]. Performance of 23% yield, 54% selectivity, and 43% methane conversion was reported. In this context, the current work proposes a modified spouted bed reactor as a new design
CR IP T
concept for improving the performance of the oxidative coupling of methane. While the use of the spouted bed technique as a chemical reactor is not novelty since this technique has been used in processes involving gas-solid reactions such as pyrolysis of polyolefins, combustion, and gasification of coal and biomass [31, 32, 33, 34, 35, 36], some significant modifications to the conventional configuration are presented in this work aiming to maximize the ethylene yield. Another valuable
AN US
contribution of the present paper is a comprehensive model-based DAE-constrained optimization strategy, which was developed and implemented to find optimal operating and design conditions of the proposed spouted bed reactor for ethylene production from methane. The results achieved in this study demonstrates the theoretical feasibility of employing the spouted bed technique to
M
enhance the performance of the oxidative coupling of methane.
ED
2. Methodology 2.1. Reactor model description
PT
The spouted bed reactor model was developed by applying mass, momentum and energy balances in one-dimension under steady-state. The proposed reactor was divided into three hydro-
CE
dynamic regions [37, 38, 39, 40]: the spout region, where gas and particle velocities are high, and the particle concentration is low. In this region, the particles and the fluid are in co-current flow;
AC
the annular region between the spout region and the column wall, which is characterized by very low velocities, and where the particle concentration is close to that found in a packed bed. In this region, the particles and the fluid are in counter-current flow, and also there is a percolation of the
6
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 1: Different regions of the spouted bed reactor
fluid from the spout into the annulus through the spout-annulus interface; and a fountain region with a much lower particle concentration than in a packed bed.
The spout and the annular region were taken into account in the conical and cylindrical part
M
of the reactor separately since a cross flow of both the fluid and particles between the spout and the annulus exists only in the conical part for the spouted bed with draft tube. In the cylindrical
ED
part, there is not such a cross flow because of the insertion of the internal draft tube, unless that a porous draft tube is being modeled. For the conventional spouted bed, this cross flow exists along
PT
the whole reactor in both the conical and cylindrical part. This division can be seen in Fig. 1 for the spouted bed with draft tube. It can be seen that both the conical part (region I) and the cylindrical
CE
part (region II) were divided into spout and annulus. The fountain region was not considered in the modeling because this region has no significant effects on the reactor performance results for very
AC
fast catalytic reactions due to the low particle concentration in this region [41, 42]. The division for the conventional spouted bed was the same, with the only difference that the internal draft tube was not taken into account in the model, allowing for a cross flow of the fluid and particles in the
7
(b)
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
(c)
(d)
Figure 2: Control volumes for balances. Global mass balance in spout (a) and annulus (b). Momentum balance of gas in spout (c) and annulus (d)
M
cylindrical part as well.
ED
The spouted bed reactor model for all cases investigated in this work was considered adiabatic, and a pseudo-homogeneous approach was adopted for the energy balance, which simplifies the model
PT
considerably. This means that the temperatures for the gas and solid phase were considered to be approximately the same. While a adiabatic reactor was assumed, it should be highlighted that
CE
in some real cases such an assumption might not hold true since the large heat of reaction of the oxidative coupling of methane could cause the increase of the temperature forming sharp peaks,
AC
which must be controlled by using cooling designs. Another assumption in the development of the reactor model was to consider as negligible the
fluid-particle mass transfer resistance, so that component concentrations on the catalyst surface and in the bulk fluid were assumed to be equal. Fig. (2) shows some selected control volumes for 8
ACCEPTED MANUSCRIPT
which the governing equations (Table 1) have been derived. The kinetic model used in this work describes the oxidative coupling of methane with ten reactions including nine chemical species. This reaction model has been reported in detail by Stansh
CR IP T
et al., 1997 [43], and is based on kinetic measurements in a fixed bed micro-catalytic reactor using La2 O3 /CaO as catalyst. The model was obtained over a wide range of reaction conditions, which increases its applicability in simulations (1 < PO2 < 20 kPa, 10 < PCH4 < 95 kPa, 973 < T < 1228 ˙ STC < 250 kg · s · m−3 ). Further details on the kinetic model and the used parameK, 0.76 < mcat /Q ters as well as the thermodynamic data for all compounds included in simulations and optimizations
AC
CE
PT
ED
M
AN US
can be found in Appendix B and Appendix C.
9
ACCEPTED MANUSCRIPT
Conical
Conservation equation Spout
Equation
Global mass balance
Gas phase
1 As
Solid phase
1 As
·
dAs dz dAs dz
Gas phase
us ·
Solid phase
vs ·
dus dz dvs dz
1 us
+
·
dus dz dvs dz
1 vs · dwi,s dz =
+
Us · ρg,s ·
Component mass balance Momentum balance
·
+ −
1 εs
·
dεs dz
1 (1−εs )
·
dP dz
βB ·(us −vs ) ρp ·(1−εs )
−
+g·
g εs
dεs dz
−
ρg,s ρp
·
=
βB ·(us −vs ) ρg,s ·εs
−1
dTs dz
[Us · ρg,s · Cpg,s + Vs · ρp · Cpp,s ] · +
Annulus
Global mass balance
Gas phase
Solid phase Component mass balance Momentum balance
Gas phase
Solid phase
10 Cylindrical
Spout
Global mass balance
M
Energy balance
Gas phase
ED
Solid phase
Component mass balance
PT
Energy balance Annulus
Global mass balance
CE
Component mass balance
AC
Energy balance
P = rr ˙ j,s · (−4Hrj,s )
j ρp ·Cpp,s d[As ·Vs ] · · [T − T ] a s As dz dρg,a dAa dua 4Ur ·ds 1 1 1 Aa · dz + ua · dz + ρg,a · dz = ua ·εa ·(d2c −d2s ) dAa dva 4Vr ·ds 1 1 Aa · dz + va · dz = va ·(1−εa )·(d2c −d2s ) d[Aa ·Ua ·ρg,a ] dw Ua · ρg,a · dzi,a = cr ˙ i,a · MWi + A1a · dz βB ·(ua +va ) g dua 1 dP ua · dz = − ρg,a ·εa · dz − εa − ρg,a ·εa ρg,a βB ·(ua +va ) a va · dv dz = ρp ·(1−εa ) + g · ρp − 1 "
AN US
Energy balance
dρg,s 4·Ur dz = − εs ·us ·ds 4·Vr (1−εs )·vs ·ds
1 ρg,s
+
cr ˙ i,s · MWi
= − ρg,s1·εs ·
=
CR IP T
Table 1: Mass, momentum and energy balance equations in the different regions of the spouted bed
Region
(1) (2) (3) (4) (5) (6) (7) (8)
· [wi,s − wi,a ]
(9) (10) (11)
ncr P
a [Va · ρp · Cpp,a + Ua · ρg,a · Cpg,a ] · dT rr ˙ j,a · (−4Hrj,a ) dz = ρp · (1 − εa )· j # ngr P d[Aa ·Ua ·ρg,a ] Cpg,a · Aa · (Ts − Ta ) + εa · rr ˙ j,a · (−4Hrj,a ) + dz
(12)
j
1 us 1 vs
·
·
Us ·
dρg,s 1 s · dε dz + ρg,s · dz dεs 1 − (1−εs ) · dz = 0 dw ρg,s · dzi,s = cr ˙ i,s · MWi
dus dz dvs dz
+
1 εs
=0
(15)
[Us · ρg,s · Cpg,s + Vs · ρp · Cpp,s ] ·
dTs dz
00
+Udt · Adt · [Tp,a − Tg,s ]
Gas phase
1 Aa
Solid phase
1 Aa
·
dAa dz dAa dz
+
1 ua
·
+ v1a · · dw Ua · ρg,a · dzi,a
dua dz dva dz
+
1 ρg,a
·
dρg,a dz
P = rr ˙ j,s · (−4Hrj,s ) j
=0
(16) (17)
=0
(18)
= cr ˙ i,a · MWi
(19)
[Va · ρp · Cpp,a + Ua · ρg,a · Cpg,a ] · 00
(13) (14)
00
dTa dz
−K·
dT2a dz2
+Udt · Adt · [Tp,a − Tg,s ] + Uw · Aw · [Tp,a − Tw ]
=
P rr ˙ j,s · (−4Hrj,s ) j
(20)
ACCEPTED MANUSCRIPT
2.2. Simulation of conventional spouted bed reactor When solving the model for the conventional spouted bed reactor, some boundary conditions have to be established to get a realistic solution of the DAE system. These boundary conditions are
CR IP T
easy to be defined at the inlet of the spout region since the inlet gas concentrations and temperature are well defined at the entrance of the bed, which gives the fluid conditions at the bottom of the spout. However, the boundary conditions at the bottom of the annulus are not that clear as there is no a feed stream of fluid in this region for conventional spouted beds, which forces the inlet gas velocity in the annulus to be zero. In consequence, the concentrations of components at the base
AN US
of the annulus cannot be assumed the same as in the inlet of the spout region. This particular condition of conventional spouted bed reactors leads to formulating complex boundary conditions to solve the model. In this way, a suitable boundary condition to be assumed in the annulus would be to satisfy an overall energy balance between the spout and the annulus at the exit of the reactor,
M
since the only knowledge of the temperature at the base of the annulus is that this depends on the temperature of solids moving downward in the annulus, but these temperatures are unknown. This
ED
boundary condition at the exit of the reactor can be simplified as Ta (H) = Ts (H). The knowledge a priori of the concentration of components at the base of the annulus region
PT
is also crucial for solving the DAE system that represents a conventional spouted bed reactor. A suitable way to estimate these concentrations would be to determine the residence time at z = 0
CE
that a particular component takes in the transition spout-annulus. This residence time can be calculated by assuming that the reactor behaves as a continuous flow stirred-tank reactor (CSTR)
AC
at an infinitesimally small control volume of the reactor at the base of the annulus as suggested by Smith et al., 1982 [44]. Thus, this residence time for a cylindrical spouted bed with conical base can be calculated as τg0 =
Hm/6·Umf
(see [44] for details). The initial concentration of all components
11
ACCEPTED MANUSCRIPT
in the annulus can be estimated by applying a mass balance for each component on the assumed CSTR as follows: ci,a0 = ci,s0 + τg0 · cr ˙ i,a (Ta0 )
CR IP T
(21)
With this set of equations, the temperature at the base of the annulus Ta0 is iterated until some convergence criterion is reached to satisfy the boundary condition at the exit of the reactor. This criterion is given by the difference between both temperatures in the annulus and the spout at z = H, which have to be minimum. Details on the implementation of this algorithm in GAMS
AN US
are presented below, and a schematic representation is shown in Fig. A.1. 1. It is assumed a temperature at the base of the annulus region Ta0
2. The system of non-linear algebraic equations given by equation 21 is solved using the CONOPT solver. This step is defined in a separate GAMS script file, which is then nested within the
M
main program. Such a nesting procedure in GAMS was done by using the “call ” and “execute” GAMS commands in the main program. These commands execute the “GAMS script
ED
file” in runtime and store the results in a “.gdx ” file, which are loaded back into the main program.
PT
3. With the guessed Ta0 , the calculated ci,a0 and, the fixed boundary conditions at the inlet of the spout Ts0 and ci,s0 , the DAE system composed by mass, momentum and energy balances
CE
in the two regions of the reactor is solved as an initial value problem using the OCFE method. The solver used to address the problem at this step is IPOPT
AC
4. A convergence criterion is checked (to avoid excessively large computation times and convergence issues, initially a difference of 5 K between the temperatures in the annulus and the spout was defined as the convergence criterion). If this criterion is not satisfied, then a
12
ACCEPTED MANUSCRIPT
new Ta0 is assumed (employing the step 5) and the algorithm continues with step 2. If this criterion is satisfied, then the algorithm continues with step 6 5. It can be seen that the steps 1-4 represent a variant of the well-known shooting methods.
CR IP T
These methods convert a boundary-value problem into an equivalent initial-value problem. Thus, the unspecified initial conditions of the system of differential equations are guessed, and the equations are integrated forward as a set of simultaneous initial value problems. At the end, the calculated final values are compared with the boundary conditions and the guessed initial conditions are corrected if necessary. This procedure is repeated until the specified
AN US
terminal values are achieved within a small convergence criterion. A trial-and-error approach can be implemented to guess the initial conditions every time. However, this is tedious and impractical because the trajectory of the right solution is unknown leading to a relatively large CPU-time-consuming. In this work, instead of using a trial-and-error approach, the false
M
position method was used to estimate the new Ta0 . This approach was applied successfully to this boundary-value problem using another separate GAMS script file and nesting it within
ED
the main program. This approach showed good numerical stability 6. Since the convergence criterion adopted in step 4 does not guarantee the equality of temper-
PT
atures at the exit, this step is required to refine the solution by solving the whole problem as a boundary value problem. In this step, the last solution found in step 3 for all of the
CE
considered variables is used as an initial guess of the new boundary value problem. This time, the system of non-linear algebraic equations given by equation 21 and the boundary condition
AC
at the exit of the reactor given by Ta (H) = Ts (H) are included in the GAMS formulation as equality constraints, and the problem is solved again by using the OCFE method, and the IPOPT solver
13
ACCEPTED MANUSCRIPT
2.3. Optimization strategy The objective function of the optimization problem considered in this work was the C2+ yield (YC2+ ), which depends on the methane conversion (XCH4 ) and the C2+ selectivity (SC2+ ). C2+
CR IP T
is referred to C2 H6 and C2 H4 together. The equality constraints of this optimization problem are composed by a highly coupled and non-linear set of DAEs distributed in space, which was discretized by using the method of orthogonal collocation on finite elements, along with the GaussRadau collocation points because of their compatibility with the NLP formulations and numerical stability properties [45]. These collocation points satisfy the continuity of the states profiles across
AN US
element boundaries automatically, and no additional constraints have to be included in the model. The optimization problem was formulated within GAMS optimization-modeling platform, and the IPOPT solver, which is based on interior-point methods [46], was chosen to solve the problem. This solver can handle quite efficiently problems with large size, a large number of equality and
M
inequality constraints and a potentially large number of degrees of freedom [47, 45]. Before activating the optimizer to carry out any optimizations, some simulation studies were
ED
conducted to determine a proper number of finite elements to get a good accuracy in a reasonable computing time, and to use the simulation results as initial guesses of all variables in optimizations.
PT
These previous simulations also provided a first insight of the behavior of the solution, so that useful ranges of all variables were established for the subsequent optimizations. Initial guesses for
CE
variables were also needed for simulations due to the number of variables involved. So, a systematic initialization routine was proposed where the problem was solved element to element until reaching
AC
the total number of finite elements. Next section details such a routine. For the simulation cases, the number of equality constraints and variables were the same (square
systems). This was accomplished in GAMS by employing a dummy objective function, which can
14
ACCEPTED MANUSCRIPT
simply be any constant value. In this way, the solver tries to find a feasible solution to satisfy all of the equality constraints. However, for the optimization cases, it is always necessary to have more variables than equality constraints, introducing some degrees of freedom into the model. Inequality
CR IP T
constraints were also included in the optimizations based on upper and lower bounds, not only for decision variables but all of the variables in the model. For simulations, however, no bounds for variables were required since the solver was much faster in this case.
2.3.1. Initialization algorithm used in simulations
AN US
When no initial guesses were available to carry out a simulation, it was necessary to divide the entire domain into a higher number of finite elements and then solve the problem element to element sequentially. Thus, the model was firstly solved only for the first finite element, and as this has a very small length, the initial conditions or boundary conditions at the inlet were used as initial guesses for all variables defined in the interior collocation points.
M
After the model was successfully solved for the first finite element, the solution obtained for the
ED
variables at the end boundary of this finite element was employed as an initial guess for solving the model in the next finite element, and so on in a loop that increases the number of finite elements every time until a maximum number of finite elements was reached. In this manner, the solution
PT
for the whole problem was obtained as a sequence of solutions for small problems defined in every
CE
single finite element. Fig. 3 shows this procedure schematically, and an algorithm is presented in Fig. A.2.
AC
2.3.2. Optimization algorithm of the spouted bed reactor Optimizations were solved by following a sequence of three stages: first, a simulation of the
model was solved by fixing the decision variables to constant values. Second, an optimization
15
ACCEPTED MANUSCRIPT
Finite element, 1
Finite element, 2
Finite element, 3
CR IP T
Ic
Finite element, 4
Figure 3: Graphical representation of the initialization routine used in simulations for a sequence of four finite
AN US
elements. Ic stands the initial conditions
was solved using the solution obtained previously in the simulation as initial guesses for all of the variables. In this stage, the decision variables were freed (to have some degrees of freedom) to very tight bounds (upper and lower bounds) for all of them; these bounds were established by
M
defining an increasing factor 4 of ±5% the constant values used in the previous simulation. Third,
ED
a series of optimizations were solved sequentially while slowly increasing the interval defined by the upper and lower bounds (defining new bounds very close to the old ones). In every optimization, initial guesses for all variables and multiplier bounds were taken from the solution of the previous
PT
optimization. This process finished when the bounds of all of the decision variables were inactive,
CE
i.e., the variables were no longer at their bounds, or these could not be moved anymore (negative values for some variables were not allowed). The obtained solution at this point was assumed to be
AC
the final optimal solution, and calculations end up. This optimization procedure is shown in Fig. A.3.
16
ACCEPTED MANUSCRIPT
3. Results and Discussion 3.1. Validation of the fluid-dynamic model A single simulation of a spouted bed without considering reaction is presented in this section in
CR IP T
order to evaluate the fluid-dynamic model. For such a purpose, the model developed in this work was adapted to represent a conical spouted bed with the aim of comparing the simulation results with those results reported by Olazar et al., 2009 [39]. This author also developed a one-dimensional model similar to the model presented here. For this reason, a conical spouted bed was chosen as a
AN US
good reference point to prove the accuracy and reliability of the proposed fluid-dynamic model.
The main goal of carrying out this simulation is to guarantee the use of a relatively simple but reliable model capable of predicting the velocity distributions of gas and particles in every region of a spouted bed as well as the local voidage variation. All of this information is fundamental to solve the spouted bed reactor model (mass balance for component and energy balances).
M
The simulated spouted bed can be seen in Fig. 4, and the used equations for the model are those presented in Table 1 for momentum and global mass balances, where only the conical part
ED
without the draft tube was taken into account. Table 2 presents all of the parameters used in the simulation. The boundary conditions are given by: at the inlet, the axial gas velocity is up to twice
PT
the minimum spouting velocity referred to inlet cross-sectional area, while axial solids velocity is
CE
zero. The pressure at the inlet is given by equation 22 [39], and the voidage is assumed to be 1. At the outlet, the pressure is atmospheric. The axial component of both gas and particles velocities are
AC
zero at the spout-annulus interface. Instead, a radial component for these velocities was considered at this point. h γ i−0.11 −0.06 · (Re0 )ms · D0−0.08 · H1.08 P0 =Patm + 1.20 · ρb · g · tan 2
(22)
where ρb is the bulk density, γ is the cone angle, D0 is the diameter of the bed inlet, H is the height 17
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 4: Graphical representation of the simulated conical spouted bed. Source: taken from Olazar et al., 2009 [39]
Table 2: Parameters used in the simulation of the conical spouted bed
Symbol
Value
Units
γ
33
◦
Patm
101325
Pa
dp
3
mm
g
9.8
Particle density
ρp
2360
Gas density
ρg
1.18
Gas viscosity
µ
1.86e-5
m s2 kg m3 kg m3 kg m·s
D0
0.04
m
Stagnant bed height
H
0.2
m
Column base diameter
Di
0.06
m
Annulus voidage
εa
0.45
−
Spout diameter
ds
0.04
m
M
Parameter
Angle of the column
ED
Atmospheric pressure Particle diameter
AC
CE
PT
Gravity
Inlet diameter
18
ACCEPTED MANUSCRIPT
of the stagnant bed, and (Re0 )ms is the Reynolds number of minimum spouting referred to D0 , which is given by [39]: 0.5
(Re0 )ms =0.126 · Ar
·
Di D0
1.68 h γ i−0.57 · tan 2
(23)
CR IP T
Ar is the Archimedes number, and Di is the diameter of the bed bottom.
Fig. 5 shows a comparison between the simulated results obtained in this work and those results reported by Olazar et al., 2009 [39], as well as a comparison with experimental data for some profiles such as the voidage profile and, gas and particle velocity profiles in the spout. The experimental
AN US
data were reported in [48, 49], and also used by Olazar et al., 2009 [39] for comparison purposes. The model was implemented in GAMS using the OCFE method. After the discretization, the final NLP is composed of a square system of 1700 variables and 1700 equality constraints. The solution of this NLP problem was achieved in GAMS using a dummy objective function. The simulation
M
time was 52 CPUs (4GB and 2.2GHz HP laptop) using the IPOPT solver. A minimum of 15 finite elements and three collocations points were required for a satisfactory
ED
convergence and accuracy. Problems of numerical stability were detected when a smaller number of finite elements was used to solve the model. As observed from Fig. 5, a simple one-dimensional
PT
model provides adequate predictions for gas velocity in spout and voidage distributions. However, the model was unable to predict the particle velocities accurately in the spout. Similar results were
CE
reported by Olazar et al., 2009 [39]. When the results obtained here are compared with the results reported by Olazar et al., 2009 [39],
AC
it can be observed a slightly quantitative difference, even although the profiles agree qualitatively. This difference is mainly caused by some errors in the Olazar et al., 2009 [39] formulation, more
specifically the reported equations for momentum balances, for further details see [50]. Also, it can
19
ACCEPTED MANUSCRIPT
30 Experimental Olazar et al. (2009) This Work
0.9
Experimental Olazar et al. (2009) This Work
28 26 24
0.8
22
us
εs
20 0.7
18
14 12
0.5
10 0.4
0
0.04
0.08
0.12
0.16
8
0.2
0
0.04
z [m]
CR IP T
16
0.6
0.08
0.12
0.16
0.2
z [m]
(a)
(b)
3.5
3
AN US
3
2.5
vs
ua
2
1.5
1.5
1
Experimental Olazar et al. (2009) This Work 0
0.04
0.08
0.12
z [m]
0.06
PT
va
0.08
0.04
CE
0.02
0.04
0.08
0
0.04
0.08
0.12
0.16
0.2
(d) 1.4 1.2 1 0.8 0.6 0.4 0.2
Olazar et al. (2009) This Work 0.12
0.16
0.2
z [m]
0
Olazar et al. (2009) This Work 0
0.04
0.08
0.12
0.16
0.2
z [m]
(e)
(f)
AC
0
Olazar et al. (2009) This Work
z [m]
Ur
0.1
0
0
0.2
ED
(c)
0.16
M
0
0.5
Figure 5: Comparison between results obtained in this work, experimental data [48, 49] and simulated results reported by Olazar et al., 2009 [39]. (a) voidage (b) gas velocity in the spout (c) solids velocity in the spout (d) gas velocity in the annulus (e) solids velocity in the annulus (f) gas velocity at spout-annulus interface
20
ACCEPTED MANUSCRIPT
be observed a reasonable agreement between the results obtained in this work and the experimental data, similar as the results reported by Olazar et. al., 2009 [39]. On the other hand, Olazar et al., 2009 [39] stated that the considerable difference between the
CR IP T
predicted data and the experimental ones for the particle velocities in spout is due to the constant spout diameter considered in the simulations. However, any further study was conducted in this work to verify such a claim. Instead, the spout diameter was assumed constant in all simulations and optimizations, since a variable spout diameter would lead to the incorporation of highly nonlinear equations into the model. While this is not a problem for simulations, it could make difficult
AN US
convergence in optimizations.
3.2. Simulation of a spouted bed reactor for the OCM
3.2.1. Comparison between experimental and simulated data for a conventional spouted bed reactor A reference case was selected in order to check the performance of the proposed model and
M
verify its predictive ability through of a comparison of simulated results and experimental data. In
ED
this fashion, data of a laboratory-scale reactor from Mleczko and Marschall, 1997 [51] were chosen as a base case due to the limited availability in the open literature of data for the OCM in spouted
PT
bed reactors. This was the only published work in considering the OCM in a spouted bed reactor under realistic operating conditions. Mleczko and Marschall, 1997 [51] investigated experimentally
CE
the oxidative coupling of methane in various types of reactors including an internally circulating fluidized-bed, a bubbling fluidized-bed, and a conventional spouted bed.
AC
Table 3 contains the specifications and conditions used in this simulation. Values for the feed gas composition and geometric parameters of the spouted bed were taken from Mleczko and Marschall, 1997 [51] while particle specifications and the kinetic model are given by Stansch et al., 1997 [43]. Details of all the species considered in the OCM and included in the simulations, as well as the 21
ACCEPTED MANUSCRIPT
Table 3: Parameters used in the simulation of the OCM in a spouted bed reactor
Symbol
Value
Units
Height of conical part
zc
0.03
m
Atmospheric pressure
Patm
101325
Pa
dp
0.355
mm
g
9.8
Particle density
ρp
3600
m s2 kg m3
Inlet diameter
d0
0.01
m
Column diameter
dc
0.05
m
Stagnant bed height
H
0.035
m
dbc
3 · Di
m
εa
0.42
−
PCH4
57.10
PO2
22.9
PN2
20.0
CH4 /O2
2.5
−
T0,s
1073-1173
K
Particle diameter Gravity
Diameter at the column base Voidage in annulus
Methane to oxygen ratio Feed gas temperature
kPa
AN US
Feed gas partial pressures
CR IP T
Parameter
properties of such components can be found in Appendix B. The considered reactor is operated in
M
adiabatic mode and the superficial velocity of the feed gas is U0 = 1.1 · Ums . The equations of the
ED
model were solved by following the algorithm presented in section 2.2 since the initial conditions for temperature and composition at the base of the annulus are unknown, so the problem has to be solved as a boundary value problem that satisfies Ta (H) = Ts (H) at the exit of the reactor.
PT
Fig. 6 shows a comparison of the experimental methane conversion, C2+ product selectivity and
CE
C2+ product yield with the results obtained in the simulations for different temperature conditions. By analyzing Fig. 6, it is noted that the proposed model for a conventional spouted bed reactor,
AC
in which the OCM process takes place, presents a reasonable agreement for the three quantities analyzed when compared with experimental data throughout the considered temperature range, even although a slight overprediction is observed in the simulated results. Nevertheless, these results suggest that the model proposed in this work is representative of the spouted bed reactor
22
ACCEPTED MANUSCRIPT
XCH4
SC2+
YC2+
40
20
10
0 1093
1113
1133
1153
1173
AN US
1073
CR IP T
[%]
30
Temperature [K]
Figure 6: Comparison between simulated (lines) and experimental data (bars) for CH4 conversion, C2+ selectivity and C2+ yield. Experimental data taken from [51]
M
behavior.
ED
3.2.2. Simulation of a spouted bed reactor with a draft tube The behavior of a spouted bed reactor with a draft tube for the OCM is also investigated in this work. Thus, a reference case at a fixed feed gas temperature of 1113 K was defined to forecast
PT
how such a modification in a conventional spouted bed affects the performance of the OCM. The
CE
remaining specifications of this reference case were the same as those listed in Table 3. Fig. 7 shows theoretical predictions of the model in both temperatures in the annulus and
AC
spout along the axis for a conventional spouted bed reactor as well as for a spouted bed with draft tube. Temperature profiles present similar tendencies in both reactors. It is appreciated that the temperature in the spout has a mild increase close to the reactor inlet, follows by a remarkable reduction. This temperature gradient occurs in a very short reactor length of just a few millimeters.
23
ACCEPTED MANUSCRIPT
1150
1050
1000
950
CR IP T
Temperature [K]
1100
Ts Ta
900
0
0.5
1
1.5
2
2.5
3
·10−2
AN US
Reactor Length [m]
Figure 7: Temperature profiles in the spout and annulus for conventional spouted bed (solid lines) and spouted bed with draft tube (dashed lines)
The first increase is caused mainly by the heat released by the exothermic reactions that take place
M
in the OCM process while the abrupt reduction is because of the cross-flow of particles coming from
ED
the annulus, which are at a lower temperature. For the remaining reactor length, the temperature profile in the spout follows a smooth increase until it reaches a flat behavior. This flat profile at
PT
the last part of the reactor length for the temperatures in the annulus and spout is caused by the efficient mixing of the solids and gas in a spouted bed operation. It is also observed that the
CE
temperature in the annulus is initially lower than in the spout but, then this starts to rise and after a certain point presents a similar behavior to the temperature in the spout until both temperatures
AC
achieved the same value. This condition was imposed as a boundary condition. Another interesting point to mention is the moderate temperature hot spot occurred just after
the reactor inlet, which is much lower than the hot spots reported in others types of reactors. For example, unacceptable hot spots formation with temperature variation between 400 K and 700 24
ACCEPTED MANUSCRIPT
K have been reported in very small FBR [52, 29, 53, 54, 14] and even of 280 K close to the gas distributor in some fluidized bed reactors [55]. These elevated temperature gradients are regularly avoided by using a large nitrogen dilution in the feed gas or operating the reactor in cooling mode.
CR IP T
In this reference case, a nitrogen dilution of 20% was used, which is much lower than those dilutions commonly found in the literature to prevent that the temperature gets out of control (usually a feed dilution with nitrogen of more than 80% is used) [29, 52].
This numerical simulation demonstrates the theoretical feasibility of a spouted bed used as a chemical reactor for the OCM to control the hot spot formation in an acceptable range, operating
AN US
the reactor in adiabatic mode and using small nitrogen quantities in the feed so that complicated and expensive reactor design for cooling purposes could be avoided.
The differences between the temperature profiles in both types of reactors are a consequence of the performance achieved in each reactor, which was higher in the conventional spouted bed
M
reactor. The higher performance accomplished in the conventional spouted bed reactor can better be appreciated in the Fig. 8, which illustrates the concentration distribution profiles in the spout
ED
and annulus of all the considered species for the two types of reactors. Profiles in both reactors have similar tendencies. However, it is evident that these profiles are
PT
lower in the spouted bed reactor with a draft tube, indicating a lower performance of the OCM in this reactor. This lower performance is caused because the insertion of a device in the spout-
CE
annulus interface reduces the fluid flow and mass transfer in the annular region considerably. In consequence, the C2+ production rates and the methane consumption are also reduced, so that
AC
it can be concluded that the insertion of the draft tube in a spouted bed reactor, in which the OCM takes place, leads to mass transport limitations between the spout and annulus causing a considerable reduction of C2+ yields.
25
ACCEPTED MANUSCRIPT
6
7
6
5
5
4
4 3 3 2 2 1
1
0
0.5
1
1.5
2
2.5
Reactor Length [m]
0
3
(a)
0
0.5
1
1.5
2
Reactor Length [m]
·10−2
AN US
0
CR IP T
m3
mol Molar Concentration
CH4 O2 CO CO2 C2 H6 C2 H4 H2 O H2
2.5
3
·10−2
(b)
Figure 8: Concentration profiles in the spout (a) and annulus (b) region for conventional spouted bed (solid lines) and spouted bed with draft tube (dashed lines)
M
An important aspect to be noted from Fig. 8 is that better C2+ yields (higher concentrations of both ethane and ethylene) were achieved in the annular part of the conventional spouted bed.
ED
This trend can be explained by the distributed feed of the components along the spout-annulus interface caused by the cross-flow of fluid from the spout into the annulus. This means that oxygen
PT
concentrations are uniformly maintained in low levels at every point along the reactor in the annular part. The fact of keeping a distributed feed of oxygen in this region of the reactor favors the reaction
CE
conditions of the OCM, which has previously been demonstrated in some packed bed membrane reactors [56, 29]. Also, it can be seen that concentration distributions of the components remained
AC
almost constant in the annular region after the first half of the reactor because the oxygen is almost completely consumed in this part of the reactor.
26
ACCEPTED MANUSCRIPT
3.2.3. Effects of feed gas temperature and methane-to-oxygen ratio in the feed on the yield A very basic sensitivity analysis is presented in this section to examine the effects of some operating conditions on the methane conversion, C2+ selectivity, and C2+ yield, and thereby, to identify
CR IP T
some key parameters for further improvement of the reactor performance through optimizations of these parameters. To this end, several numerical simulations of the reference case were performed by varying one variable at a time to provide a first insight on the influence these parameters on the reactor performance, and on the optimal regions for them, in which the reactor reaches its highest performance. This does not represent a proper sensitivity analysis as more sophisticated methods
AN US
such as Monte Carlo simulations are needed, neither does not represent an adequate optimization, but help us to get some initial insights of the effect of some parameters on the reaction performance. For this analysis, the feed gas temperature and the methane-to-oxygen ratio were varied but keeping an undiluted feed and the other variables constant. The temperature was varied from 973 K
M
to 1173 K due to the range of applicability of the kinetic model, and the methane-to-oxygen from 2 to 10, respectively. The reason for using a minimum value of 2 for the methane-to-oxygen ratio
ED
is because some preliminary simulations indicated a total oxidation of the hydrocarbons to carbon dioxides with lower methane-to-oxygen ratios. A similar conclusion was reported by Tye et al., 2002
PT
[57]. For all of these reaction conditions, an oxygen consumption almost completely was achieved at the exit of the reactor. The effect of the feed gas temperature and the methane-to-oxygen ratio
CE
on the methane conversion, C2+ selectivity, and C2+ yield is displayed in Figures 9, 10 and 11. A series of several plot types are presented in each case to allow for a better understanding
AC
of the influence of these variables on the specified quantity. First, it is shown the classical 2-D representation plot for each independent variable (temperature and methane-to-oxygen ratio in the feed) to see the effect of each of these variables separately, and then a 3-D surface plot and a 2-D
27
ACCEPTED MANUSCRIPT
XCH4 [%]
30 25 20 15
=2 = 2.5 =3 =4 =5 =6 =7 =8 =9 = 10
40
T = 973 [K] T = 993 [K] T = 1013 [K] T = 1033 [K] T = 1053 [K] T = 1073 [K] T = 1093 [K] T = 1113 [K] T = 1133 [K] T = 1153 [K] T = 1173 [K]
35 30 25 20 15
10
10
5 973
5 1003
1033 1063 1093 1123 Temperature [K]
1153
2
(a)
4
5 6 7 CH4 /O2 ratio
8
9
9
35
8
30
AN US
7
25 15
25
6 5
20
4
5 3
4 5 6 CH4 /O 2
7 8 ratio
9
10
1153 ] 1113 [K re 1073 tu a r 1033 e 993 mp Te
15
3 2
993
1033 1073 1113 Temperature [K]
10 1153
(d)
M
(c)
XCH4 [%]
CH4 /O2 ratio
35
2
10
(b)
10
XCH4 [%]
3
CR IP T
CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2
35
XCH4 [%]
40
ED
Figure 9: Effect of reaction temperature and methane to oxygen ratio in the feed over CH4 conversion. (a) 2-D plot of the effect of T over CH4 conversion at several levels of CH4 /O2 (b) 2-D plot of the effect of CH4 /O2 over CH4 conversion at several levels of T (c) 3-D surface plot of the combined effect of T and CH4 /O2 over CH4 conversion
PT
(d) 2-D contour plot of the combined effect of T and CH4 /O2 over CH4 conversion
CE
contour plot are also given as alternatives to see the combined effect and the interaction between these two independent variables.
AC
Fig. 9 shows the variations of the methane conversion with the feed gas temperature and methane-to-oxygen ratio. More specifically, Fig. 9 (a) presents the effect of the temperature on the methane conversion at different levels of methane-to-oxygen ratio, while Fig. 9 (b) presents the effect of the methane-to-oxygen ratio in the feed on the methane conversion at various levels of 28
ACCEPTED MANUSCRIPT
70
SC2+ [%]
60 50 40
=2 = 2.5 =3 =4 =5 =6 =7 =8 =9 = 10
80
60 50 40
30 20 973
T = 973 [K] T = 993 [K] T = 1013 [K] T = 1033 [K] T = 1053 [K] T = 1073 [K] T = 1093 [K] T = 1113 [K] T = 1133 [K] T = 1153 [K] T = 1173 [K]
70
30 20 1003
1033 1063 1093 1123 Temperature [K]
1153
2
(a)
8
9
10
60
8
50 40 30 20 10 2
3
4
5
C
7
6
9
8
2 /O
t ra
H4
io
50
6
45
5
40 35
4
30
3 2
993
1033 1073 1113 Temperature [K]
25 1153
(d)
M
(c)
10
55
7
AN US
60
SC2+ [%]
CH4 /O2 ratio
70
SC2+ [%]
5 6 7 CH4 /O2 ratio
65
9 80
1033 1073 1113 1153 Temper ature [K ]
4
(b)
10
993
3
CR IP T
CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2
SC2+ [%]
80
ED
Figure 10: Effect of reaction temperature and methane to oxygen ratio in the feed over C2+ selectivity. (a) 2-D plot of the effect of T over C2+ selectivity at several levels of CH4 /O2 (b) 2-D plot of the effect of CH4 /O2 over C2+ selectivity at several levels of T (c) 3-D surface plot of the combined effect of T and CH4 /O2 over C2+ selectivity
PT
(d) 2-D contour plot of the combined effect of T and CH4 /O2 over C2+ selectivity
CE
temperature. As can be noticed from Fig. 9 (a), the methane conversion shows the same trend with increasing temperature for all methane-to-oxygen ratios, passing through a maximum value, which
AC
is much more notable and higher at the lower levels of methane-to-oxygen ratio due to a higher oxygen concentration at these conditions. Fig. 9 (b) shows that the methane conversion decreases almost exponentially when the methane-to-oxygen ratio is increased at a given temperature. It is also noted in Fig. 9 (c-d) that the methane conversion has an optimal point throughout the range of 29
ACCEPTED MANUSCRIPT
YC2+ [%]
13 11 9 7
=2 = 2.5 =3 =4 =5 =6 =7 =8 =9 = 10
15
T = 973 [K] T = 993 [K] T = 1013 [K] T = 1033 [K] T = 1053 [K] T = 1073 [K] T = 1093 [K] T = 1113 [K] T = 1133 [K] T = 1153 [K] T = 1173 [K]
13 11 9 7
5
5
973
1003
1033 1063 1093 1123 Temperature [K]
1153
2
(a)
4
5 6 7 CH4 /O2 ratio
8
9
13
9
12
8
11
AN US
7
11 9 7
10
6
9
5
8 7
4
5 3
4 5 6 CH4 /O 2
7 8 ratio
9
10
1153 ] 1113 [K re 1073 tu a r 1033 e 993 mp Te
993
1033 1073 1113 Temperature [K]
5 1153
(d)
M
(c)
6
3 2
YC2+ [%]
CH4 /O2 ratio
13
2
10
(b)
10
YC2+ [%]
3
CR IP T
CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2
YC2+ [%]
15
ED
Figure 11: Effect of reaction temperature and methane to oxygen ratio in the feed over C2+ yield. (a) 2-D plot of the effect of T over C2+ yield at several levels of CH4 /O2 (b) 2-D plot of the effect of CH4 /O2 over C2+ yield at several levels of T (c) 3-D surface plot of the combined effect of T and CH4 /O2 over C2+ yield (d) 2-D contour plot
PT
of the combined effect of T and CH4 /O2 over C2+ yield
CE
temperatures and methane-to-oxygen ratios investigated. This optimal point (methane conversion of 39.15 %) is reached on the minimum methane-to-oxygen ratio used and a temperature of 1106 K.
AC
For C2+ selectivity similar analysis can be done by analyzing Fig. 10, which shows the effect of feed gas temperature and the methane-to-oxygen ratio on the C2+ selectivity. In this case, it is also noticed that the selectivity profile shows a maximum point over the temperature range for all of the methane-to-oxygen ratio levels (see Fig. 10 (a)). This behavior is steeper at the lower levels 30
ACCEPTED MANUSCRIPT
of methane-to-oxygen ratio. Fig. 10 (b) shows that C2+ selectivity increases almost linearly with an increase of the methane-to-oxygen ratio for a given temperature. The effect that the methaneto-oxygen ratio has on the selectivity is opposite to its effect on the conversion, which decreases
CR IP T
with increasing the methane-to-oxygen ratio. Figures 9 (c-d) illustrate an optimal point for C2+ selectivity (70.58 %), which is located at a temperature of 1112 K and on the highest methaneto-oxygen ratio used. However, it is worth noting that the lowest methane conversions were also found at this point, making evident that an intermediate point is required to provide the highest conversions and at the same time a high selectivity.
AN US
Fig. 11 exhibits that C2+ yield has a similar behavior to the methane conversion when varying the feed gas temperature and the methane-to-oxygen ratio but with a much sharper trend, that is, the C2+ yield enhances with increasing temperature until reaching a maximum value and after this point, the bigger the temperature, the lower the C2+ yield. This behavior was the same for all
M
methane-to-oxygen ratios, and again as in the case of the methane conversion, the maximum temperature behavior was much more prominent when decreasing the methane-to-oxygen ratio. This
ED
fact shows that the highest C2+ yields were reached over the points of higher methane conversions but not in those of higher ethylene selectivities. The optimal point achieved for C2+ yield was of
PT
14 % at a feed gas temperature of 1118 K and a methane-to-oxygen ration of 2. The low methane conversions achieved by increasing the methane-to-oxygen ratio is explained
CE
by the kinetic mechanism and stoichiometry (see Appendix C), which suggest that high oxygen concentrations in the feed favor the methane conversion. On the other hand, the rise of the C2+
AC
selectivity with the feed gas temperature before reaching its maximum value results from the bigger activation energy of the selective, primary reaction step (equation C.2) compared to the nonselective ones (equations C.1 and C.2) [14].
31
ACCEPTED MANUSCRIPT
The drop of the C2+ selectivity after its maximum value is caused likely because some unwanted reactions such as oxidation and steam reforming of ethylene (equations C.6 and C.8) are favored at higher temperatures. Moreover, the rise of the C2+ selectivity with increasing the methane-to-
CR IP T
oxygen ratio is because the reaction order with respect to oxygen is lower for the selective, primary reaction step (equation C.2) than for the non-selective ones (equations C.1 and C.2), while the reaction order with respect to methane is higher for the former [58].
All of these aforementioned trends were found in the simulations of the spouted bed reactor and agreed with experimental observations as well as with several simulation work conclusions
AN US
commonly reported in the literature for various types of reactor and catalysts for the OCM [59, 60, 57, 61, 58, 14].
3.3. Optimization of a spouted bed reactor for the OCM
Several operational and geometrical parameters of the spouted bed reactor can be optimized,
M
and this section starts with the optimization including only two decision variables or free variables
ED
that can be adjusted during the optimization process to enhance the performance of the reactor. These variables are the feed gas temperature and the methane-to-oxygen ratio. These two variables
PT
were chosen because a basic sensitivity analysis was done before, varying each of them using simple simulations. So, a first insight is known about where the optimal point can be located whether the
CE
remaining variables are kept constant. The addition of others decision variables to the optimization study is explored in the second part of this section. Finally, the last part presents the optimization
AC
results for a modified spouted bed reactor with a draft tube, in which feed streams of methane and oxygen are split by employing a secondary feed zone of fluid in the annular region of the reactor. In this way, a distributed feed of oxygen is accomplished in the annular zone as a great solution to the low C2+ yield achieved in the conventional configuration of the reactor. 32
ACCEPTED MANUSCRIPT
3.3.1. Optimization of a conventional spouted bed reactor Inequality constraints for decision variables of this first optimization were as follows: 973 ≤ Tg,0 ≤ 1173 and 1 ≤ PCH4/PO2 ≤ 10. In addition to these bounds for the feed gas temperature and the methane-
CR IP T
to-oxygen ratio in the feed, the following bounds were introduced to avoid evaluation errors in some constraints: • All concentrations should be smaller than 20
mol m3
• XCH4 , SC2+ , and YC2+ should be between 0 and 1
AN US
• Partial pressures should be smaller than 106 Pa
• All specific heat capacities should be greater than 10 and smaller than 106 • All mixture averaged molar masses should be greater than 1
J mol·K
kg mol
M
• All dynamic viscosities should be greater than 10−8 and smaller than 1 Pa · s W m·K
ED
• All thermal conductivities should be greater than 10−6
The result of this optimization is shown in Fig. 12 alongside the optimal point found by simply
PT
changing each of these variables at a time as performed in section 3.2.3. It is noted that both results are consistent. A slightly quantitative difference is observed regarding the optimal temperature and
CE
the maximum yield found; this small difference is expected because the simulations were performed in a temperature range of 973 − 1173 K but not every point in this range was included in the
AC
simulations. This small difference can also be caused by the numerical precision of the solver in the different stages in simulations and optimizations since different solver options were adjusted for each case to achieve convergence.
33
ACCEPTED MANUSCRIPT
1118.45
Optimizations results 1121.5
20
14.8
4
5
Temperature
CH4 O2
ratio
1033
2
1003
1
973
0
AN US
C2+ Yield
3
CR IP T
2 2
10
0
1093
1063
[%]
14
5
CH4 /O2 ratio
15
1123
Temperature [K]
Simulations results 25
Figure 12: Comparison between optimal yield found by simulations (blue bars) and optimizations (red bars)
3.3.2. Improving the reactor performance by adding decision variables in optimizations Table 4 contains a list of some key parameters of the conventional spouted bed that may be
M
changed within given bounds. These should not be understood as an exhaustive list of all possible decision variables in a spouted bed. Several other potential decision variables exist. Among those
ED
are the inlet pressures and the type of catalyst. For matters of problem size and complexities in the resulting optimization problem, these were disregarded here. The lower and upper bounds were
PT
chosen after a systematic literature search about the effect of these variables on the OCM process and the spouted bed operation.
CE
For example, bounds for the feed gas temperature and methane-to-oxygen ratio in the feed were selected according to the validity range of the kinetic model [43]. Bounds for dilution with
AC
nitrogen in the feed were defined based on Tiemersma, 2010 [56] and Jaˇso et al., 2010 [62]. Bounds for catalyst properties (density and diameter of catalyst) were chosen according to Esche et al., 2012 [29] and Lenhart, 2010 [63]. Bounds on hydrodynamic parameters were selected taking into
34
ACCEPTED MANUSCRIPT
Table 4: Decision variables for the conventional spouted bed reactor
Parameter
Symbol
Initial value
Lower bound
Upper bound
Units
Operating parameters T0,g
1000
973
1173
[K]
5
1
10
[−]
40
0
80
[%]
U0 Ums
1.1
1
2
[−]
1500
3700
0.3
2
0.05
0.5
[−]
0.07
[m]
CH4 O2
Methane-to-oxygen ratio in the feed
0
N2 |0
Dilution with nitrogen in the feed Feed gas superficial velocity
Design parameters Particle density
ρp
3600
Particle diameter
dp
0.355
Ratio of the fluid inlet diameter to the
d0 dc
0.2
AN US
cylinder diameter
CR IP T
Feed gas temperature
Stagnant bed height
H
0.035
0.02
h
kg m3
i
[mm]
account the operating and geometrical conditions that must be met to achieve a stable spouted bed operation [64].
M
Extra inequality constraints should be considered when some geometrical parameters of a spouted bed device are used as decision variables. While changing these geometrical parameters
ED
during an optimization, some criteria must be satisfied to guarantee a stable operation of spouted
d0 ≤ dc
Umf Ut
1/2
≈
d0 n/2 ≤ εmf dc
(24)
CE
PT
bed. According to Epstein and Grace, 2011 [64], these constraints are as follows:
where d0 is the inlet diameter, dc is the column diameter, Umf is superficial gas velocity at minimum
AC
fluidization, Ut is the terminal gas velocity, and εmf is the voidage at minimum fluidization. n is given in function to the dimensionless Archimedes number as follows [64]: 4.8 − n = 0.047 · Ar0.57 n − 2.4 35
(25)
ACCEPTED MANUSCRIPT
Other two inequality constraints used regarding to the geometrical parameters were:
d0 ≤ 30 dp
where dp is the particle diameter.
CR IP T
dc 30 ≥ n/2 dp εmf
(26)
(27)
AN US
When the above constraints are introduced in the optimization problem, this latter became harder to solve, presenting some instabilities problems and locally infeasible points. The results presented in this section were not obtained by adding all of the decision variables at once to the optimization problem. Instead, these were added to the optimization slowly increasing the degrees of freedom sequentially until all the decision variables were freed. The decision variables were added
M
to the optimization as follows: first, the only two decision variables were the feed gas temperature
ED
and the methane-to-oxygen ratio in the feed as can be seen in previous section. Afterward, the remaining decision variables were added in the following sequence: dilution with nitrogen in the
PT
feed, feed gas superficial velocity, ratio of the fluid inlet diameter to the cylinder diameter, stagnant bed height, particle diameter and particle density. This order choice was arbitrary and, certainly,
CE
the order of adding the decision variables could affect the search for the final optimum, and this can lead to a local optimum, but any further study was conducted in this work to investigate such
AC
an effect on the optimal results. A high number of intermediate steps were performed to achieve the final result for the opti-
mization process, but only five of these intermediate steps will be presented for practical purposes and for giving a general idea of how this process was done until an optimal solution was reached. 36
ACCEPTED MANUSCRIPT
Table 5: Progress of the decision variables in several stages of the optimization process
T0,g LB
CV
UB
LB
1
-
1000
-
2
990
1010UB
3
985
4 5
CV
UB
LB
CV
UB
LB
CV
UB
-
5
-
-
40
-
-
1.1
-
1010
4.0
4LB
6.0
30
50UB
50
1.0
1.2UB
1.2
1040UB
1040
3.0
3LB
7.0
20
52F
975
1080UB
1080
2.0
2.8F
8.0
10
54F
973
1109F
1173
1.0
2.9F
10
0
54F
ρp
No
U0 Ums
N2 |0
d0 dc
dp
CV
UB
LB
CV
UB
LB
CV
1
-
3600
-
-
0.355
-
-
0.2
2
2500
3650UB
3650
0.34
0.5UB
0.5
0.15
3
2000
3660UB
3660
0.33
0.92F
1.0
4
1000
3680UB
3680
0.31
1.0F
5
1500
3700UB
3700
0.3
1.0F
60
1.0
1.2F
1.6
70
1.0
1.38F
1.8
80
1.0
1.38F
2.0
H
UB
LB
CV
UB
-
-
0.035
-
0.3UB
0.3
0.030
0.04UB
0.04
0.1
0.35UB
0.35
0.025
0.05UB
0.05
1.5
0.07
0.35F
0.45
0.220
0.06UB
0.06
2.0
0.05
0.35F
0.5
0.02
0.65F
0.07
AN US
LB
CR IP T
No
CH4 O2 0
M
LB=lower bound; UB=upper bound; CV=current value; F=free variable without bounds active
ED
Whenever a decision variable is at its current lower or upper bound, or at neither of them an index LB, UB or F indicates this respectively. Table 5 shows results for values of decision variables
PT
examined at several stages of the whole optimization procedure until the optimal value was found. Table 6 presents the progress of the methane conversion, C2+ selectivity and C2+ yield at each of
CE
these intermediate steps.
It can be observed that the optimal feed gas temperature and methane-to-oxygen ratio were
AC
slightly different from the optimal values found before, but this differences can be caused due to the interaction of these two variables with others of the included decision variables. Also, there is a trend of using a relatively high amount of dilution with nitrogen, the optimal point was 54%, but
37
ACCEPTED MANUSCRIPT
Table 6: Evolution of the methane conversion, C2+ selectivity and C2+ yield in several stages of the optimization process
C2+ yield [%]
C2+ selectivity [%]
Conversion of CH4 [%]
1
7.38
70.3
10.5
2
10.12
67.0
15.1
3
14.3
64.4
4
18.75
62.1
5
21.86
58.3
CR IP T
No
22.2
30.2 37.5
also some optimizations were carried out by excluding this decision variable, and the yield does not
AN US
drop very much. It was also noted that the optimal inlet superficial velocity is higher than that initially used.
In relation to the catalyst properties, it was found that the optimal particle diameter was in the range of particles classified as Geldart G, which is easier spoutable (1-5 mm). However, an
M
assumption in the model was to disregard fluid-solid mass transfer resistance in the catalyst, and this assumption could cause over-predicted results, which could be more remarkable for larger
ED
particle diameter as those considered in this work. In this way, it would be important to consider an effectiveness factor to account the mass transfer resistance.
PT
More difficult to interpret and less pronounced are effects of catalyst density. Increasing this seems to improve the performance of the reactor, but it is unknown to what extent the catalyst
CE
density can be increased beyond the specified 3600 kg/m3 , since to modify the catalyst density, it is needed to modify the catalyst itself. A value of 3700 was assumed here to be the upper limit,
AC
and at the final stage of the optimization process, the upper bound was still active for this decision variable. Optimal points for geometrical parameters were also found successfully, but the range of these decision variables was not so broad due to the difficulties arise when the additional inequality
38
ACCEPTED MANUSCRIPT
constraints were added to the optimization. It has to be said that between steps four and five (although more intermediate steps were done in between them) the changes made to the variable bounds had to be chosen very carefully,
CR IP T
and the increases consequentially became even smaller to avoid infeasibilities. Finally, in Table 6 can be verified that the optimal C2+ yield found in a conventional spouted bed reactor increase considerably when others decision variables are included in the optimization compared with that previously calculated with only two decision variables.
AN US
3.3.3. Improving the reactor performance by modifying the conventional spouted bed
This section explores the possibility of inserting a porous membrane as a draft tube in the conventional spouted bed configuration to increase the fluid flow rate from the spout region to the annular region. Also, a secondary feed zone of fluid with low inlet velocity in the annular region of the reactor is contemplated. For this configuration, the oxygen is fed into the spout and the
M
methane into the annulus. In this way, it is possible to achieve a distributive feeding of oxygen in
ED
the annulus where the reaction will take place.
Two concepts of reactors will be then combined in order to enhance the performance of the OCM process. The flexibility of a membrane reactor concept will be combined with the excellent
PT
mass and heat transfer properties found in spouted bed reactors. A sketch of this reactor concept
CE
is shown in Fig. 13. It is important to highlight that, although the proposed concept will be theoretically studied here with the aid of optimizations to find useful operational and geometrical
AC
parameters to operate and design the proposed reactor, future work to prove the practical viability of this reactor concept will be necessary. A slight modification needs to be done to the balance equations to account for the modification
included in the spouted bed reactor (membrane as a draft tube). This modification includes the
39
CH4
AN US
O2
CR IP T
ACCEPTED MANUSCRIPT
Figure 13: Spouted-fluid bed with porous draft tube
mass flux in the annular region due to the fluid flow through the membrane tube. Ideally, dense membranes are preferred, but here a porous membrane is considered because the high complexity
M
in modeling dense membranes. The porous membrane modeled here is a γ − Al2 O3 membrane. This particular membrane has been applied in simulation studies of packed bed membrane reactors
ED
[56, 65].
Generally, in a porous membrane all of the components can permeate through the membrane
PT
because this type of membrane is not highly selective, but in this work, it was considered in advance and without any reasonable theoretical fundamental that the only component able to permeate the
CE
membrane is oxygen. In this way, if the pore diameter of these membrane tubes is small enough, the diffusion of oxygen through the membrane can be described by the Dusty-Gas-Model, which
AC
incorporates a contribution of Knudsen diffusion and viscous flow [56].
JO2
1 = R·T
4 K0 3
s
8·R·T P + B0 π · MWO2 µg 40
!
·
4Pi,s−a m (rdt + δm ) · ln rdtr+δ dt
(28)
ACCEPTED MANUSCRIPT
where JO2 is the molar flux of oxygen through the membrane, K0 and B0 are parameters of the model, 4Pi,s−a is the partial pressure difference between the spout and annulus of component i, in this case, oxygen, rdt is the radius of the draft tube, and δm is the membrane thickness.
CR IP T
As it can be seen with this model, the distribution of oxygen through the spout-annulus interface is governed by the partial pressure difference (driving force) and the membrane permeability. Because the exact performance of the membrane is of less interest, and the effect of this on the spouted bed reactor is more important, the above model was considered adequate as a first attempt
AN US
to model a spouted bed reactor with a membrane draft tube, so a constant flux of oxygen was used (equation 28). However, it is necessary to evaluate the performance of the membrane as well as the influence of using a variable oxygen flux to extend the application of this concept. Heat and mass transfer correlations describing transport from the catalyst bed to the membrane wall were taken
M
from literature [56].
Besides the decision variables considered in the optimization case for the conventional spouted
ED
bed reactor, for this optimization, others decision variables were included. These new decision variables were: the draft tube diameter with its lower and upper bounds measured in relation to
PT
the inlet diameter (1.5 · d0 ≤ ddt ≤ 0.5 · d0 ), the distance from the inlet of the draft tube with its lower and upper bounds (0 ≤ h0−dt ≤ H), the height of the draft tube with its lower and upper
CE
bounds (0 ≤ hdt ≤ H), and the superficial velocity of the auxiliary feed in the annulus, which was measured with respect to the minimum fluidization velocity to guarantee a stable spouted bed
AC
operation [64]. The limits of this were 1 ≤ U0,a /Umf ≤ 5. An additional constraint was introduced in this case to avoid that the height draft tube and the
distance from the inlet of the draft tube exceed the stagnant bed height (hdt + h0−dt ≤ H). In this
41
ACCEPTED MANUSCRIPT
Table 7: Results for the optimization of the spouted bed reactor with a porous draft tube Reactor performance
Optimal decision variables
Methane conversion
XCH4
45.4
[%]
Feed gas temperature for both annulus and spout
C2+ selectivity
SC2+
73.2
[%]
Methane-to-oxygen ration in the feed
C2+ yield
YC2+
33.2
[%]
Dilution with nitrogen in the feed in the annulus
Particle density Particle diameter
1119.5
[K]
1.5
[−]
N2 |0
51
[%]
U0 Ums
1.45
ρp
3700
dp
1.2
0
[−] h i kg m3
CR IP T
Feed gas superficial velocity
T0,g
CH4 O2
[mm]
Ratio of the fluid inlet diameter to the cylinder diameter
d0 dc
0.4
[−]
Stagnant bed height
H
0.07
[m]
ddt
1.05 · d0
[m]
h0−dt
0.5
[cm]
hdt
3
[cm]
U0,a Umf
3.5
[−]
Draft tube diameter
Distance from the feed inlet of the draft tube Height of the draft tube
AN US
Superficial velocity of the auxiliary feed in the annulus
optimization study, the feed gas temperature in the annulus and the spout were always the same (T0,g ), but could be different increasing then the degrees of freedom. Another decision variables as membrane properties and feed gas temperature in the annulus can also be included, but this is
M
subject to future work. Table 7 presents the results of the optimization process at the last stage, i.e., when the bounds were inactive for decision variables.
ED
Fig. 14 shows the optimal performance of the three different spouted bed reactors considered, the conventional spouted bed (CSB), the spouted bed with a draft tube (SBDT) and the spouted-
PT
fluid bed with a membrane draft tube (SFBMDT). It is clear that the reactor concept proposed here is superior in performance than the conventional spouted bed reactor optimized in the previous
CE
section, and the spouted bed with draft tube was the reactor with lower performance. The proposal of using a porous membrane as a draft tube seems to be a helpful solution, at least in theory, to
AC
enhance the reactor performance. However, further investigation must be conducted on this reactor concept to prove its validity.
42
ACCEPTED MANUSCRIPT
80
Conversion XCH4 h i Selectivity SC2+ h i Yield YC2+
70 60
40 30 20 10 0 CSB
SBDT
SFBMDT
AN US
Temperature [K]
CR IP T
[%]
50
Figure 14: Comparison between performance of the three spouted bed reactors optimized. CSB - Conventional spouted bed. SBDT - Spouted bed with a draft tube. SFBMDT - spouted-fluid bed with a membrane draft tube
M
4. Conclusions
A simultaneous approach for solving a DAE-constrained optimization problem was used and
ED
applied to a spouted bed reactor for ethylene production from methane to improve the performance of the oxidative coupling of methane. This methodology was able to provide reliable results when
PT
some predicted data were compared with results reported in the literature, verifying the efficiency of the full discretization approach for a DAE-constraint optimization problem. However, it is
CE
important to mention that the model was unable to predict well the solids velocities in the spout, this can be caused by the constant diameter assumed in spout, but some further investigation by
AC
using more rigorous modelling including CFD models must be performed to verify this, since radial variations of velocities or shear stress, which were not taken into account in the modeling, could also affect the solids velocities in the spout.
43
ACCEPTED MANUSCRIPT
On the other hand, the OCFE method implemented for discretization was computationally efficient in terms of CPU time. Therefore, the optimization approach implemented in this work has great potential for design applications and can be successfully applied to theoretical studies of
CR IP T
highly non-linear and complex engineering problems to improve chemical processes. Regarding the non-conventional spouted bed proposed, it can be concluded that the theoretical feasibility of combining several reactor concepts for a spouted bed operation to enhance the performance of the oxidative coupling of methane was demonstrated by means of detailed numerical simulations and optimizations, provided that some modifications of the conventional configuration
AN US
be carried out. Among these modifications, it can be highlighted, the use a porous membrane as a draft tube and the addition of a secondary feed zone of fluid with low inlet velocity in the annular region. These two modifications allow for a distributive feeding of oxygen in the annulus where the reaction takes place. Also, it was demonstrated the feasibility of using a spouted bed as a chemical
M
reactor for the OCM to control the hot spot formation in an acceptable range, operating the reactor in adiabatic mode and using small nitrogen quantities in the feed so that complicated and expensive
PT
Acknowledgments
ED
reactor design for cooling purposes could be avoided.
Financial support of this research was provided by S˜ ao Paulo Research Foundation (FAPESP),
CE
grant numbers 2011/22800-3 and 2013/14896-6, and by National Council for Scientific and Tech-
AC
nological Development (CNPq), grant number 141009/2012-8.
44
ACCEPTED MANUSCRIPT
Nomenclature
Description
BVP
Boundary Value Problem
DAE
Differential–Algebraic Equation
DV
Decision Variables
FBR
Fixed Bed Reactor
GAMS
General Algebraic Modeling System
IPOPT
Interior Point OPTimizer
IVP
Initial Value Problem
LB
Lower Bounds
ncr
number of catalytic reactions
ngr
number of gas reactions
NLP
Nonlinear Programming Problem
OCFE
Orthogonal Collocation on Finite Elements
OCM
Oxidative Coupling of Methane
UB
Upper Bounds
˙ G
PT
˙ Q rr ˙
0
AC
A
00
reaction rate per unit volume in gas reactions mol · m−3 · s−1 reaction rate per mass of catalyst in catalytic reactions mol · g−1 · s−1 mass flow rate of particles kg · s−1 surface area of solids per volume of bed m2 specific surface areas of solids m2 cross–sectional area m2
CE
rr ˙
A
M
component rate mol · m−3 · s−1 mass flow rate of gas phase kg · s−1 volumetric flow rate m3 · s−1
cr ˙
˙ W
ED
Roman Symbols
AN US
Symbol
A
B0 and K0
CR IP T
Acronyms
parameters of the Dusty-Gas- Model
45
ACCEPTED MANUSCRIPT
CD
drag coefficient heat capacity J · mol−1 · K−1 concentration mol · m−3
Cp c
Ea,j F g h J k01,3 k02 K0H2 O K0i KH2 O Ki l
activation energy of reaction j J · mol−1
force acting over differential element [N] acceleration of gravity m · s−2 heat transfer coefficient W · m−2 · K−1 molar flux mol · m−2 · s−1
preexponential factor of rate coefficient, k1,3 mol · bar1/2 · g−1 · h−1 preexponential factor of rate coefficient, k2 mol · bar−1 · g−1 · h−1
preexponential factor of dissociative adsorption constant KH2 O preexponential factor of adsorption constant, KCH4 ,CO,H2 bar−1
dissociative adsorption constant of H2 O
adsorption constants for CH4 , CO and H2 bar−1
length [m]
Pi
partial pressure of component i [Pa]
P
pressure [Pa]
M
m
molecular weight kg · mol−1
MW
ED
mass [kg]
ideal gas constant J · mol−1 · K−1
r
radius [m]
S
PT
R
selectivity in C2 hydrocarbons
T
temperature [K]
X
methane conversion
Y
yield in C2 hydrocarbons
u V
AC
v
CE
w
superficial gas velocity m · s−1 average gas velocity m · s−1 superficial particle velocity m · s−1 average particle velocity m · s−1
U
CR IP T
diameter [m]
AN US
d
mass fraction
46
ACCEPTED MANUSCRIPT
Greek Symbols
δ
differential distance [m]
γ
cone angle [◦ ]
thermal conductivity W · m−1 · K−1
λ
dynamic viscosity [Pa · s] density kg · m−3
µ ρ
residence time [s]
4G0f 4Hads,i 4Hr
4H0f 4z
Gibbs energy of formation J · mol−1 adsorption enthalpy of component i J · mol−1 enthalpy of reaction J · m−3 · s−1 standard enthalpy of formation J · mol−1 height of differential element [m]
ε
voidage or porosity
Subscripts inlet and initial condition
a
annulus
atm
atmospheric
B
Gidaspow model
b
bouyant
bc
base of the column
c
column
cat
catalyst
D
drag
ED
PT
AC
g
CE
dt
M
0
dp
AN US
τ
horizontally projected particle
draft tube gas
g-p
gas–particles
i
components
j
reactions
m
maximum
CR IP T
fluid–particle friction coefficient in drag model kg · m−3 · s−1
β
47
membrane
mf
minimum fluidization
ms
minimum spouting
o
reference state
p
particles
p-g
particles–gas
r
radial
s
spout
s-a
spout–annulus interface
STC
standard temperature condition
w
wall
AC
CE
PT
ED
M
AN US
m
CR IP T
ACCEPTED MANUSCRIPT
48
ACCEPTED MANUSCRIPT
Appendix A. Flowchart of algorithms used in simulations and optimizations Start – Define operating conditions – Define geometric parameters
Guess Ta0
CR IP T
– Calculate thermophysical properties
– Ocfe Method
Solve ci,a0 = ci,s0 + τg0 · cr ˙ i,a
– Conopt Solver
AN US
Set up the initial value pro-
– Gams Script
blem → Ts0 , Ta0 , ci,s0 and ci,a0 – False position method (new Ta0 ) – Gams Script
Solve the Dae system composed of
– Ocfe Method
mass, momentum and energy balances
– Ipopt Solver
Calculate the temperature difference
M
at the reactor exit: Ta (H) − Ts (H)
CE
PT
ED
No
abs (Ta (H) − Ts (H)) ≤ 5 K
Yes
Solution of the Ivp Initial guesses of
Set up the Bvp → Ts0 , ci,s0 as initial
variables for Bvp
conditions and; ci,a0 = ci,s0 + τg0 · cr ˙ i,a and Ta (H) = Ts (H) used as equality constraints
– Ocfe Method
AC
– Gams Script
Solve the Bvp
– Ipopt Solver Solution of the Bvp
Figure A.1: Flowchart of the GAMS algorithm to solve numerically a conventional spouted bed reactor model
49
ACCEPTED MANUSCRIPT
Start Define № of interior collocation points (Icp) Divide the domain into n Fe with n = 1, 2..., j
CR IP T
Select j = 1
Use initial conditions as initial guesses for
Increase n to a higher value
all state variables on the Icp for the jth Fe Solve the problem for j = 1
No
AN US
Optimal solution found?
Yes
Do j = j + 1
Save solution for the (j − 1)th Fe
M
Use solution of the end bound for the
(j − 1)th Fe as initial guess for all state variables on the Icp for the jth Fe
ED
Keep as initial guesses for all state varia-
bles for the (j − 1)th Fe its own solutions
PT
Solve again the problem for the jth Fe
AC
CE
No
End of the program
Optimal solution found?
Yes No
j≤n
Yes
Figure A.2: Flowchart of the algorithm for the initialization routine used in simulations
50
ACCEPTED MANUSCRIPT
Start Define decision variables (Dv)
Activate all Dv Set up maximum and minimum values for upper and lower bounds of the Dv ⇒ Ubmax and Lbmin
Lbmin ≤ Dvk ≤ Ubmax k = 1, 2..., No of Dv
CR IP T
Assign constant values to all Dv
Define realistic bounds for the rest of variables included in the model
Set up tight bounds for all Dv (4 = ±5% of
Run a fully func-
the fixed values used in the previous simulation)
tional simulation
Activate the optimizer by adding a scalar objective function ⇒ YC2+
Valid solution found?
AN US
No
Ubcurrent = constant value of Dv (1 + 4) Lbcurrent = constant value of Dv (1 − 4)
Yes
Set up initial guesses for all
Save current solu-
variables and multipliers
tion for all variables
Activate solver options
Increase Ubmax and Lbmin for
PT
ED
No
CE
Optimal solution found
AC
End of the program
Ipopt.opt file
Run an optimization
M
all Dv violating the condition
to speed up convergence
Ubcurrent ≤ Ubmax Lbcurrent ≥ Lbmin Yes Save current solution for all variables and/or multiplier bounds
No
Check if bounds of all Dv are active Yes Ubcurrent = Ubcurrent + 4 Lbcurrent = Lbcurrent − 4
Figure A.3: Flowchart of the algorithm to carry out optimizations on the spouted bed reactor model
51
ACCEPTED MANUSCRIPT
Appendix B. Thermodynamic properties
Table B.1: Chemical species used in this work. Molecular weights were taken from [66]
Index
Formula
Name
MW [kg/kmol]
1
CH4
Methane
2
O2
Oxygen
3
CO
Carbon monoxide
4
CO2
Carbon dioxide
5
C 2 H6
Ethane
6
C 2 H4
Ethylene
7
H2 O
Water
8
H2
Hydrogen
2.016
9
N2
Nitrogen
28.01
CR IP T
16.04 31.40
28.01
44.01
30.07
28.05
AN US
18.02
Table B.2: Standard enthalpy and Gibbs energy of formation. Source: [67]
CH4
4H0f [kJ/mol]
M
Component
-74.86
-50.49
0.00
0.00
CO
-110.62
-137.15
CO2
-393.77
-394.37
C2 H 6
-84.74
-31.92
C2 H 4
52.33
-68.44
-242.00
-228.59
H2
0.00
0.00
N2
0.00
0.00
PT
ED
O2
H2 O
AC
CE
4G0f [kJ/mol]
52
ACCEPTED MANUSCRIPT
Table B.3: Ideal gas heat capacity coefficients. Source: [68]
B
C
D
CH4
19.251
5.213e-02
1.197e-05
-1.132e-08
O2
28.106
-3.680e-06
1.746e-05
-1.065e-08
CO
30.869
-1.285e-02
2.789e-05
-1.272e-08
CO2
19.795
7.344e-02
-5.602e-05
C 2 H6
5.409
1.781e-01
-6.938e-05
C 2 H4
3.806
1.566e-01
-8.348e-05
H2 O
32.243
1.924e-03
1.056e-05
H2
27.143
9.274e-03
-1.381e-05
N2
31.150
-1.357e-02
2.680e-05
2
CR IP T
A
1.715e-08
8.713e-09
1.755e-08
-3.596e-09
7.645e-09
-1.168e-08
AN US
Component
Cpi = A + B · T + C · T + D · T
3
with Cpi in
J mol·K
and T in K
Table B.4: Dynamic viscosity coefficients. Source: [69]
A
B
C
D
E
-7.759e-07
5.048e-08
-4.310e-11
3.118e-14
-9.810e-18
O2
-1.026e-06
9.263e-08
-8.066e-11
5.113e-14
-1.295e-17
CO
1.384e-07
7.431e-08
-6.300e-11
3.948e-14
-1.032e-17
CO2
-1.802e06
6.599e-08
-3.711e-11
1.586e-14
-3.000e-18
-4.537e-07
3.554e-08
-9.658e-12
0.00
0.00
-6.216e-07
3.970e-08
-1.259e-11
0.00
0.00
-1.072e-06
3.525e-08
3.575e-12
0.00
0.00
1.802e-06
2.717e-08
-1.340e-11
5.850e-15
-1.040e-18
-1.020e-07
7.479e-08
-5.904e-11
3.230e-14
-6.730e-18
C 2 H6 C 2 H4
H2
CE
N2
PT
H2 O
ED
CH4
M
Component
AC
µi = A + B · T + C · T2 + D · T3 + E · T4 with µi in Pa · s and T in K
53
ACCEPTED MANUSCRIPT
Table B.5: Thermal conductivity coefficients. Source: [69]
B
C
D
E
8.150e-03
8.000e-06
3.515e-07
-3.387e-10
1.409e-13
O2
-1.290e-03
1.070e-04
-5.263e-08
2.568e-11
-5.040e-15
CO
-7.800e-04
1.030e-04
-6.759e-08
3.945e-11
-9.470e-15
CO2
-3.882e-03
5.300e-05
7.146e-08
-7.031e-11
1.809e-14
C 2 H6
-9.100e-04
9.000e-06
2.674e-07
-1.656e-10
4.833e-14
C 2 H4
3.250e-03
-1.000e-05
2.710e-07
-1.568e-10
5.086e-14
H2 O
4.600e-04
4.600e-05
5.115e-08
0.00
0.00
H2
6.500e-04
7.670e-04
-6.871e-07
5.065e-10
-1.385e-13
N2
-1.300e-04
1.010e-04
3.361e-11
-7.100e-15
CH4
2
3
-6.065e-08
CR IP T
A
AN US
Component
λi = A + B · T + C · T + D · T + E · T
4
with λi in
W m·K
and T in K
Appendix C. Kinetic model for the oxidative coupling of methane
M
The kinetic model published by Stansh et al., 1997 [43] was used in this work. This model was developed using La2 O3 /CaO as catalyst and is composed of nine catalytic reactions and a single
ED
non-catalytic homogeneous reaction in gas phase, which involve nine chemical species. The kinetic scheme considers the following reactions: CH4 + 2O2 −→ CO2 + 2H2 O
PT
1 2 O2
2CH4 +
−→ C2 H6 + H2 O
CH4 + O2 −→ CO + H2 O + H2 → CO2
CE
CO +
1 2 O2
C2 H 6 +
1 2 O2
−→ C2 H4 + H2 O
C2 H4 + 2O2 −→ 2CO + 2H2 O
AC
C2 H6 −→ C2 H4 + H2
4H0298 = −802.9 kJ/mol CH4
(C.1)
= −177.0 kJ/mol CH4
(C.2)
4H0298
4H0298 = −277.8 kJ/mol CH4 4H0298 4H0298
C2 H4 + 2H2 O −→ 2CO + 4H2
4H0298
CO2 + H2 * ) CO + H2 O
4H0298
CO + H2 O * ) CO2 + H2
= −283.2 kJ/mol CO
(C.4)
= −104.9 kJ/mol C2 H6
(C.5)
4H0298 = −757.6 kJ/mol C2 H4
4H0298
(C.6)
= +137 kJ/mol C2 H6
(C.7)
= +210 kJ/mol C2 H4
(C.8)
4H0298 = −41 kJ/mol CO
54
(C.3)
= +41 kJ/mol CO2
(C.9) (C.10)
ACCEPTED MANUSCRIPT
Some corrections for the indices in reaction rate definitions of each step in the mechanism reported by Stansh et al., 1997 [43] have been reported as incorrect by several researchers [53, 29, 27]. The corrected reactions rates are as follows:
rr ˙j=
n
m
1+
KjCO2
· e−
4H
j ad,CO2 R·T
· PCO2
!2
=⇒
j = 1, 3 − 6
CR IP T
j Ea
kj0 · e− R·T · POj2 · PC j
n2 ∆Had,O E2 a 2 2 k20 · e− R·T · KO2 · e− R·T · PO2 · Pm CH4
AN US
rr ˙2= 2 n2 ∆H2 ∆Had,O ad,CO2 2 1 + KO2 · e− R·T · PO2 + K2CO2 · e− R·T · PCO2 E7 a
7 rr ˙ 7 = k70 · e− R·T · Pm C2 H6
E8 a
n8 8 rr ˙ 8 = k80 · e− R·T · Pm C2 H4 · PH2 O E9 a
M
n9 9 rr ˙ 9 = k90 · e− R·T · Pm CO · PH2 O E10 a
ED
n10 − R·T 10 rr ˙ 10 = k10 · Pm 0 ·e CO2 · PH2
the C in PC is the respective species providing carbon in the reactants. Values for kinetic parameters
AC
CE
PT
kj0 , Eja , KjCO2 , KO2 , 4Hjad,CO2 , 4Had,O2 , mj and nj are given in Table C.6.
55
ACCEPTED MANUSCRIPT
Table C.6: Kinetic parameters for oxidative coupling of methane. Source: [43]
[kJ/mol]
(B.1)
kj0 mol/g/s/Pa−(m+n)
4Hjad,CO2
48
KjCO2 −1 Pa
2.00e-06
2.50e-13
(B.2)
2.32e+01
182
(B.3)
5.20e-07
(B.4)
4Had,O2
mj
nj
[kJ/mol]
[−]
[−]
-175
0.0
0.0
0.24
0.76
8.30e-14
-186
2.30e-12
-124
1.00
0.40
68
3.60e-14
-187
1.10e-04
104
4.00e-13
-168
(B.5)
1.70e-01
157
4.50e-13
-166
(B.6)
6.00e-02
166
1.60e-13
-211
(B.7)
1.20e+07
226
0.0
0.0
CR IP T
Eja
KO2 −1 Pa
Reaction
a
[kJ/mol]
9.30e+03
300
0.0
0.0
(B.9)
1.90e-04
173
0.0
0.0
(B.10) 2.60e-02 3 −1 mol/s/m /Pa
220
0.0
0.0
a
0.57
0.85
0.0
0.0
1.00
0.55
0.0
0.0
0.95
0.37
0.0
0.0
1.00
0.96
0.0
0.0
0.0
0.0
0.0
0.0
0.97
0.00
0.0
0.0
1.00
1.00
0.0
0.0
1.00
1.00
M
References
0.0
AN US
(B.8)
0.0
[1] J. H. Lunsford, Catalytic conversion of methane to more useful chemicals and fuels: a challenge
ED
for the 21st century, Catalysis Today 63 (2000) 165–174. [2] R. Montenegro, S. Shi Koo Pan, G´ as natural como mat´eria-prima para a produ¸c˜ aode eteno
PT
no estado do Rio de Janeiro, BNDES 12 (2000) 135–60.
CE
[3] H. R. Godini, S. Xiao, S. Jaˇso, S. St¨ unkel, D. Salerno, N. X. Son, S. Song, G. Wozny, Technoeconomic analysis of integrating the methane oxidative coupling and methane reforming pro-
AC
cesses, Fuel Processing Technology 106 (2013) 684 – 694. [4] C. He, F. You, Toward more cost-effective and greener chemicals production from shale gas by integrating with bioethanol dehydration: Novel process design and simulation-based optimization, AIChE Journal 61 (2015) 1209–1232. 56
ACCEPTED MANUSCRIPT
[5] G. Hammer, T. L¨ ubcke, R. Kettner, M. R. Pillarella, H. Recknagel, A. Commichau, H.-J. Neumann, B. Paczynska-Lahme, Natural Gas, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2006, pp. –. URL: http://dx.doi.org/10.1002/14356007.a17_073.pub2.
CR IP T
doi:10.1002/14356007.a17_073.pub2. [6] J. Rostrup-Nielsen, Production of synthesis gas, Catalysis Today 18 (1993) 305 – 324.
[7] K. R. Hall, A new gas to liquids (GTL) or gas to ethylene (GTE) technology, Catalysis Today 106 (2005) 243–246.
AN US
[8] A. Holmen, Direct conversion of methane to fuels and chemicals, Catalysis Today 142 (2009) 2–8.
[9] B. Vora, J. Q. Chen, A. Bozzano, B. Glover, P. Barger, Various routes to methane utilization– SAPO-34 catalysis offers the best option, Catalysis Today 141 (2009) 77 – 83.
M
[10] O. W. Perez-Lopez, T. M. Farias, C. P. Correa, Ethylene production by natural gas partial
ED
oxidation over MgO-based catalysts, in: 4th Mercosur Congress on Process Systems Engineering and 2nd Mercosur Congress on Chemical Engineering, volume 2nd, Universidade Federal
PT
do Rio Grande do Sul, Costa Verde, RJ - Basil, 2005, pp. 1–10. [11] F. Dautzenberg, J. Schlatter, J. Fox, J. Rostrup-Nielsen, L. Christiansen, Catalyst and reactor
CE
requirements for the oxidative coupling of methane, Catalysis Today 13 (1992) 503 – 509. [12] M. Stoukides, Solid-electrolyte membrane reactors: Current experience and future outlook,
AC
Catalysis Reviews 42 (2000) 1–70.
[13] L. Olivier, S. Haag, H. Pennemann, C. Hofmann, C. Mirodatos, A. C. van Veen, High-
57
ACCEPTED MANUSCRIPT
temperature parallel screening of catalysts for the oxidative coupling of methane, Catalysis Today 137 (2008) 80 – 89. [14] L. Mleczko, M. Baerns, Catalytic oxidative coupling of methane—reaction engineering aspects
CR IP T
and process schemes, Fuel Processing Technology 42 (1995) 217 – 248. [15] E. V. Kondratenko, M. Baerns, Oxidative Coupling of Methane, Wiley–VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2008, pp. 3010–3023. URL: http://dx.doi.org/10.1002/ 9783527610044.hetcat0152. doi:10.1002/9783527610044.hetcat0152.
AN US
[16] A. L. Tonkovich, R. W. Carr, R. Aris, Enhanced C2 yields from methane oxidative coupling by means of a separative chemical reactor, Science 262 (1993) 221–223. [17] Y. Jiang, I. V. Yentekakis, C. G. Vayenas, Methane to ethylene with 85 percent yield in a gas recycle electrocatalytic reactor-separator, Science 264 (1994) 1563–1566.
M
[18] A. V. Kruglov, M. C. Bjorklund, R. W. Curr, Optimization of the simulated countercur-
ED
rent moving-bed chromatographic reactor for the oxidative coupling of methane, Chemical Engineering Science 51 (1996) 2945–2950.
PT
[19] N. Lape˜ na-Rey, P. H. Middleton, The selective oxidation of methane to ethane and ethylene in a solid oxide electrolyte reactor, Applied Catalysis A: General 240 (2003) 207–222.
CE
[20] V. R. Choudhary, S. T. Chaudhari, A. M. Rajput, V. H. Rane, Beneficial effect of oxygen distribution on methane conversion and C2 –selectivity in oxidative coupling of methane to C2 –
AC
hydrocarbons over lanthanum-promoted magnesium oxide, J. Chem. Soc., Chem. Commun. 0 (1989) 1526–1527.
58
ACCEPTED MANUSCRIPT
[21] Y. Lu, A. G. Dixon, W. R. Moser, Y. H. Ma, Analysis and optimization of cross-flow reactors with staged feed policies—isothermal operation with parallel–series, irreversible reaction systems, Chemical Engineering Science 52 (1997) 1349 – 1363.
CR IP T
[22] C. Hamel, S. Thomas, K. Sch¨ adlich, A. Seidel-Morgenstern, Theoretical analysis of reactant dosing concepts to perform parallel–series reactions, Chemical Engineering Science 58 (2003) 4483 – 4492.
[23] Y. Lu, A. G. Dixon, W. R. Moser, Y. H. Ma, Oxidative coupling of methane in a modified
AN US
γ–alumina membrane reactor, Chemical Engineering Science 55 (2000) 4901 – 4912.
[24] S. Bhatia, C. Y. Thien, A. R. Mohamed, Oxidative coupling of methane (OCM) in a catalytic membrane reactor and comparison of its performance with other catalytic reactors, Chemical Engineering Journal 148 (2009) 525 – 532.
M
[25] A. L. Y. Tonkovich, R. W. Carr, A simulated countercurrent moving-bed chromatographic reactor for the oxidative coupling of methane: Experimental results, Chemical Engineering
ED
Science 49 (1994) 4647 – 4656.
PT
[26] H. R. Godini, H. Arellano-Garcia, M. Omidkhah, R. Karimzadeh, G. Wozny, Model-based analysis of reactor feeding policies for methane oxidative coupling, Industrial & Engineering
CE
Chemistry Research 49 (2010) 3544–3552. [27] E. Lenhart, E. Esche, H. Arellano-Garcia, L. Biegler, Oxidative coupling of methane: Optimal
AC
operating policies for a membrane reactor network, Chemie Ingenieur Technik 84 (2012) 1980– 1988.
[28] N. Holst, S. Jaˇso, H. R. Godini, S. Gl¨ oser, H. Arellano-Garcia, G. Wozny, J. Steinbach, Two-
59
ACCEPTED MANUSCRIPT
dimensional model for oxidative coupling of methane in a packed-bed membrane reactor, Chemical Engineering & Technology 35 (2012) 294–301. [29] E. Esche, H. Arellano-Garcia, L. Biegler, G. Wozny, Two-dimensional modeling of a packed-bed
CR IP T
membrane reactor for the oxidative coupling of methane, Chemical Engineering Transactions 29 (2012) 1537 – 1542.
[30] E. Esche, H. Arellano-Garcia, L. T. Biegler, Optimal operation of a membrane reactor network, AIChE Journal 60 (2014) 170–180.
AN US
[31] J. Makibar, A. R. Fernandez-Akarregi, L. D´ıaz, G. Lopez, M. Olazar, Pilot scale conical spouted bed pyrolysis reactor: Draft tube selection and hydrodynamic performance, Powder Technology 219 (2012) 49–58.
[32] M. Amutio, G. Lopez, M. Artetxe, G. Elordi, M. Olazar, J. Bilbao, Influence of temperature
M
on biomass pyrolysis in a conical spouted bed reactor, Resources, Conservation and Recycling
ED
59 (2012) 23–31.
[33] G. Alwan, M. Aldahhan, F. Mehdi, Multiobjective Optimization of a Spouted Bed Reactor:
PT
Improving the Stability and Operability, LAP Lambert Academic Publishing, 2012. [34] H. B. Vuthaluru, D. Zhang, Effect of coal blending on particle agglomeration and defluidisation
CE
during spouted-bed combustion of low-rank coals, Fuel Processing Technology 70 (2001) 41 – 51.
AC
[35] B. H. Song, A. P. Watkinson, Three-stage well-mixed reactor model for a pressurized coal gasifier, The Canadian Journal of Chemical Engineering 78 (2000) 143–155.
60
ACCEPTED MANUSCRIPT
[36] M. Olazar, J. M. Arandes, G. Zabala, A. T. Aguayo, J. Bilbao, Design and operation of a catalytic polymerization reactor in a dilute spouted bed regime, Industrial & Engineering Chemistry Research 36 (1997) 1637–1643.
CR IP T
[37] B. D. Hook, H. Littman, M. H. Morgan, Y. Arkun, A priori modelling of an adiabatic spouted bed catalytic reactor, The Canadian Journal of Chemical Engineering 70 (1992) 966–982.
[38] R. S. Krzywanski, N. Epstein, B. D. Bowen, Multi-dimensional model of a spouted bed, The Canadian Journal of Chemical Engineering 70 (1992) 858–872.
AN US
[39] M. Olazar, G. Lopez, H. Altzibar, A. Barona, J. Bilbao, One-dimensional modelling of conical spouted beds, Chemical Engineering and Processing: Process Intensification 48 (2009) 1264– 1269.
[40] A. Niksiar, M. Sohrabi, A novel hydrodynamic model for conical spouted beds based on
M
streamtube modeling approach, Powder Technology 267 (2014) 371 – 380.
ED
[41] P. N. Kechagiopoulos, S. S. Voutetakis, A. A. Lemonidou, Cold flow experimental study and computer simulations of a compact spouted bed reactor, Chemical Engineering and Processing:
PT
Process Intensification 82 (2014) 137 – 149. [42] I. Sanchez, G. Mazza, G. Flamant, D. Gauthier, A streamtube non-isothermal spouted-bed
CE
reactor mathematical model, Chemical Engineering Science 55 (2000) 193–202. [43] Z. Stansch, L. Mleczko, M. Baerns, Comprehensive kinetics of oxidative coupling of methane
AC
over the La2 O3 /CaO catalyst, Industrial & Engineering Chemistry Research 36 (1997) 2568–
2579.
61
ACCEPTED MANUSCRIPT
[44] K. J. Smith, Y. Arkun, H. Littman, Studies on modeling and control of spouted bed reactors—i: Reactor modeling, Chemical Engineering Science 37 (1982) 567–579. [45] L. Biegler, Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical
CR IP T
Processes, MPS-SIAM Series on Optimization, Society for Industrial and Applied Mathematics, 2010.
[46] A. W¨achter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical Programming 106 (2006) 25–57.
AN US
[47] L. T. Biegler, An overview of simultaneous strategies for dynamic optimization, Chemical Engineering and Processing: Process Intensification 46 (2007) 1043 – 1053. [48] M. Olazar, M. J. San-Jos´e, S. Alvarez, A. Morales, J. Bilbao, Measurement of particle velocities
Research 37 (1998) 4520–4527.
M
in conical spouted beds using an optical fiber probe, Industrial & Engineering Chemistry
ED
[49] M. J. San-Jos´e, M. Olazar, S. Alvarez, J. Bilbao, Local bed voidage in conical spouted beds, Industrial & Engineering Chemistry Research 37 (1998) 2553–2558.
PT
[50] A. Niksiar, M. Sohrabi, A. Rahimi, A correction on a one-dimensional model for conical spouted beds published in chemical engineering and processing, 48 (2009) 1264-1269, Chemical
CE
Engineering and Processing: Process Intensification 70 (2013) 289 – 291. [51] L. Mleczko, K.-J. Marschall, Performance of an internally circulating fluidized–bed reactor for
AC
the catalytic oxidative coupling of methane, The Canadian Journal of Chemical Engineering 75 (1997) 610–619.
62
ACCEPTED MANUSCRIPT
[52] H. R. Godini, S. Jaˇso, S. Xiao, H. Arellano-Garcia, M. Omidkhah, G. Wozny, Methane oxidative coupling: Synthesis of membrane reactor networks, Industrial & Engineering Chemistry Research 51 (2012) 7747–7761.
CR IP T
[53] S. Jaˇso, Modeling and Design of the Fluidized Bed Reactor for the Oxidative Coupling of Methane, Ph.d. thesis, Department of Chemical Engineering, Berlin, Germany, 2012.
[54] D. Schweer, L. Meeczko, M. Baerns, Ocm in a fixed-bed reactor: limits and perspectives, Catalysis Today 21 (1994) 357 – 369.
AN US
[55] L. Mleczko, U. Pannek, M. Rothaemel, M. Baerns, Oxidative coupling of methane over a la2o3/cao catalyst. optimization of reaction conditions in a bubbling fluidized-bed reactor, The Canadian Journal of Chemical Engineering 74 (1996) 279–287.
[56] T. P. Tiemersma, Integrated autothermal reactor concepts for oxidative coupling and reforming
M
of methane, Ph.d. thesis, Chemical Engineering Department, Enschede, The Netherlands, 2010.
ED
[57] C. T. Tye, A. R. Mohamed, S. Bhatia, Modeling of catalytic reactor for oxidative coupling of methane using La2 O3 /CaO catalyst, Chemical Engineering Journal 87 (2002) 49–59.
PT
[58] U. Pannek, L. Mleczko, Comprehensive model of oxidative coupling of methane in a fluidizedbed reactor, Chemical Engineering Science 51 (1996) 3575 – 3590.
CE
[59] C. Y. Thien, A. R. Mohamed, S. Bhatia, Process optimization of oxidative coupling of methane for ethylene production using response surface methodology, Journal of Chemical Technology
AC
& Biotechnology 82 (2007) 81–91.
[60] N. A. S. Amin, S. E. Pheng, Influence of process variables and optimization of ethylene yield
63
ACCEPTED MANUSCRIPT
in oxidative coupling of methane over Li/MgO catalyst, Chemical Engineering Journal 116 (2006) 187–195. [61] U. Pannek, L. Mleczko, Reaction engineering simulations of oxidative coupling of methane in
CR IP T
a circulating fluidized-bed reactor, Chemical Engineering & Technology 21 (1998) 811–821. [62] S. Jaˇso, H. Godini, H. Arellano-Garcia, M. Omidkhah, G. Wozny, Analysis of attainable reactor performance for the oxidative methane coupling process, Chemical Engineering Science 65 (2010) 6341 – 6352.
AN US
[63] E. Lenhart, Optimal Design and Operation of a Membrane Reactor Network, Studienarbeit, Fachgebiet Dynamik und Betrieb technischer Anlagen & Carnegie Mellon University – Chemical Engineering, Berlin, Germany, 2010.
[64] N. Epstein, J. Grace, Spouted and Spout-Fluid Beds: Fundamentals and Applications, Cam-
M
bridge University Press, Cambridge, 2011. URL: http://books.google.com.br/books?id=
ED
yMfe4GTWRNsC.
[65] W. Kiatkittipong, T. Tagawa, S. Goto, S. Assabumrungrat, K. Silpasup, P. Praserthdam,
PT
Comparative study of oxidative coupling of methane modeling in various types of reactor, Chemical Engineering Journal 115 (2005) 63–71.
CE
[66] Dippr Diadem, Design institute for physical properties (Dippr) – information and data evaluation manager, 2000.
AC
[67] R. Sinnott, Chemical Engineering Design, volume 6 of Coulson & Richardson’s Chemical Engineering Series, 4 ed., Elsevier Butterworth–Heinemann, Oxford, United Kingdom, 2005.
64
ACCEPTED MANUSCRIPT
[68] R. Reid, J. Prausnitz, B. Poling, The Properties of Gases and Liquids, Chemical engineering series, 4 ed., McGraw–Hill, New York, USA, 1987. [69] V. Verfahrenstechnik, Chemieingenieurwesen, V. Gesellschaft (Eds.), VDI Heat Atlas, 2 ed.,
AC
CE
PT
ED
M
AN US
CR IP T
Springer, Rodenbach, Germany, 2010.
65
Large-Scale DAE-Constrained Optimization Applied to a Modified Spouted Bed Reactor for Ethylene Production from Methane D.M. Yancy-Caballero, L.T. Biegler, R. Guirardello PII: DOI: Reference:
S0098-1354(18)30194-7 10.1016/j.compchemeng.2018.03.017 CACE 6056
To appear in:
Computers and Chemical Engineering
Received date: Revised date: Accepted date:
13 October 2017 16 March 2018 18 March 2018
Please cite this article as: D.M. Yancy-Caballero, L.T. Biegler, R. Guirardello, Large-Scale DAEConstrained Optimization Applied to a Modified Spouted Bed Reactor for Ethylene Production from Methane, Computers and Chemical Engineering (2018), doi: 10.1016/j.compchemeng.2018.03.017
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.
ACCEPTED MANUSCRIPT
Highlights • A spouted bed reactor was proposed to enhance the performance of the oxidative coupling of methane.
CR IP T
• An optimization strategy was developed to demonstrate the theoretical feasibility of the proposed reactor.
• The use a porous membrane as a draft tube in a spouted bed improves the heat and mass
AC
CE
PT
ED
M
AN US
transfer.
1
ACCEPTED MANUSCRIPT
Large-Scale DAE-Constrained Optimization Applied to a Modified Spouted Bed Reactor for Ethylene Production from Methane D. M. Yancy-Caballeroa , L. T. Bieglerc , R. Guirardellob,∗ a Brazilian
CR IP T
Bioethanol Science and Technology Laboratory (CTBE), Brazilian Center for Research in Energy and Materials (CNPEM), R. Giuseppe M´ aximo Scolfaro, 10000, 13083-970, Campinas-SP, Brazil. b School of Chemical Engineering, University of Campinas, Av. Albert Einstein 500, 13083-852, Campinas-SP, Brazil c Department of Chemical Engineering, Carnegie Mellon University, Doherty Hall, 5000 Forbes Avenue Pittsburgh-PA, United States
AN US
Abstract
In this paper, a modified spouted bed reactor is proposed to enhance the yield of the oxidative coupling of methane (OCM). Optimization techniques are used to carry out a theoretical analysis of ethylene production via OCM and define some optimal operating conditions of the reacting system. A model-based DAE-constrained optimization strategy is proposed and applied to the OCM process
M
to illustrate the computational capability of the proposed formulation, and the theoretical feasibility of the proposed reactor. The model developed for the reactor is a one-dimensional model composed
ED
of material, energy, and momentum balances. This model along with the kinetic model constitute a non-linear and differential-algebraic system, which is discretized using orthogonal collocation on
PT
finite elements with continuous profiles approximated by Lagrange polynomials. The resulting
CE
algebraic collocation equations are written as equality constraints in the optimization problem, which is solved with the IPOPT solver within the optimization-modeling platform. An initialization
AC
routine based on simulations was carried out to guarantee convergence in optimizations. Results from simulations and optimizations showed the potential of combining different reactor concepts to ∗ Corresponding
author: Tel: +55 19 3521 3955; fax: +55 19 3521 3965 Email address: [email protected]; P.O. Box 6066. (R. Guirardello)
Preprint submitted to Chemical Engineering Journal
March 19, 2018
ACCEPTED MANUSCRIPT
improve the ethylene production from natural gas via oxidative coupling of methane. Keywords: Dynamic optimization, oxidative coupling of methane, spouted bed reactor, modeling, simulation, ethylene, methane
CR IP T
1. Introduction Natural gas is the third most important energy source in the world nowadays, but only a portion of the extracted natural gas is transported for use since most of it is burned or re-injected into the reservoir [1, 2]. The use of natural gas can be improved if local conversion techniques are developed,
AN US
considering that the natural gas price depends mainly on the transportation. Therefore, conversion processes of methane, which is the principal component of natural gas, are encouraged for natural gas valorization. Several studies have been published on the efficient use of natural gas in chemical and petrochemical industry, taking advantage in this manner of the enormous potential of the world
M
reserves of this natural source of hydrocarbons. Efforts have focused especially in process design problems [3, 4].
ED
The amount of methane used as raw material in the chemical industry represents between 57% of the total consumption [5]. In the long term, with the rapid depletion of crude oil and the
PT
environmental impacts of the oil-derived compounds, the natural gas can become the primary source of energy in the 21st century.
CE
Direct and indirect conversion methods are well known for the conversion of the methane to high value-added chemical products. The indirect use of methane is a more attractive way compared to
AC
the direct one since higher yields are achieved, being these the preferred processes for commercial applications of producing fuels in large-scale [6]. However, these methods are complex processes that require an intermediate step for producing syngas, which consumes enormous amounts of energy causing high production costs [7]. In this sense, direct conversion methods are promising because 3
ACCEPTED MANUSCRIPT
the methane is directly converted without requiring the syngas production step. Partial oxidation and oxidative coupling of methane (OCM) are the principal routes of direct conversion, where methanol, formaldehyde or ethylene can be produced from methane in just one step, but these
CR IP T
processes are challenging because the produced compounds are more active than methane, causing reduced yields [8, 9]. Among the direct conversion methods, the oxidative coupling of methane has received particular attention due to a higher yield than other direct methods. Furthermore, this can be carried out at lower temperatures than the conventional processes for ethylene production [10], and because of it is an exothermic process, the produced energy may be used for purposes of
AN US
steam production.
Two obstacles have hampered the industrial application of the OCM process: the reduced yield of the reaction and the large amounts of energy released [11, 12, 13]. From a practical standpoint, the first limitation can be addressed with the development of more selective catalysts and the
M
advances in the chemical reaction engineering. However, the second limitation, relate to the large heat released, can only be resolved by improving and developing new designs of chemical reactors
ED
since the thermodynamics and fluid dynamics also influence this restriction, and not only the reaction kinetics [14, 15].
PT
Several reactor types and configurations have been proposed for the OCM process, some of the most studied are fixed bed reactors, fluidized bed reactors, membrane reactors, countercurrent
CE
moving bed chromatographic reactors, and electrochemical reactors. All of these types of reactors have its advantages and disadvantages, but only the first three ones have been explored and
AC
developed in detail [14]. Experimental and simulation works have been reported in the literature trying to understand the phenomenological behavior of these reactor types [16, 17, 14, 18, 19, 15]. A distributed oxygen feed along the reactor is just one of the suggested proposals to increase the
4
ACCEPTED MANUSCRIPT
yield of the OCM [20]. In this manner, the application of a membrane reactor to carry out the OCM was proposed to allow flows of small amounts of oxygen into the reaction zone to activate the methane and enable
CR IP T
the coupling of methyl radicals that are formed, limiting at the same time undesired complete combustion reactions. Different membrane types have been used for this purpose (porous membranes, dense membranes, catalytic membranes and perovskite-type membranes). Theoretical studies have reported yields higher than 30% in conventional membrane reactors [21, 22], while experimental ones have reported yields of 28% in porous membrane reactors [23] and around 35% in catalytic
AN US
membrane reactors [24]. Another innovative type of reactor proposed for the OCM is the countercurrent moving bed chromatographic reactor. The highest yields achieved in the oxidative coupling of methane have been reported in this kind of reactor (>50 %) [16, 25]. However, this type of reactor still has to be developed in detail before any large-scale application.
M
Fluidized bed reactors have been considered by some researchers as the best type of reactors for the OCM due to the huge heat released during the reaction, which is more complicated to be
ED
controlled in other types of reactors like the fixed bed reactors [11, 14]. Furthermore, the ability of these reactors to operate isothermally, allow for continuous recirculation and the ease of replacing
PT
the deactivated catalyst, become the fluidized bed reactor a potentially attractive reactor design for the OCM, but most of the selective catalysts have shown agglomerations in fluidized bed reactor
CE
causing problems to get a stable fluidization. Most of the studies about OCM in fluidized bed reactors reported yields less than 20%, which limits the applicability of this type of reactor.
AC
Recently, the Wozny’s research group at Berlin Institute of Technology proposed a reactor
network based on membrane and plug flow/fixed-bed reactors for improving the performance of the OCM process. They investigated feeding policies for an optimal operation of the process combining
5
ACCEPTED MANUSCRIPT
several types of reactors [26, 27, 28, 29, 30]. Performance of 23% yield, 54% selectivity, and 43% methane conversion was reported. In this context, the current work proposes a modified spouted bed reactor as a new design
CR IP T
concept for improving the performance of the oxidative coupling of methane. While the use of the spouted bed technique as a chemical reactor is not novelty since this technique has been used in processes involving gas-solid reactions such as pyrolysis of polyolefins, combustion, and gasification of coal and biomass [31, 32, 33, 34, 35, 36], some significant modifications to the conventional configuration are presented in this work aiming to maximize the ethylene yield. Another valuable
AN US
contribution of the present paper is a comprehensive model-based DAE-constrained optimization strategy, which was developed and implemented to find optimal operating and design conditions of the proposed spouted bed reactor for ethylene production from methane. The results achieved in this study demonstrates the theoretical feasibility of employing the spouted bed technique to
M
enhance the performance of the oxidative coupling of methane.
ED
2. Methodology 2.1. Reactor model description
PT
The spouted bed reactor model was developed by applying mass, momentum and energy balances in one-dimension under steady-state. The proposed reactor was divided into three hydro-
CE
dynamic regions [37, 38, 39, 40]: the spout region, where gas and particle velocities are high, and the particle concentration is low. In this region, the particles and the fluid are in co-current flow;
AC
the annular region between the spout region and the column wall, which is characterized by very low velocities, and where the particle concentration is close to that found in a packed bed. In this region, the particles and the fluid are in counter-current flow, and also there is a percolation of the
6
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 1: Different regions of the spouted bed reactor
fluid from the spout into the annulus through the spout-annulus interface; and a fountain region with a much lower particle concentration than in a packed bed.
The spout and the annular region were taken into account in the conical and cylindrical part
M
of the reactor separately since a cross flow of both the fluid and particles between the spout and the annulus exists only in the conical part for the spouted bed with draft tube. In the cylindrical
ED
part, there is not such a cross flow because of the insertion of the internal draft tube, unless that a porous draft tube is being modeled. For the conventional spouted bed, this cross flow exists along
PT
the whole reactor in both the conical and cylindrical part. This division can be seen in Fig. 1 for the spouted bed with draft tube. It can be seen that both the conical part (region I) and the cylindrical
CE
part (region II) were divided into spout and annulus. The fountain region was not considered in the modeling because this region has no significant effects on the reactor performance results for very
AC
fast catalytic reactions due to the low particle concentration in this region [41, 42]. The division for the conventional spouted bed was the same, with the only difference that the internal draft tube was not taken into account in the model, allowing for a cross flow of the fluid and particles in the
7
(b)
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
(c)
(d)
Figure 2: Control volumes for balances. Global mass balance in spout (a) and annulus (b). Momentum balance of gas in spout (c) and annulus (d)
M
cylindrical part as well.
ED
The spouted bed reactor model for all cases investigated in this work was considered adiabatic, and a pseudo-homogeneous approach was adopted for the energy balance, which simplifies the model
PT
considerably. This means that the temperatures for the gas and solid phase were considered to be approximately the same. While a adiabatic reactor was assumed, it should be highlighted that
CE
in some real cases such an assumption might not hold true since the large heat of reaction of the oxidative coupling of methane could cause the increase of the temperature forming sharp peaks,
AC
which must be controlled by using cooling designs. Another assumption in the development of the reactor model was to consider as negligible the
fluid-particle mass transfer resistance, so that component concentrations on the catalyst surface and in the bulk fluid were assumed to be equal. Fig. (2) shows some selected control volumes for 8
ACCEPTED MANUSCRIPT
which the governing equations (Table 1) have been derived. The kinetic model used in this work describes the oxidative coupling of methane with ten reactions including nine chemical species. This reaction model has been reported in detail by Stansh
CR IP T
et al., 1997 [43], and is based on kinetic measurements in a fixed bed micro-catalytic reactor using La2 O3 /CaO as catalyst. The model was obtained over a wide range of reaction conditions, which increases its applicability in simulations (1 < PO2 < 20 kPa, 10 < PCH4 < 95 kPa, 973 < T < 1228 ˙ STC < 250 kg · s · m−3 ). Further details on the kinetic model and the used parameK, 0.76 < mcat /Q ters as well as the thermodynamic data for all compounds included in simulations and optimizations
AC
CE
PT
ED
M
AN US
can be found in Appendix B and Appendix C.
9
ACCEPTED MANUSCRIPT
Conical
Conservation equation Spout
Equation
Global mass balance
Gas phase
1 As
Solid phase
1 As
·
dAs dz dAs dz
Gas phase
us ·
Solid phase
vs ·
dus dz dvs dz
1 us
+
·
dus dz dvs dz
1 vs · dwi,s dz =
+
Us · ρg,s ·
Component mass balance Momentum balance
·
+ −
1 εs
·
dεs dz
1 (1−εs )
·
dP dz
βB ·(us −vs ) ρp ·(1−εs )
−
+g·
g εs
dεs dz
−
ρg,s ρp
·
=
βB ·(us −vs ) ρg,s ·εs
−1
dTs dz
[Us · ρg,s · Cpg,s + Vs · ρp · Cpp,s ] · +
Annulus
Global mass balance
Gas phase
Solid phase Component mass balance Momentum balance
Gas phase
Solid phase
10 Cylindrical
Spout
Global mass balance
M
Energy balance
Gas phase
ED
Solid phase
Component mass balance
PT
Energy balance Annulus
Global mass balance
CE
Component mass balance
AC
Energy balance
P = rr ˙ j,s · (−4Hrj,s )
j ρp ·Cpp,s d[As ·Vs ] · · [T − T ] a s As dz dρg,a dAa dua 4Ur ·ds 1 1 1 Aa · dz + ua · dz + ρg,a · dz = ua ·εa ·(d2c −d2s ) dAa dva 4Vr ·ds 1 1 Aa · dz + va · dz = va ·(1−εa )·(d2c −d2s ) d[Aa ·Ua ·ρg,a ] dw Ua · ρg,a · dzi,a = cr ˙ i,a · MWi + A1a · dz βB ·(ua +va ) g dua 1 dP ua · dz = − ρg,a ·εa · dz − εa − ρg,a ·εa ρg,a βB ·(ua +va ) a va · dv dz = ρp ·(1−εa ) + g · ρp − 1 "
AN US
Energy balance
dρg,s 4·Ur dz = − εs ·us ·ds 4·Vr (1−εs )·vs ·ds
1 ρg,s
+
cr ˙ i,s · MWi
= − ρg,s1·εs ·
=
CR IP T
Table 1: Mass, momentum and energy balance equations in the different regions of the spouted bed
Region
(1) (2) (3) (4) (5) (6) (7) (8)
· [wi,s − wi,a ]
(9) (10) (11)
ncr P
a [Va · ρp · Cpp,a + Ua · ρg,a · Cpg,a ] · dT rr ˙ j,a · (−4Hrj,a ) dz = ρp · (1 − εa )· j # ngr P d[Aa ·Ua ·ρg,a ] Cpg,a · Aa · (Ts − Ta ) + εa · rr ˙ j,a · (−4Hrj,a ) + dz
(12)
j
1 us 1 vs
·
·
Us ·
dρg,s 1 s · dε dz + ρg,s · dz dεs 1 − (1−εs ) · dz = 0 dw ρg,s · dzi,s = cr ˙ i,s · MWi
dus dz dvs dz
+
1 εs
=0
(15)
[Us · ρg,s · Cpg,s + Vs · ρp · Cpp,s ] ·
dTs dz
00
+Udt · Adt · [Tp,a − Tg,s ]
Gas phase
1 Aa
Solid phase
1 Aa
·
dAa dz dAa dz
+
1 ua
·
+ v1a · · dw Ua · ρg,a · dzi,a
dua dz dva dz
+
1 ρg,a
·
dρg,a dz
P = rr ˙ j,s · (−4Hrj,s ) j
=0
(16) (17)
=0
(18)
= cr ˙ i,a · MWi
(19)
[Va · ρp · Cpp,a + Ua · ρg,a · Cpg,a ] · 00
(13) (14)
00
dTa dz
−K·
dT2a dz2
+Udt · Adt · [Tp,a − Tg,s ] + Uw · Aw · [Tp,a − Tw ]
=
P rr ˙ j,s · (−4Hrj,s ) j
(20)
ACCEPTED MANUSCRIPT
2.2. Simulation of conventional spouted bed reactor When solving the model for the conventional spouted bed reactor, some boundary conditions have to be established to get a realistic solution of the DAE system. These boundary conditions are
CR IP T
easy to be defined at the inlet of the spout region since the inlet gas concentrations and temperature are well defined at the entrance of the bed, which gives the fluid conditions at the bottom of the spout. However, the boundary conditions at the bottom of the annulus are not that clear as there is no a feed stream of fluid in this region for conventional spouted beds, which forces the inlet gas velocity in the annulus to be zero. In consequence, the concentrations of components at the base
AN US
of the annulus cannot be assumed the same as in the inlet of the spout region. This particular condition of conventional spouted bed reactors leads to formulating complex boundary conditions to solve the model. In this way, a suitable boundary condition to be assumed in the annulus would be to satisfy an overall energy balance between the spout and the annulus at the exit of the reactor,
M
since the only knowledge of the temperature at the base of the annulus is that this depends on the temperature of solids moving downward in the annulus, but these temperatures are unknown. This
ED
boundary condition at the exit of the reactor can be simplified as Ta (H) = Ts (H). The knowledge a priori of the concentration of components at the base of the annulus region
PT
is also crucial for solving the DAE system that represents a conventional spouted bed reactor. A suitable way to estimate these concentrations would be to determine the residence time at z = 0
CE
that a particular component takes in the transition spout-annulus. This residence time can be calculated by assuming that the reactor behaves as a continuous flow stirred-tank reactor (CSTR)
AC
at an infinitesimally small control volume of the reactor at the base of the annulus as suggested by Smith et al., 1982 [44]. Thus, this residence time for a cylindrical spouted bed with conical base can be calculated as τg0 =
Hm/6·Umf
(see [44] for details). The initial concentration of all components
11
ACCEPTED MANUSCRIPT
in the annulus can be estimated by applying a mass balance for each component on the assumed CSTR as follows: ci,a0 = ci,s0 + τg0 · cr ˙ i,a (Ta0 )
CR IP T
(21)
With this set of equations, the temperature at the base of the annulus Ta0 is iterated until some convergence criterion is reached to satisfy the boundary condition at the exit of the reactor. This criterion is given by the difference between both temperatures in the annulus and the spout at z = H, which have to be minimum. Details on the implementation of this algorithm in GAMS
AN US
are presented below, and a schematic representation is shown in Fig. A.1. 1. It is assumed a temperature at the base of the annulus region Ta0
2. The system of non-linear algebraic equations given by equation 21 is solved using the CONOPT solver. This step is defined in a separate GAMS script file, which is then nested within the
M
main program. Such a nesting procedure in GAMS was done by using the “call ” and “execute” GAMS commands in the main program. These commands execute the “GAMS script
ED
file” in runtime and store the results in a “.gdx ” file, which are loaded back into the main program.
PT
3. With the guessed Ta0 , the calculated ci,a0 and, the fixed boundary conditions at the inlet of the spout Ts0 and ci,s0 , the DAE system composed by mass, momentum and energy balances
CE
in the two regions of the reactor is solved as an initial value problem using the OCFE method. The solver used to address the problem at this step is IPOPT
AC
4. A convergence criterion is checked (to avoid excessively large computation times and convergence issues, initially a difference of 5 K between the temperatures in the annulus and the spout was defined as the convergence criterion). If this criterion is not satisfied, then a
12
ACCEPTED MANUSCRIPT
new Ta0 is assumed (employing the step 5) and the algorithm continues with step 2. If this criterion is satisfied, then the algorithm continues with step 6 5. It can be seen that the steps 1-4 represent a variant of the well-known shooting methods.
CR IP T
These methods convert a boundary-value problem into an equivalent initial-value problem. Thus, the unspecified initial conditions of the system of differential equations are guessed, and the equations are integrated forward as a set of simultaneous initial value problems. At the end, the calculated final values are compared with the boundary conditions and the guessed initial conditions are corrected if necessary. This procedure is repeated until the specified
AN US
terminal values are achieved within a small convergence criterion. A trial-and-error approach can be implemented to guess the initial conditions every time. However, this is tedious and impractical because the trajectory of the right solution is unknown leading to a relatively large CPU-time-consuming. In this work, instead of using a trial-and-error approach, the false
M
position method was used to estimate the new Ta0 . This approach was applied successfully to this boundary-value problem using another separate GAMS script file and nesting it within
ED
the main program. This approach showed good numerical stability 6. Since the convergence criterion adopted in step 4 does not guarantee the equality of temper-
PT
atures at the exit, this step is required to refine the solution by solving the whole problem as a boundary value problem. In this step, the last solution found in step 3 for all of the
CE
considered variables is used as an initial guess of the new boundary value problem. This time, the system of non-linear algebraic equations given by equation 21 and the boundary condition
AC
at the exit of the reactor given by Ta (H) = Ts (H) are included in the GAMS formulation as equality constraints, and the problem is solved again by using the OCFE method, and the IPOPT solver
13
ACCEPTED MANUSCRIPT
2.3. Optimization strategy The objective function of the optimization problem considered in this work was the C2+ yield (YC2+ ), which depends on the methane conversion (XCH4 ) and the C2+ selectivity (SC2+ ). C2+
CR IP T
is referred to C2 H6 and C2 H4 together. The equality constraints of this optimization problem are composed by a highly coupled and non-linear set of DAEs distributed in space, which was discretized by using the method of orthogonal collocation on finite elements, along with the GaussRadau collocation points because of their compatibility with the NLP formulations and numerical stability properties [45]. These collocation points satisfy the continuity of the states profiles across
AN US
element boundaries automatically, and no additional constraints have to be included in the model. The optimization problem was formulated within GAMS optimization-modeling platform, and the IPOPT solver, which is based on interior-point methods [46], was chosen to solve the problem. This solver can handle quite efficiently problems with large size, a large number of equality and
M
inequality constraints and a potentially large number of degrees of freedom [47, 45]. Before activating the optimizer to carry out any optimizations, some simulation studies were
ED
conducted to determine a proper number of finite elements to get a good accuracy in a reasonable computing time, and to use the simulation results as initial guesses of all variables in optimizations.
PT
These previous simulations also provided a first insight of the behavior of the solution, so that useful ranges of all variables were established for the subsequent optimizations. Initial guesses for
CE
variables were also needed for simulations due to the number of variables involved. So, a systematic initialization routine was proposed where the problem was solved element to element until reaching
AC
the total number of finite elements. Next section details such a routine. For the simulation cases, the number of equality constraints and variables were the same (square
systems). This was accomplished in GAMS by employing a dummy objective function, which can
14
ACCEPTED MANUSCRIPT
simply be any constant value. In this way, the solver tries to find a feasible solution to satisfy all of the equality constraints. However, for the optimization cases, it is always necessary to have more variables than equality constraints, introducing some degrees of freedom into the model. Inequality
CR IP T
constraints were also included in the optimizations based on upper and lower bounds, not only for decision variables but all of the variables in the model. For simulations, however, no bounds for variables were required since the solver was much faster in this case.
2.3.1. Initialization algorithm used in simulations
AN US
When no initial guesses were available to carry out a simulation, it was necessary to divide the entire domain into a higher number of finite elements and then solve the problem element to element sequentially. Thus, the model was firstly solved only for the first finite element, and as this has a very small length, the initial conditions or boundary conditions at the inlet were used as initial guesses for all variables defined in the interior collocation points.
M
After the model was successfully solved for the first finite element, the solution obtained for the
ED
variables at the end boundary of this finite element was employed as an initial guess for solving the model in the next finite element, and so on in a loop that increases the number of finite elements every time until a maximum number of finite elements was reached. In this manner, the solution
PT
for the whole problem was obtained as a sequence of solutions for small problems defined in every
CE
single finite element. Fig. 3 shows this procedure schematically, and an algorithm is presented in Fig. A.2.
AC
2.3.2. Optimization algorithm of the spouted bed reactor Optimizations were solved by following a sequence of three stages: first, a simulation of the
model was solved by fixing the decision variables to constant values. Second, an optimization
15
ACCEPTED MANUSCRIPT
Finite element, 1
Finite element, 2
Finite element, 3
CR IP T
Ic
Finite element, 4
Figure 3: Graphical representation of the initialization routine used in simulations for a sequence of four finite
AN US
elements. Ic stands the initial conditions
was solved using the solution obtained previously in the simulation as initial guesses for all of the variables. In this stage, the decision variables were freed (to have some degrees of freedom) to very tight bounds (upper and lower bounds) for all of them; these bounds were established by
M
defining an increasing factor 4 of ±5% the constant values used in the previous simulation. Third,
ED
a series of optimizations were solved sequentially while slowly increasing the interval defined by the upper and lower bounds (defining new bounds very close to the old ones). In every optimization, initial guesses for all variables and multiplier bounds were taken from the solution of the previous
PT
optimization. This process finished when the bounds of all of the decision variables were inactive,
CE
i.e., the variables were no longer at their bounds, or these could not be moved anymore (negative values for some variables were not allowed). The obtained solution at this point was assumed to be
AC
the final optimal solution, and calculations end up. This optimization procedure is shown in Fig. A.3.
16
ACCEPTED MANUSCRIPT
3. Results and Discussion 3.1. Validation of the fluid-dynamic model A single simulation of a spouted bed without considering reaction is presented in this section in
CR IP T
order to evaluate the fluid-dynamic model. For such a purpose, the model developed in this work was adapted to represent a conical spouted bed with the aim of comparing the simulation results with those results reported by Olazar et al., 2009 [39]. This author also developed a one-dimensional model similar to the model presented here. For this reason, a conical spouted bed was chosen as a
AN US
good reference point to prove the accuracy and reliability of the proposed fluid-dynamic model.
The main goal of carrying out this simulation is to guarantee the use of a relatively simple but reliable model capable of predicting the velocity distributions of gas and particles in every region of a spouted bed as well as the local voidage variation. All of this information is fundamental to solve the spouted bed reactor model (mass balance for component and energy balances).
M
The simulated spouted bed can be seen in Fig. 4, and the used equations for the model are those presented in Table 1 for momentum and global mass balances, where only the conical part
ED
without the draft tube was taken into account. Table 2 presents all of the parameters used in the simulation. The boundary conditions are given by: at the inlet, the axial gas velocity is up to twice
PT
the minimum spouting velocity referred to inlet cross-sectional area, while axial solids velocity is
CE
zero. The pressure at the inlet is given by equation 22 [39], and the voidage is assumed to be 1. At the outlet, the pressure is atmospheric. The axial component of both gas and particles velocities are
AC
zero at the spout-annulus interface. Instead, a radial component for these velocities was considered at this point. h γ i−0.11 −0.06 · (Re0 )ms · D0−0.08 · H1.08 P0 =Patm + 1.20 · ρb · g · tan 2
(22)
where ρb is the bulk density, γ is the cone angle, D0 is the diameter of the bed inlet, H is the height 17
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 4: Graphical representation of the simulated conical spouted bed. Source: taken from Olazar et al., 2009 [39]
Table 2: Parameters used in the simulation of the conical spouted bed
Symbol
Value
Units
γ
33
◦
Patm
101325
Pa
dp
3
mm
g
9.8
Particle density
ρp
2360
Gas density
ρg
1.18
Gas viscosity
µ
1.86e-5
m s2 kg m3 kg m3 kg m·s
D0
0.04
m
Stagnant bed height
H
0.2
m
Column base diameter
Di
0.06
m
Annulus voidage
εa
0.45
−
Spout diameter
ds
0.04
m
M
Parameter
Angle of the column
ED
Atmospheric pressure Particle diameter
AC
CE
PT
Gravity
Inlet diameter
18
ACCEPTED MANUSCRIPT
of the stagnant bed, and (Re0 )ms is the Reynolds number of minimum spouting referred to D0 , which is given by [39]: 0.5
(Re0 )ms =0.126 · Ar
·
Di D0
1.68 h γ i−0.57 · tan 2
(23)
CR IP T
Ar is the Archimedes number, and Di is the diameter of the bed bottom.
Fig. 5 shows a comparison between the simulated results obtained in this work and those results reported by Olazar et al., 2009 [39], as well as a comparison with experimental data for some profiles such as the voidage profile and, gas and particle velocity profiles in the spout. The experimental
AN US
data were reported in [48, 49], and also used by Olazar et al., 2009 [39] for comparison purposes. The model was implemented in GAMS using the OCFE method. After the discretization, the final NLP is composed of a square system of 1700 variables and 1700 equality constraints. The solution of this NLP problem was achieved in GAMS using a dummy objective function. The simulation
M
time was 52 CPUs (4GB and 2.2GHz HP laptop) using the IPOPT solver. A minimum of 15 finite elements and three collocations points were required for a satisfactory
ED
convergence and accuracy. Problems of numerical stability were detected when a smaller number of finite elements was used to solve the model. As observed from Fig. 5, a simple one-dimensional
PT
model provides adequate predictions for gas velocity in spout and voidage distributions. However, the model was unable to predict the particle velocities accurately in the spout. Similar results were
CE
reported by Olazar et al., 2009 [39]. When the results obtained here are compared with the results reported by Olazar et al., 2009 [39],
AC
it can be observed a slightly quantitative difference, even although the profiles agree qualitatively. This difference is mainly caused by some errors in the Olazar et al., 2009 [39] formulation, more
specifically the reported equations for momentum balances, for further details see [50]. Also, it can
19
ACCEPTED MANUSCRIPT
30 Experimental Olazar et al. (2009) This Work
0.9
Experimental Olazar et al. (2009) This Work
28 26 24
0.8
22
us
εs
20 0.7
18
14 12
0.5
10 0.4
0
0.04
0.08
0.12
0.16
8
0.2
0
0.04
z [m]
CR IP T
16
0.6
0.08
0.12
0.16
0.2
z [m]
(a)
(b)
3.5
3
AN US
3
2.5
vs
ua
2
1.5
1.5
1
Experimental Olazar et al. (2009) This Work 0
0.04
0.08
0.12
z [m]
0.06
PT
va
0.08
0.04
CE
0.02
0.04
0.08
0
0.04
0.08
0.12
0.16
0.2
(d) 1.4 1.2 1 0.8 0.6 0.4 0.2
Olazar et al. (2009) This Work 0.12
0.16
0.2
z [m]
0
Olazar et al. (2009) This Work 0
0.04
0.08
0.12
0.16
0.2
z [m]
(e)
(f)
AC
0
Olazar et al. (2009) This Work
z [m]
Ur
0.1
0
0
0.2
ED
(c)
0.16
M
0
0.5
Figure 5: Comparison between results obtained in this work, experimental data [48, 49] and simulated results reported by Olazar et al., 2009 [39]. (a) voidage (b) gas velocity in the spout (c) solids velocity in the spout (d) gas velocity in the annulus (e) solids velocity in the annulus (f) gas velocity at spout-annulus interface
20
ACCEPTED MANUSCRIPT
be observed a reasonable agreement between the results obtained in this work and the experimental data, similar as the results reported by Olazar et. al., 2009 [39]. On the other hand, Olazar et al., 2009 [39] stated that the considerable difference between the
CR IP T
predicted data and the experimental ones for the particle velocities in spout is due to the constant spout diameter considered in the simulations. However, any further study was conducted in this work to verify such a claim. Instead, the spout diameter was assumed constant in all simulations and optimizations, since a variable spout diameter would lead to the incorporation of highly nonlinear equations into the model. While this is not a problem for simulations, it could make difficult
AN US
convergence in optimizations.
3.2. Simulation of a spouted bed reactor for the OCM
3.2.1. Comparison between experimental and simulated data for a conventional spouted bed reactor A reference case was selected in order to check the performance of the proposed model and
M
verify its predictive ability through of a comparison of simulated results and experimental data. In
ED
this fashion, data of a laboratory-scale reactor from Mleczko and Marschall, 1997 [51] were chosen as a base case due to the limited availability in the open literature of data for the OCM in spouted
PT
bed reactors. This was the only published work in considering the OCM in a spouted bed reactor under realistic operating conditions. Mleczko and Marschall, 1997 [51] investigated experimentally
CE
the oxidative coupling of methane in various types of reactors including an internally circulating fluidized-bed, a bubbling fluidized-bed, and a conventional spouted bed.
AC
Table 3 contains the specifications and conditions used in this simulation. Values for the feed gas composition and geometric parameters of the spouted bed were taken from Mleczko and Marschall, 1997 [51] while particle specifications and the kinetic model are given by Stansch et al., 1997 [43]. Details of all the species considered in the OCM and included in the simulations, as well as the 21
ACCEPTED MANUSCRIPT
Table 3: Parameters used in the simulation of the OCM in a spouted bed reactor
Symbol
Value
Units
Height of conical part
zc
0.03
m
Atmospheric pressure
Patm
101325
Pa
dp
0.355
mm
g
9.8
Particle density
ρp
3600
m s2 kg m3
Inlet diameter
d0
0.01
m
Column diameter
dc
0.05
m
Stagnant bed height
H
0.035
m
dbc
3 · Di
m
εa
0.42
−
PCH4
57.10
PO2
22.9
PN2
20.0
CH4 /O2
2.5
−
T0,s
1073-1173
K
Particle diameter Gravity
Diameter at the column base Voidage in annulus
Methane to oxygen ratio Feed gas temperature
kPa
AN US
Feed gas partial pressures
CR IP T
Parameter
properties of such components can be found in Appendix B. The considered reactor is operated in
M
adiabatic mode and the superficial velocity of the feed gas is U0 = 1.1 · Ums . The equations of the
ED
model were solved by following the algorithm presented in section 2.2 since the initial conditions for temperature and composition at the base of the annulus are unknown, so the problem has to be solved as a boundary value problem that satisfies Ta (H) = Ts (H) at the exit of the reactor.
PT
Fig. 6 shows a comparison of the experimental methane conversion, C2+ product selectivity and
CE
C2+ product yield with the results obtained in the simulations for different temperature conditions. By analyzing Fig. 6, it is noted that the proposed model for a conventional spouted bed reactor,
AC
in which the OCM process takes place, presents a reasonable agreement for the three quantities analyzed when compared with experimental data throughout the considered temperature range, even although a slight overprediction is observed in the simulated results. Nevertheless, these results suggest that the model proposed in this work is representative of the spouted bed reactor
22
ACCEPTED MANUSCRIPT
XCH4
SC2+
YC2+
40
20
10
0 1093
1113
1133
1153
1173
AN US
1073
CR IP T
[%]
30
Temperature [K]
Figure 6: Comparison between simulated (lines) and experimental data (bars) for CH4 conversion, C2+ selectivity and C2+ yield. Experimental data taken from [51]
M
behavior.
ED
3.2.2. Simulation of a spouted bed reactor with a draft tube The behavior of a spouted bed reactor with a draft tube for the OCM is also investigated in this work. Thus, a reference case at a fixed feed gas temperature of 1113 K was defined to forecast
PT
how such a modification in a conventional spouted bed affects the performance of the OCM. The
CE
remaining specifications of this reference case were the same as those listed in Table 3. Fig. 7 shows theoretical predictions of the model in both temperatures in the annulus and
AC
spout along the axis for a conventional spouted bed reactor as well as for a spouted bed with draft tube. Temperature profiles present similar tendencies in both reactors. It is appreciated that the temperature in the spout has a mild increase close to the reactor inlet, follows by a remarkable reduction. This temperature gradient occurs in a very short reactor length of just a few millimeters.
23
ACCEPTED MANUSCRIPT
1150
1050
1000
950
CR IP T
Temperature [K]
1100
Ts Ta
900
0
0.5
1
1.5
2
2.5
3
·10−2
AN US
Reactor Length [m]
Figure 7: Temperature profiles in the spout and annulus for conventional spouted bed (solid lines) and spouted bed with draft tube (dashed lines)
The first increase is caused mainly by the heat released by the exothermic reactions that take place
M
in the OCM process while the abrupt reduction is because of the cross-flow of particles coming from
ED
the annulus, which are at a lower temperature. For the remaining reactor length, the temperature profile in the spout follows a smooth increase until it reaches a flat behavior. This flat profile at
PT
the last part of the reactor length for the temperatures in the annulus and spout is caused by the efficient mixing of the solids and gas in a spouted bed operation. It is also observed that the
CE
temperature in the annulus is initially lower than in the spout but, then this starts to rise and after a certain point presents a similar behavior to the temperature in the spout until both temperatures
AC
achieved the same value. This condition was imposed as a boundary condition. Another interesting point to mention is the moderate temperature hot spot occurred just after
the reactor inlet, which is much lower than the hot spots reported in others types of reactors. For example, unacceptable hot spots formation with temperature variation between 400 K and 700 24
ACCEPTED MANUSCRIPT
K have been reported in very small FBR [52, 29, 53, 54, 14] and even of 280 K close to the gas distributor in some fluidized bed reactors [55]. These elevated temperature gradients are regularly avoided by using a large nitrogen dilution in the feed gas or operating the reactor in cooling mode.
CR IP T
In this reference case, a nitrogen dilution of 20% was used, which is much lower than those dilutions commonly found in the literature to prevent that the temperature gets out of control (usually a feed dilution with nitrogen of more than 80% is used) [29, 52].
This numerical simulation demonstrates the theoretical feasibility of a spouted bed used as a chemical reactor for the OCM to control the hot spot formation in an acceptable range, operating
AN US
the reactor in adiabatic mode and using small nitrogen quantities in the feed so that complicated and expensive reactor design for cooling purposes could be avoided.
The differences between the temperature profiles in both types of reactors are a consequence of the performance achieved in each reactor, which was higher in the conventional spouted bed
M
reactor. The higher performance accomplished in the conventional spouted bed reactor can better be appreciated in the Fig. 8, which illustrates the concentration distribution profiles in the spout
ED
and annulus of all the considered species for the two types of reactors. Profiles in both reactors have similar tendencies. However, it is evident that these profiles are
PT
lower in the spouted bed reactor with a draft tube, indicating a lower performance of the OCM in this reactor. This lower performance is caused because the insertion of a device in the spout-
CE
annulus interface reduces the fluid flow and mass transfer in the annular region considerably. In consequence, the C2+ production rates and the methane consumption are also reduced, so that
AC
it can be concluded that the insertion of the draft tube in a spouted bed reactor, in which the OCM takes place, leads to mass transport limitations between the spout and annulus causing a considerable reduction of C2+ yields.
25
ACCEPTED MANUSCRIPT
6
7
6
5
5
4
4 3 3 2 2 1
1
0
0.5
1
1.5
2
2.5
Reactor Length [m]
0
3
(a)
0
0.5
1
1.5
2
Reactor Length [m]
·10−2
AN US
0
CR IP T
m3
mol Molar Concentration
CH4 O2 CO CO2 C2 H6 C2 H4 H2 O H2
2.5
3
·10−2
(b)
Figure 8: Concentration profiles in the spout (a) and annulus (b) region for conventional spouted bed (solid lines) and spouted bed with draft tube (dashed lines)
M
An important aspect to be noted from Fig. 8 is that better C2+ yields (higher concentrations of both ethane and ethylene) were achieved in the annular part of the conventional spouted bed.
ED
This trend can be explained by the distributed feed of the components along the spout-annulus interface caused by the cross-flow of fluid from the spout into the annulus. This means that oxygen
PT
concentrations are uniformly maintained in low levels at every point along the reactor in the annular part. The fact of keeping a distributed feed of oxygen in this region of the reactor favors the reaction
CE
conditions of the OCM, which has previously been demonstrated in some packed bed membrane reactors [56, 29]. Also, it can be seen that concentration distributions of the components remained
AC
almost constant in the annular region after the first half of the reactor because the oxygen is almost completely consumed in this part of the reactor.
26
ACCEPTED MANUSCRIPT
3.2.3. Effects of feed gas temperature and methane-to-oxygen ratio in the feed on the yield A very basic sensitivity analysis is presented in this section to examine the effects of some operating conditions on the methane conversion, C2+ selectivity, and C2+ yield, and thereby, to identify
CR IP T
some key parameters for further improvement of the reactor performance through optimizations of these parameters. To this end, several numerical simulations of the reference case were performed by varying one variable at a time to provide a first insight on the influence these parameters on the reactor performance, and on the optimal regions for them, in which the reactor reaches its highest performance. This does not represent a proper sensitivity analysis as more sophisticated methods
AN US
such as Monte Carlo simulations are needed, neither does not represent an adequate optimization, but help us to get some initial insights of the effect of some parameters on the reaction performance. For this analysis, the feed gas temperature and the methane-to-oxygen ratio were varied but keeping an undiluted feed and the other variables constant. The temperature was varied from 973 K
M
to 1173 K due to the range of applicability of the kinetic model, and the methane-to-oxygen from 2 to 10, respectively. The reason for using a minimum value of 2 for the methane-to-oxygen ratio
ED
is because some preliminary simulations indicated a total oxidation of the hydrocarbons to carbon dioxides with lower methane-to-oxygen ratios. A similar conclusion was reported by Tye et al., 2002
PT
[57]. For all of these reaction conditions, an oxygen consumption almost completely was achieved at the exit of the reactor. The effect of the feed gas temperature and the methane-to-oxygen ratio
CE
on the methane conversion, C2+ selectivity, and C2+ yield is displayed in Figures 9, 10 and 11. A series of several plot types are presented in each case to allow for a better understanding
AC
of the influence of these variables on the specified quantity. First, it is shown the classical 2-D representation plot for each independent variable (temperature and methane-to-oxygen ratio in the feed) to see the effect of each of these variables separately, and then a 3-D surface plot and a 2-D
27
ACCEPTED MANUSCRIPT
XCH4 [%]
30 25 20 15
=2 = 2.5 =3 =4 =5 =6 =7 =8 =9 = 10
40
T = 973 [K] T = 993 [K] T = 1013 [K] T = 1033 [K] T = 1053 [K] T = 1073 [K] T = 1093 [K] T = 1113 [K] T = 1133 [K] T = 1153 [K] T = 1173 [K]
35 30 25 20 15
10
10
5 973
5 1003
1033 1063 1093 1123 Temperature [K]
1153
2
(a)
4
5 6 7 CH4 /O2 ratio
8
9
9
35
8
30
AN US
7
25 15
25
6 5
20
4
5 3
4 5 6 CH4 /O 2
7 8 ratio
9
10
1153 ] 1113 [K re 1073 tu a r 1033 e 993 mp Te
15
3 2
993
1033 1073 1113 Temperature [K]
10 1153
(d)
M
(c)
XCH4 [%]
CH4 /O2 ratio
35
2
10
(b)
10
XCH4 [%]
3
CR IP T
CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2
35
XCH4 [%]
40
ED
Figure 9: Effect of reaction temperature and methane to oxygen ratio in the feed over CH4 conversion. (a) 2-D plot of the effect of T over CH4 conversion at several levels of CH4 /O2 (b) 2-D plot of the effect of CH4 /O2 over CH4 conversion at several levels of T (c) 3-D surface plot of the combined effect of T and CH4 /O2 over CH4 conversion
PT
(d) 2-D contour plot of the combined effect of T and CH4 /O2 over CH4 conversion
CE
contour plot are also given as alternatives to see the combined effect and the interaction between these two independent variables.
AC
Fig. 9 shows the variations of the methane conversion with the feed gas temperature and methane-to-oxygen ratio. More specifically, Fig. 9 (a) presents the effect of the temperature on the methane conversion at different levels of methane-to-oxygen ratio, while Fig. 9 (b) presents the effect of the methane-to-oxygen ratio in the feed on the methane conversion at various levels of 28
ACCEPTED MANUSCRIPT
70
SC2+ [%]
60 50 40
=2 = 2.5 =3 =4 =5 =6 =7 =8 =9 = 10
80
60 50 40
30 20 973
T = 973 [K] T = 993 [K] T = 1013 [K] T = 1033 [K] T = 1053 [K] T = 1073 [K] T = 1093 [K] T = 1113 [K] T = 1133 [K] T = 1153 [K] T = 1173 [K]
70
30 20 1003
1033 1063 1093 1123 Temperature [K]
1153
2
(a)
8
9
10
60
8
50 40 30 20 10 2
3
4
5
C
7
6
9
8
2 /O
t ra
H4
io
50
6
45
5
40 35
4
30
3 2
993
1033 1073 1113 Temperature [K]
25 1153
(d)
M
(c)
10
55
7
AN US
60
SC2+ [%]
CH4 /O2 ratio
70
SC2+ [%]
5 6 7 CH4 /O2 ratio
65
9 80
1033 1073 1113 1153 Temper ature [K ]
4
(b)
10
993
3
CR IP T
CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2
SC2+ [%]
80
ED
Figure 10: Effect of reaction temperature and methane to oxygen ratio in the feed over C2+ selectivity. (a) 2-D plot of the effect of T over C2+ selectivity at several levels of CH4 /O2 (b) 2-D plot of the effect of CH4 /O2 over C2+ selectivity at several levels of T (c) 3-D surface plot of the combined effect of T and CH4 /O2 over C2+ selectivity
PT
(d) 2-D contour plot of the combined effect of T and CH4 /O2 over C2+ selectivity
CE
temperature. As can be noticed from Fig. 9 (a), the methane conversion shows the same trend with increasing temperature for all methane-to-oxygen ratios, passing through a maximum value, which
AC
is much more notable and higher at the lower levels of methane-to-oxygen ratio due to a higher oxygen concentration at these conditions. Fig. 9 (b) shows that the methane conversion decreases almost exponentially when the methane-to-oxygen ratio is increased at a given temperature. It is also noted in Fig. 9 (c-d) that the methane conversion has an optimal point throughout the range of 29
ACCEPTED MANUSCRIPT
YC2+ [%]
13 11 9 7
=2 = 2.5 =3 =4 =5 =6 =7 =8 =9 = 10
15
T = 973 [K] T = 993 [K] T = 1013 [K] T = 1033 [K] T = 1053 [K] T = 1073 [K] T = 1093 [K] T = 1113 [K] T = 1133 [K] T = 1153 [K] T = 1173 [K]
13 11 9 7
5
5
973
1003
1033 1063 1093 1123 Temperature [K]
1153
2
(a)
4
5 6 7 CH4 /O2 ratio
8
9
13
9
12
8
11
AN US
7
11 9 7
10
6
9
5
8 7
4
5 3
4 5 6 CH4 /O 2
7 8 ratio
9
10
1153 ] 1113 [K re 1073 tu a r 1033 e 993 mp Te
993
1033 1073 1113 Temperature [K]
5 1153
(d)
M
(c)
6
3 2
YC2+ [%]
CH4 /O2 ratio
13
2
10
(b)
10
YC2+ [%]
3
CR IP T
CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2 CH4 /O2
YC2+ [%]
15
ED
Figure 11: Effect of reaction temperature and methane to oxygen ratio in the feed over C2+ yield. (a) 2-D plot of the effect of T over C2+ yield at several levels of CH4 /O2 (b) 2-D plot of the effect of CH4 /O2 over C2+ yield at several levels of T (c) 3-D surface plot of the combined effect of T and CH4 /O2 over C2+ yield (d) 2-D contour plot
PT
of the combined effect of T and CH4 /O2 over C2+ yield
CE
temperatures and methane-to-oxygen ratios investigated. This optimal point (methane conversion of 39.15 %) is reached on the minimum methane-to-oxygen ratio used and a temperature of 1106 K.
AC
For C2+ selectivity similar analysis can be done by analyzing Fig. 10, which shows the effect of feed gas temperature and the methane-to-oxygen ratio on the C2+ selectivity. In this case, it is also noticed that the selectivity profile shows a maximum point over the temperature range for all of the methane-to-oxygen ratio levels (see Fig. 10 (a)). This behavior is steeper at the lower levels 30
ACCEPTED MANUSCRIPT
of methane-to-oxygen ratio. Fig. 10 (b) shows that C2+ selectivity increases almost linearly with an increase of the methane-to-oxygen ratio for a given temperature. The effect that the methaneto-oxygen ratio has on the selectivity is opposite to its effect on the conversion, which decreases
CR IP T
with increasing the methane-to-oxygen ratio. Figures 9 (c-d) illustrate an optimal point for C2+ selectivity (70.58 %), which is located at a temperature of 1112 K and on the highest methaneto-oxygen ratio used. However, it is worth noting that the lowest methane conversions were also found at this point, making evident that an intermediate point is required to provide the highest conversions and at the same time a high selectivity.
AN US
Fig. 11 exhibits that C2+ yield has a similar behavior to the methane conversion when varying the feed gas temperature and the methane-to-oxygen ratio but with a much sharper trend, that is, the C2+ yield enhances with increasing temperature until reaching a maximum value and after this point, the bigger the temperature, the lower the C2+ yield. This behavior was the same for all
M
methane-to-oxygen ratios, and again as in the case of the methane conversion, the maximum temperature behavior was much more prominent when decreasing the methane-to-oxygen ratio. This
ED
fact shows that the highest C2+ yields were reached over the points of higher methane conversions but not in those of higher ethylene selectivities. The optimal point achieved for C2+ yield was of
PT
14 % at a feed gas temperature of 1118 K and a methane-to-oxygen ration of 2. The low methane conversions achieved by increasing the methane-to-oxygen ratio is explained
CE
by the kinetic mechanism and stoichiometry (see Appendix C), which suggest that high oxygen concentrations in the feed favor the methane conversion. On the other hand, the rise of the C2+
AC
selectivity with the feed gas temperature before reaching its maximum value results from the bigger activation energy of the selective, primary reaction step (equation C.2) compared to the nonselective ones (equations C.1 and C.2) [14].
31
ACCEPTED MANUSCRIPT
The drop of the C2+ selectivity after its maximum value is caused likely because some unwanted reactions such as oxidation and steam reforming of ethylene (equations C.6 and C.8) are favored at higher temperatures. Moreover, the rise of the C2+ selectivity with increasing the methane-to-
CR IP T
oxygen ratio is because the reaction order with respect to oxygen is lower for the selective, primary reaction step (equation C.2) than for the non-selective ones (equations C.1 and C.2), while the reaction order with respect to methane is higher for the former [58].
All of these aforementioned trends were found in the simulations of the spouted bed reactor and agreed with experimental observations as well as with several simulation work conclusions
AN US
commonly reported in the literature for various types of reactor and catalysts for the OCM [59, 60, 57, 61, 58, 14].
3.3. Optimization of a spouted bed reactor for the OCM
Several operational and geometrical parameters of the spouted bed reactor can be optimized,
M
and this section starts with the optimization including only two decision variables or free variables
ED
that can be adjusted during the optimization process to enhance the performance of the reactor. These variables are the feed gas temperature and the methane-to-oxygen ratio. These two variables
PT
were chosen because a basic sensitivity analysis was done before, varying each of them using simple simulations. So, a first insight is known about where the optimal point can be located whether the
CE
remaining variables are kept constant. The addition of others decision variables to the optimization study is explored in the second part of this section. Finally, the last part presents the optimization
AC
results for a modified spouted bed reactor with a draft tube, in which feed streams of methane and oxygen are split by employing a secondary feed zone of fluid in the annular region of the reactor. In this way, a distributed feed of oxygen is accomplished in the annular zone as a great solution to the low C2+ yield achieved in the conventional configuration of the reactor. 32
ACCEPTED MANUSCRIPT
3.3.1. Optimization of a conventional spouted bed reactor Inequality constraints for decision variables of this first optimization were as follows: 973 ≤ Tg,0 ≤ 1173 and 1 ≤ PCH4/PO2 ≤ 10. In addition to these bounds for the feed gas temperature and the methane-
CR IP T
to-oxygen ratio in the feed, the following bounds were introduced to avoid evaluation errors in some constraints: • All concentrations should be smaller than 20
mol m3
• XCH4 , SC2+ , and YC2+ should be between 0 and 1
AN US
• Partial pressures should be smaller than 106 Pa
• All specific heat capacities should be greater than 10 and smaller than 106 • All mixture averaged molar masses should be greater than 1
J mol·K
kg mol
M
• All dynamic viscosities should be greater than 10−8 and smaller than 1 Pa · s W m·K
ED
• All thermal conductivities should be greater than 10−6
The result of this optimization is shown in Fig. 12 alongside the optimal point found by simply
PT
changing each of these variables at a time as performed in section 3.2.3. It is noted that both results are consistent. A slightly quantitative difference is observed regarding the optimal temperature and
CE
the maximum yield found; this small difference is expected because the simulations were performed in a temperature range of 973 − 1173 K but not every point in this range was included in the
AC
simulations. This small difference can also be caused by the numerical precision of the solver in the different stages in simulations and optimizations since different solver options were adjusted for each case to achieve convergence.
33
ACCEPTED MANUSCRIPT
1118.45
Optimizations results 1121.5
20
14.8
4
5
Temperature
CH4 O2
ratio
1033
2
1003
1
973
0
AN US
C2+ Yield
3
CR IP T
2 2
10
0
1093
1063
[%]
14
5
CH4 /O2 ratio
15
1123
Temperature [K]
Simulations results 25
Figure 12: Comparison between optimal yield found by simulations (blue bars) and optimizations (red bars)
3.3.2. Improving the reactor performance by adding decision variables in optimizations Table 4 contains a list of some key parameters of the conventional spouted bed that may be
M
changed within given bounds. These should not be understood as an exhaustive list of all possible decision variables in a spouted bed. Several other potential decision variables exist. Among those
ED
are the inlet pressures and the type of catalyst. For matters of problem size and complexities in the resulting optimization problem, these were disregarded here. The lower and upper bounds were
PT
chosen after a systematic literature search about the effect of these variables on the OCM process and the spouted bed operation.
CE
For example, bounds for the feed gas temperature and methane-to-oxygen ratio in the feed were selected according to the validity range of the kinetic model [43]. Bounds for dilution with
AC
nitrogen in the feed were defined based on Tiemersma, 2010 [56] and Jaˇso et al., 2010 [62]. Bounds for catalyst properties (density and diameter of catalyst) were chosen according to Esche et al., 2012 [29] and Lenhart, 2010 [63]. Bounds on hydrodynamic parameters were selected taking into
34
ACCEPTED MANUSCRIPT
Table 4: Decision variables for the conventional spouted bed reactor
Parameter
Symbol
Initial value
Lower bound
Upper bound
Units
Operating parameters T0,g
1000
973
1173
[K]
5
1
10
[−]
40
0
80
[%]
U0 Ums
1.1
1
2
[−]
1500
3700
0.3
2
0.05
0.5
[−]
0.07
[m]
CH4 O2
Methane-to-oxygen ratio in the feed
0
N2 |0
Dilution with nitrogen in the feed Feed gas superficial velocity
Design parameters Particle density
ρp
3600
Particle diameter
dp
0.355
Ratio of the fluid inlet diameter to the
d0 dc
0.2
AN US
cylinder diameter
CR IP T
Feed gas temperature
Stagnant bed height
H
0.035
0.02
h
kg m3
i
[mm]
account the operating and geometrical conditions that must be met to achieve a stable spouted bed operation [64].
M
Extra inequality constraints should be considered when some geometrical parameters of a spouted bed device are used as decision variables. While changing these geometrical parameters
ED
during an optimization, some criteria must be satisfied to guarantee a stable operation of spouted
d0 ≤ dc
Umf Ut
1/2
≈
d0 n/2 ≤ εmf dc
(24)
CE
PT
bed. According to Epstein and Grace, 2011 [64], these constraints are as follows:
where d0 is the inlet diameter, dc is the column diameter, Umf is superficial gas velocity at minimum
AC
fluidization, Ut is the terminal gas velocity, and εmf is the voidage at minimum fluidization. n is given in function to the dimensionless Archimedes number as follows [64]: 4.8 − n = 0.047 · Ar0.57 n − 2.4 35
(25)
ACCEPTED MANUSCRIPT
Other two inequality constraints used regarding to the geometrical parameters were:
d0 ≤ 30 dp
where dp is the particle diameter.
CR IP T
dc 30 ≥ n/2 dp εmf
(26)
(27)
AN US
When the above constraints are introduced in the optimization problem, this latter became harder to solve, presenting some instabilities problems and locally infeasible points. The results presented in this section were not obtained by adding all of the decision variables at once to the optimization problem. Instead, these were added to the optimization slowly increasing the degrees of freedom sequentially until all the decision variables were freed. The decision variables were added
M
to the optimization as follows: first, the only two decision variables were the feed gas temperature
ED
and the methane-to-oxygen ratio in the feed as can be seen in previous section. Afterward, the remaining decision variables were added in the following sequence: dilution with nitrogen in the
PT
feed, feed gas superficial velocity, ratio of the fluid inlet diameter to the cylinder diameter, stagnant bed height, particle diameter and particle density. This order choice was arbitrary and, certainly,
CE
the order of adding the decision variables could affect the search for the final optimum, and this can lead to a local optimum, but any further study was conducted in this work to investigate such
AC
an effect on the optimal results. A high number of intermediate steps were performed to achieve the final result for the opti-
mization process, but only five of these intermediate steps will be presented for practical purposes and for giving a general idea of how this process was done until an optimal solution was reached. 36
ACCEPTED MANUSCRIPT
Table 5: Progress of the decision variables in several stages of the optimization process
T0,g LB
CV
UB
LB
1
-
1000
-
2
990
1010UB
3
985
4 5
CV
UB
LB
CV
UB
LB
CV
UB
-
5
-
-
40
-
-
1.1
-
1010
4.0
4LB
6.0
30
50UB
50
1.0
1.2UB
1.2
1040UB
1040
3.0
3LB
7.0
20
52F
975
1080UB
1080
2.0
2.8F
8.0
10
54F
973
1109F
1173
1.0
2.9F
10
0
54F
ρp
No
U0 Ums
N2 |0
d0 dc
dp
CV
UB
LB
CV
UB
LB
CV
1
-
3600
-
-
0.355
-
-
0.2
2
2500
3650UB
3650
0.34
0.5UB
0.5
0.15
3
2000
3660UB
3660
0.33
0.92F
1.0
4
1000
3680UB
3680
0.31
1.0F
5
1500
3700UB
3700
0.3
1.0F
60
1.0
1.2F
1.6
70
1.0
1.38F
1.8
80
1.0
1.38F
2.0
H
UB
LB
CV
UB
-
-
0.035
-
0.3UB
0.3
0.030
0.04UB
0.04
0.1
0.35UB
0.35
0.025
0.05UB
0.05
1.5
0.07
0.35F
0.45
0.220
0.06UB
0.06
2.0
0.05
0.35F
0.5
0.02
0.65F
0.07
AN US
LB
CR IP T
No
CH4 O2 0
M
LB=lower bound; UB=upper bound; CV=current value; F=free variable without bounds active
ED
Whenever a decision variable is at its current lower or upper bound, or at neither of them an index LB, UB or F indicates this respectively. Table 5 shows results for values of decision variables
PT
examined at several stages of the whole optimization procedure until the optimal value was found. Table 6 presents the progress of the methane conversion, C2+ selectivity and C2+ yield at each of
CE
these intermediate steps.
It can be observed that the optimal feed gas temperature and methane-to-oxygen ratio were
AC
slightly different from the optimal values found before, but this differences can be caused due to the interaction of these two variables with others of the included decision variables. Also, there is a trend of using a relatively high amount of dilution with nitrogen, the optimal point was 54%, but
37
ACCEPTED MANUSCRIPT
Table 6: Evolution of the methane conversion, C2+ selectivity and C2+ yield in several stages of the optimization process
C2+ yield [%]
C2+ selectivity [%]
Conversion of CH4 [%]
1
7.38
70.3
10.5
2
10.12
67.0
15.1
3
14.3
64.4
4
18.75
62.1
5
21.86
58.3
CR IP T
No
22.2
30.2 37.5
also some optimizations were carried out by excluding this decision variable, and the yield does not
AN US
drop very much. It was also noted that the optimal inlet superficial velocity is higher than that initially used.
In relation to the catalyst properties, it was found that the optimal particle diameter was in the range of particles classified as Geldart G, which is easier spoutable (1-5 mm). However, an
M
assumption in the model was to disregard fluid-solid mass transfer resistance in the catalyst, and this assumption could cause over-predicted results, which could be more remarkable for larger
ED
particle diameter as those considered in this work. In this way, it would be important to consider an effectiveness factor to account the mass transfer resistance.
PT
More difficult to interpret and less pronounced are effects of catalyst density. Increasing this seems to improve the performance of the reactor, but it is unknown to what extent the catalyst
CE
density can be increased beyond the specified 3600 kg/m3 , since to modify the catalyst density, it is needed to modify the catalyst itself. A value of 3700 was assumed here to be the upper limit,
AC
and at the final stage of the optimization process, the upper bound was still active for this decision variable. Optimal points for geometrical parameters were also found successfully, but the range of these decision variables was not so broad due to the difficulties arise when the additional inequality
38
ACCEPTED MANUSCRIPT
constraints were added to the optimization. It has to be said that between steps four and five (although more intermediate steps were done in between them) the changes made to the variable bounds had to be chosen very carefully,
CR IP T
and the increases consequentially became even smaller to avoid infeasibilities. Finally, in Table 6 can be verified that the optimal C2+ yield found in a conventional spouted bed reactor increase considerably when others decision variables are included in the optimization compared with that previously calculated with only two decision variables.
AN US
3.3.3. Improving the reactor performance by modifying the conventional spouted bed
This section explores the possibility of inserting a porous membrane as a draft tube in the conventional spouted bed configuration to increase the fluid flow rate from the spout region to the annular region. Also, a secondary feed zone of fluid with low inlet velocity in the annular region of the reactor is contemplated. For this configuration, the oxygen is fed into the spout and the
M
methane into the annulus. In this way, it is possible to achieve a distributive feeding of oxygen in
ED
the annulus where the reaction will take place.
Two concepts of reactors will be then combined in order to enhance the performance of the OCM process. The flexibility of a membrane reactor concept will be combined with the excellent
PT
mass and heat transfer properties found in spouted bed reactors. A sketch of this reactor concept
CE
is shown in Fig. 13. It is important to highlight that, although the proposed concept will be theoretically studied here with the aid of optimizations to find useful operational and geometrical
AC
parameters to operate and design the proposed reactor, future work to prove the practical viability of this reactor concept will be necessary. A slight modification needs to be done to the balance equations to account for the modification
included in the spouted bed reactor (membrane as a draft tube). This modification includes the
39
CH4
AN US
O2
CR IP T
ACCEPTED MANUSCRIPT
Figure 13: Spouted-fluid bed with porous draft tube
mass flux in the annular region due to the fluid flow through the membrane tube. Ideally, dense membranes are preferred, but here a porous membrane is considered because the high complexity
M
in modeling dense membranes. The porous membrane modeled here is a γ − Al2 O3 membrane. This particular membrane has been applied in simulation studies of packed bed membrane reactors
ED
[56, 65].
Generally, in a porous membrane all of the components can permeate through the membrane
PT
because this type of membrane is not highly selective, but in this work, it was considered in advance and without any reasonable theoretical fundamental that the only component able to permeate the
CE
membrane is oxygen. In this way, if the pore diameter of these membrane tubes is small enough, the diffusion of oxygen through the membrane can be described by the Dusty-Gas-Model, which
AC
incorporates a contribution of Knudsen diffusion and viscous flow [56].
JO2
1 = R·T
4 K0 3
s
8·R·T P + B0 π · MWO2 µg 40
!
·
4Pi,s−a m (rdt + δm ) · ln rdtr+δ dt
(28)
ACCEPTED MANUSCRIPT
where JO2 is the molar flux of oxygen through the membrane, K0 and B0 are parameters of the model, 4Pi,s−a is the partial pressure difference between the spout and annulus of component i, in this case, oxygen, rdt is the radius of the draft tube, and δm is the membrane thickness.
CR IP T
As it can be seen with this model, the distribution of oxygen through the spout-annulus interface is governed by the partial pressure difference (driving force) and the membrane permeability. Because the exact performance of the membrane is of less interest, and the effect of this on the spouted bed reactor is more important, the above model was considered adequate as a first attempt
AN US
to model a spouted bed reactor with a membrane draft tube, so a constant flux of oxygen was used (equation 28). However, it is necessary to evaluate the performance of the membrane as well as the influence of using a variable oxygen flux to extend the application of this concept. Heat and mass transfer correlations describing transport from the catalyst bed to the membrane wall were taken
M
from literature [56].
Besides the decision variables considered in the optimization case for the conventional spouted
ED
bed reactor, for this optimization, others decision variables were included. These new decision variables were: the draft tube diameter with its lower and upper bounds measured in relation to
PT
the inlet diameter (1.5 · d0 ≤ ddt ≤ 0.5 · d0 ), the distance from the inlet of the draft tube with its lower and upper bounds (0 ≤ h0−dt ≤ H), the height of the draft tube with its lower and upper
CE
bounds (0 ≤ hdt ≤ H), and the superficial velocity of the auxiliary feed in the annulus, which was measured with respect to the minimum fluidization velocity to guarantee a stable spouted bed
AC
operation [64]. The limits of this were 1 ≤ U0,a /Umf ≤ 5. An additional constraint was introduced in this case to avoid that the height draft tube and the
distance from the inlet of the draft tube exceed the stagnant bed height (hdt + h0−dt ≤ H). In this
41
ACCEPTED MANUSCRIPT
Table 7: Results for the optimization of the spouted bed reactor with a porous draft tube Reactor performance
Optimal decision variables
Methane conversion
XCH4
45.4
[%]
Feed gas temperature for both annulus and spout
C2+ selectivity
SC2+
73.2
[%]
Methane-to-oxygen ration in the feed
C2+ yield
YC2+
33.2
[%]
Dilution with nitrogen in the feed in the annulus
Particle density Particle diameter
1119.5
[K]
1.5
[−]
N2 |0
51
[%]
U0 Ums
1.45
ρp
3700
dp
1.2
0
[−] h i kg m3
CR IP T
Feed gas superficial velocity
T0,g
CH4 O2
[mm]
Ratio of the fluid inlet diameter to the cylinder diameter
d0 dc
0.4
[−]
Stagnant bed height
H
0.07
[m]
ddt
1.05 · d0
[m]
h0−dt
0.5
[cm]
hdt
3
[cm]
U0,a Umf
3.5
[−]
Draft tube diameter
Distance from the feed inlet of the draft tube Height of the draft tube
AN US
Superficial velocity of the auxiliary feed in the annulus
optimization study, the feed gas temperature in the annulus and the spout were always the same (T0,g ), but could be different increasing then the degrees of freedom. Another decision variables as membrane properties and feed gas temperature in the annulus can also be included, but this is
M
subject to future work. Table 7 presents the results of the optimization process at the last stage, i.e., when the bounds were inactive for decision variables.
ED
Fig. 14 shows the optimal performance of the three different spouted bed reactors considered, the conventional spouted bed (CSB), the spouted bed with a draft tube (SBDT) and the spouted-
PT
fluid bed with a membrane draft tube (SFBMDT). It is clear that the reactor concept proposed here is superior in performance than the conventional spouted bed reactor optimized in the previous
CE
section, and the spouted bed with draft tube was the reactor with lower performance. The proposal of using a porous membrane as a draft tube seems to be a helpful solution, at least in theory, to
AC
enhance the reactor performance. However, further investigation must be conducted on this reactor concept to prove its validity.
42
ACCEPTED MANUSCRIPT
80
Conversion XCH4 h i Selectivity SC2+ h i Yield YC2+
70 60
40 30 20 10 0 CSB
SBDT
SFBMDT
AN US
Temperature [K]
CR IP T
[%]
50
Figure 14: Comparison between performance of the three spouted bed reactors optimized. CSB - Conventional spouted bed. SBDT - Spouted bed with a draft tube. SFBMDT - spouted-fluid bed with a membrane draft tube
M
4. Conclusions
A simultaneous approach for solving a DAE-constrained optimization problem was used and
ED
applied to a spouted bed reactor for ethylene production from methane to improve the performance of the oxidative coupling of methane. This methodology was able to provide reliable results when
PT
some predicted data were compared with results reported in the literature, verifying the efficiency of the full discretization approach for a DAE-constraint optimization problem. However, it is
CE
important to mention that the model was unable to predict well the solids velocities in the spout, this can be caused by the constant diameter assumed in spout, but some further investigation by
AC
using more rigorous modelling including CFD models must be performed to verify this, since radial variations of velocities or shear stress, which were not taken into account in the modeling, could also affect the solids velocities in the spout.
43
ACCEPTED MANUSCRIPT
On the other hand, the OCFE method implemented for discretization was computationally efficient in terms of CPU time. Therefore, the optimization approach implemented in this work has great potential for design applications and can be successfully applied to theoretical studies of
CR IP T
highly non-linear and complex engineering problems to improve chemical processes. Regarding the non-conventional spouted bed proposed, it can be concluded that the theoretical feasibility of combining several reactor concepts for a spouted bed operation to enhance the performance of the oxidative coupling of methane was demonstrated by means of detailed numerical simulations and optimizations, provided that some modifications of the conventional configuration
AN US
be carried out. Among these modifications, it can be highlighted, the use a porous membrane as a draft tube and the addition of a secondary feed zone of fluid with low inlet velocity in the annular region. These two modifications allow for a distributive feeding of oxygen in the annulus where the reaction takes place. Also, it was demonstrated the feasibility of using a spouted bed as a chemical
M
reactor for the OCM to control the hot spot formation in an acceptable range, operating the reactor in adiabatic mode and using small nitrogen quantities in the feed so that complicated and expensive
PT
Acknowledgments
ED
reactor design for cooling purposes could be avoided.
Financial support of this research was provided by S˜ ao Paulo Research Foundation (FAPESP),
CE
grant numbers 2011/22800-3 and 2013/14896-6, and by National Council for Scientific and Tech-
AC
nological Development (CNPq), grant number 141009/2012-8.
44
ACCEPTED MANUSCRIPT
Nomenclature
Description
BVP
Boundary Value Problem
DAE
Differential–Algebraic Equation
DV
Decision Variables
FBR
Fixed Bed Reactor
GAMS
General Algebraic Modeling System
IPOPT
Interior Point OPTimizer
IVP
Initial Value Problem
LB
Lower Bounds
ncr
number of catalytic reactions
ngr
number of gas reactions
NLP
Nonlinear Programming Problem
OCFE
Orthogonal Collocation on Finite Elements
OCM
Oxidative Coupling of Methane
UB
Upper Bounds
˙ G
PT
˙ Q rr ˙
0
AC
A
00
reaction rate per unit volume in gas reactions mol · m−3 · s−1 reaction rate per mass of catalyst in catalytic reactions mol · g−1 · s−1 mass flow rate of particles kg · s−1 surface area of solids per volume of bed m2 specific surface areas of solids m2 cross–sectional area m2
CE
rr ˙
A
M
component rate mol · m−3 · s−1 mass flow rate of gas phase kg · s−1 volumetric flow rate m3 · s−1
cr ˙
˙ W
ED
Roman Symbols
AN US
Symbol
A
B0 and K0
CR IP T
Acronyms
parameters of the Dusty-Gas- Model
45
ACCEPTED MANUSCRIPT
CD
drag coefficient heat capacity J · mol−1 · K−1 concentration mol · m−3
Cp c
Ea,j F g h J k01,3 k02 K0H2 O K0i KH2 O Ki l
activation energy of reaction j J · mol−1
force acting over differential element [N] acceleration of gravity m · s−2 heat transfer coefficient W · m−2 · K−1 molar flux mol · m−2 · s−1
preexponential factor of rate coefficient, k1,3 mol · bar1/2 · g−1 · h−1 preexponential factor of rate coefficient, k2 mol · bar−1 · g−1 · h−1
preexponential factor of dissociative adsorption constant KH2 O preexponential factor of adsorption constant, KCH4 ,CO,H2 bar−1
dissociative adsorption constant of H2 O
adsorption constants for CH4 , CO and H2 bar−1
length [m]
Pi
partial pressure of component i [Pa]
P
pressure [Pa]
M
m
molecular weight kg · mol−1
MW
ED
mass [kg]
ideal gas constant J · mol−1 · K−1
r
radius [m]
S
PT
R
selectivity in C2 hydrocarbons
T
temperature [K]
X
methane conversion
Y
yield in C2 hydrocarbons
u V
AC
v
CE
w
superficial gas velocity m · s−1 average gas velocity m · s−1 superficial particle velocity m · s−1 average particle velocity m · s−1
U
CR IP T
diameter [m]
AN US
d
mass fraction
46
ACCEPTED MANUSCRIPT
Greek Symbols
δ
differential distance [m]
γ
cone angle [◦ ]
thermal conductivity W · m−1 · K−1
λ
dynamic viscosity [Pa · s] density kg · m−3
µ ρ
residence time [s]
4G0f 4Hads,i 4Hr
4H0f 4z
Gibbs energy of formation J · mol−1 adsorption enthalpy of component i J · mol−1 enthalpy of reaction J · m−3 · s−1 standard enthalpy of formation J · mol−1 height of differential element [m]
ε
voidage or porosity
Subscripts inlet and initial condition
a
annulus
atm
atmospheric
B
Gidaspow model
b
bouyant
bc
base of the column
c
column
cat
catalyst
D
drag
ED
PT
AC
g
CE
dt
M
0
dp
AN US
τ
horizontally projected particle
draft tube gas
g-p
gas–particles
i
components
j
reactions
m
maximum
CR IP T
fluid–particle friction coefficient in drag model kg · m−3 · s−1
β
47
membrane
mf
minimum fluidization
ms
minimum spouting
o
reference state
p
particles
p-g
particles–gas
r
radial
s
spout
s-a
spout–annulus interface
STC
standard temperature condition
w
wall
AC
CE
PT
ED
M
AN US
m
CR IP T
ACCEPTED MANUSCRIPT
48
ACCEPTED MANUSCRIPT
Appendix A. Flowchart of algorithms used in simulations and optimizations Start – Define operating conditions – Define geometric parameters
Guess Ta0
CR IP T
– Calculate thermophysical properties
– Ocfe Method
Solve ci,a0 = ci,s0 + τg0 · cr ˙ i,a
– Conopt Solver
AN US
Set up the initial value pro-
– Gams Script
blem → Ts0 , Ta0 , ci,s0 and ci,a0 – False position method (new Ta0 ) – Gams Script
Solve the Dae system composed of
– Ocfe Method
mass, momentum and energy balances
– Ipopt Solver
Calculate the temperature difference
M
at the reactor exit: Ta (H) − Ts (H)
CE
PT
ED
No
abs (Ta (H) − Ts (H)) ≤ 5 K
Yes
Solution of the Ivp Initial guesses of
Set up the Bvp → Ts0 , ci,s0 as initial
variables for Bvp
conditions and; ci,a0 = ci,s0 + τg0 · cr ˙ i,a and Ta (H) = Ts (H) used as equality constraints
– Ocfe Method
AC
– Gams Script
Solve the Bvp
– Ipopt Solver Solution of the Bvp
Figure A.1: Flowchart of the GAMS algorithm to solve numerically a conventional spouted bed reactor model
49
ACCEPTED MANUSCRIPT
Start Define № of interior collocation points (Icp) Divide the domain into n Fe with n = 1, 2..., j
CR IP T
Select j = 1
Use initial conditions as initial guesses for
Increase n to a higher value
all state variables on the Icp for the jth Fe Solve the problem for j = 1
No
AN US
Optimal solution found?
Yes
Do j = j + 1
Save solution for the (j − 1)th Fe
M
Use solution of the end bound for the
(j − 1)th Fe as initial guess for all state variables on the Icp for the jth Fe
ED
Keep as initial guesses for all state varia-
bles for the (j − 1)th Fe its own solutions
PT
Solve again the problem for the jth Fe
AC
CE
No
End of the program
Optimal solution found?
Yes No
j≤n
Yes
Figure A.2: Flowchart of the algorithm for the initialization routine used in simulations
50
ACCEPTED MANUSCRIPT
Start Define decision variables (Dv)
Activate all Dv Set up maximum and minimum values for upper and lower bounds of the Dv ⇒ Ubmax and Lbmin
Lbmin ≤ Dvk ≤ Ubmax k = 1, 2..., No of Dv
CR IP T
Assign constant values to all Dv
Define realistic bounds for the rest of variables included in the model
Set up tight bounds for all Dv (4 = ±5% of
Run a fully func-
the fixed values used in the previous simulation)
tional simulation
Activate the optimizer by adding a scalar objective function ⇒ YC2+
Valid solution found?
AN US
No
Ubcurrent = constant value of Dv (1 + 4) Lbcurrent = constant value of Dv (1 − 4)
Yes
Set up initial guesses for all
Save current solu-
variables and multipliers
tion for all variables
Activate solver options
Increase Ubmax and Lbmin for
PT
ED
No
CE
Optimal solution found
AC
End of the program
Ipopt.opt file
Run an optimization
M
all Dv violating the condition
to speed up convergence
Ubcurrent ≤ Ubmax Lbcurrent ≥ Lbmin Yes Save current solution for all variables and/or multiplier bounds
No
Check if bounds of all Dv are active Yes Ubcurrent = Ubcurrent + 4 Lbcurrent = Lbcurrent − 4
Figure A.3: Flowchart of the algorithm to carry out optimizations on the spouted bed reactor model
51
ACCEPTED MANUSCRIPT
Appendix B. Thermodynamic properties
Table B.1: Chemical species used in this work. Molecular weights were taken from [66]
Index
Formula
Name
MW [kg/kmol]
1
CH4
Methane
2
O2
Oxygen
3
CO
Carbon monoxide
4
CO2
Carbon dioxide
5
C 2 H6
Ethane
6
C 2 H4
Ethylene
7
H2 O
Water
8
H2
Hydrogen
2.016
9
N2
Nitrogen
28.01
CR IP T
16.04 31.40
28.01
44.01
30.07
28.05
AN US
18.02
Table B.2: Standard enthalpy and Gibbs energy of formation. Source: [67]
CH4
4H0f [kJ/mol]
M
Component
-74.86
-50.49
0.00
0.00
CO
-110.62
-137.15
CO2
-393.77
-394.37
C2 H 6
-84.74
-31.92
C2 H 4
52.33
-68.44
-242.00
-228.59
H2
0.00
0.00
N2
0.00
0.00
PT
ED
O2
H2 O
AC
CE
4G0f [kJ/mol]
52
ACCEPTED MANUSCRIPT
Table B.3: Ideal gas heat capacity coefficients. Source: [68]
B
C
D
CH4
19.251
5.213e-02
1.197e-05
-1.132e-08
O2
28.106
-3.680e-06
1.746e-05
-1.065e-08
CO
30.869
-1.285e-02
2.789e-05
-1.272e-08
CO2
19.795
7.344e-02
-5.602e-05
C 2 H6
5.409
1.781e-01
-6.938e-05
C 2 H4
3.806
1.566e-01
-8.348e-05
H2 O
32.243
1.924e-03
1.056e-05
H2
27.143
9.274e-03
-1.381e-05
N2
31.150
-1.357e-02
2.680e-05
2
CR IP T
A
1.715e-08
8.713e-09
1.755e-08
-3.596e-09
7.645e-09
-1.168e-08
AN US
Component
Cpi = A + B · T + C · T + D · T
3
with Cpi in
J mol·K
and T in K
Table B.4: Dynamic viscosity coefficients. Source: [69]
A
B
C
D
E
-7.759e-07
5.048e-08
-4.310e-11
3.118e-14
-9.810e-18
O2
-1.026e-06
9.263e-08
-8.066e-11
5.113e-14
-1.295e-17
CO
1.384e-07
7.431e-08
-6.300e-11
3.948e-14
-1.032e-17
CO2
-1.802e06
6.599e-08
-3.711e-11
1.586e-14
-3.000e-18
-4.537e-07
3.554e-08
-9.658e-12
0.00
0.00
-6.216e-07
3.970e-08
-1.259e-11
0.00
0.00
-1.072e-06
3.525e-08
3.575e-12
0.00
0.00
1.802e-06
2.717e-08
-1.340e-11
5.850e-15
-1.040e-18
-1.020e-07
7.479e-08
-5.904e-11
3.230e-14
-6.730e-18
C 2 H6 C 2 H4
H2
CE
N2
PT
H2 O
ED
CH4
M
Component
AC
µi = A + B · T + C · T2 + D · T3 + E · T4 with µi in Pa · s and T in K
53
ACCEPTED MANUSCRIPT
Table B.5: Thermal conductivity coefficients. Source: [69]
B
C
D
E
8.150e-03
8.000e-06
3.515e-07
-3.387e-10
1.409e-13
O2
-1.290e-03
1.070e-04
-5.263e-08
2.568e-11
-5.040e-15
CO
-7.800e-04
1.030e-04
-6.759e-08
3.945e-11
-9.470e-15
CO2
-3.882e-03
5.300e-05
7.146e-08
-7.031e-11
1.809e-14
C 2 H6
-9.100e-04
9.000e-06
2.674e-07
-1.656e-10
4.833e-14
C 2 H4
3.250e-03
-1.000e-05
2.710e-07
-1.568e-10
5.086e-14
H2 O
4.600e-04
4.600e-05
5.115e-08
0.00
0.00
H2
6.500e-04
7.670e-04
-6.871e-07
5.065e-10
-1.385e-13
N2
-1.300e-04
1.010e-04
3.361e-11
-7.100e-15
CH4
2
3
-6.065e-08
CR IP T
A
AN US
Component
λi = A + B · T + C · T + D · T + E · T
4
with λi in
W m·K
and T in K
Appendix C. Kinetic model for the oxidative coupling of methane
M
The kinetic model published by Stansh et al., 1997 [43] was used in this work. This model was developed using La2 O3 /CaO as catalyst and is composed of nine catalytic reactions and a single
ED
non-catalytic homogeneous reaction in gas phase, which involve nine chemical species. The kinetic scheme considers the following reactions: CH4 + 2O2 −→ CO2 + 2H2 O
PT
1 2 O2
2CH4 +
−→ C2 H6 + H2 O
CH4 + O2 −→ CO + H2 O + H2 → CO2
CE
CO +
1 2 O2
C2 H 6 +
1 2 O2
−→ C2 H4 + H2 O
C2 H4 + 2O2 −→ 2CO + 2H2 O
AC
C2 H6 −→ C2 H4 + H2
4H0298 = −802.9 kJ/mol CH4
(C.1)
= −177.0 kJ/mol CH4
(C.2)
4H0298
4H0298 = −277.8 kJ/mol CH4 4H0298 4H0298
C2 H4 + 2H2 O −→ 2CO + 4H2
4H0298
CO2 + H2 * ) CO + H2 O
4H0298
CO + H2 O * ) CO2 + H2
= −283.2 kJ/mol CO
(C.4)
= −104.9 kJ/mol C2 H6
(C.5)
4H0298 = −757.6 kJ/mol C2 H4
4H0298
(C.6)
= +137 kJ/mol C2 H6
(C.7)
= +210 kJ/mol C2 H4
(C.8)
4H0298 = −41 kJ/mol CO
54
(C.3)
= +41 kJ/mol CO2
(C.9) (C.10)
ACCEPTED MANUSCRIPT
Some corrections for the indices in reaction rate definitions of each step in the mechanism reported by Stansh et al., 1997 [43] have been reported as incorrect by several researchers [53, 29, 27]. The corrected reactions rates are as follows:
rr ˙j=
n
m
1+
KjCO2
· e−
4H
j ad,CO2 R·T
· PCO2
!2
=⇒
j = 1, 3 − 6
CR IP T
j Ea
kj0 · e− R·T · POj2 · PC j
n2 ∆Had,O E2 a 2 2 k20 · e− R·T · KO2 · e− R·T · PO2 · Pm CH4
AN US
rr ˙2= 2 n2 ∆H2 ∆Had,O ad,CO2 2 1 + KO2 · e− R·T · PO2 + K2CO2 · e− R·T · PCO2 E7 a
7 rr ˙ 7 = k70 · e− R·T · Pm C2 H6
E8 a
n8 8 rr ˙ 8 = k80 · e− R·T · Pm C2 H4 · PH2 O E9 a
M
n9 9 rr ˙ 9 = k90 · e− R·T · Pm CO · PH2 O E10 a
ED
n10 − R·T 10 rr ˙ 10 = k10 · Pm 0 ·e CO2 · PH2
the C in PC is the respective species providing carbon in the reactants. Values for kinetic parameters
AC
CE
PT
kj0 , Eja , KjCO2 , KO2 , 4Hjad,CO2 , 4Had,O2 , mj and nj are given in Table C.6.
55
ACCEPTED MANUSCRIPT
Table C.6: Kinetic parameters for oxidative coupling of methane. Source: [43]
[kJ/mol]
(B.1)
kj0 mol/g/s/Pa−(m+n)
4Hjad,CO2
48
KjCO2 −1 Pa
2.00e-06
2.50e-13
(B.2)
2.32e+01
182
(B.3)
5.20e-07
(B.4)
4Had,O2
mj
nj
[kJ/mol]
[−]
[−]
-175
0.0
0.0
0.24
0.76
8.30e-14
-186
2.30e-12
-124
1.00
0.40
68
3.60e-14
-187
1.10e-04
104
4.00e-13
-168
(B.5)
1.70e-01
157
4.50e-13
-166
(B.6)
6.00e-02
166
1.60e-13
-211
(B.7)
1.20e+07
226
0.0
0.0
CR IP T
Eja
KO2 −1 Pa
Reaction
a
[kJ/mol]
9.30e+03
300
0.0
0.0
(B.9)
1.90e-04
173
0.0
0.0
(B.10) 2.60e-02 3 −1 mol/s/m /Pa
220
0.0
0.0
a
0.57
0.85
0.0
0.0
1.00
0.55
0.0
0.0
0.95
0.37
0.0
0.0
1.00
0.96
0.0
0.0
0.0
0.0
0.0
0.0
0.97
0.00
0.0
0.0
1.00
1.00
0.0
0.0
1.00
1.00
M
References
0.0
AN US
(B.8)
0.0
[1] J. H. Lunsford, Catalytic conversion of methane to more useful chemicals and fuels: a challenge
ED
for the 21st century, Catalysis Today 63 (2000) 165–174. [2] R. Montenegro, S. Shi Koo Pan, G´ as natural como mat´eria-prima para a produ¸c˜ aode eteno
PT
no estado do Rio de Janeiro, BNDES 12 (2000) 135–60.
CE
[3] H. R. Godini, S. Xiao, S. Jaˇso, S. St¨ unkel, D. Salerno, N. X. Son, S. Song, G. Wozny, Technoeconomic analysis of integrating the methane oxidative coupling and methane reforming pro-
AC
cesses, Fuel Processing Technology 106 (2013) 684 – 694. [4] C. He, F. You, Toward more cost-effective and greener chemicals production from shale gas by integrating with bioethanol dehydration: Novel process design and simulation-based optimization, AIChE Journal 61 (2015) 1209–1232. 56
ACCEPTED MANUSCRIPT
[5] G. Hammer, T. L¨ ubcke, R. Kettner, M. R. Pillarella, H. Recknagel, A. Commichau, H.-J. Neumann, B. Paczynska-Lahme, Natural Gas, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2006, pp. –. URL: http://dx.doi.org/10.1002/14356007.a17_073.pub2.
CR IP T
doi:10.1002/14356007.a17_073.pub2. [6] J. Rostrup-Nielsen, Production of synthesis gas, Catalysis Today 18 (1993) 305 – 324.
[7] K. R. Hall, A new gas to liquids (GTL) or gas to ethylene (GTE) technology, Catalysis Today 106 (2005) 243–246.
AN US
[8] A. Holmen, Direct conversion of methane to fuels and chemicals, Catalysis Today 142 (2009) 2–8.
[9] B. Vora, J. Q. Chen, A. Bozzano, B. Glover, P. Barger, Various routes to methane utilization– SAPO-34 catalysis offers the best option, Catalysis Today 141 (2009) 77 – 83.
M
[10] O. W. Perez-Lopez, T. M. Farias, C. P. Correa, Ethylene production by natural gas partial
ED
oxidation over MgO-based catalysts, in: 4th Mercosur Congress on Process Systems Engineering and 2nd Mercosur Congress on Chemical Engineering, volume 2nd, Universidade Federal
PT
do Rio Grande do Sul, Costa Verde, RJ - Basil, 2005, pp. 1–10. [11] F. Dautzenberg, J. Schlatter, J. Fox, J. Rostrup-Nielsen, L. Christiansen, Catalyst and reactor
CE
requirements for the oxidative coupling of methane, Catalysis Today 13 (1992) 503 – 509. [12] M. Stoukides, Solid-electrolyte membrane reactors: Current experience and future outlook,
AC
Catalysis Reviews 42 (2000) 1–70.
[13] L. Olivier, S. Haag, H. Pennemann, C. Hofmann, C. Mirodatos, A. C. van Veen, High-
57
ACCEPTED MANUSCRIPT
temperature parallel screening of catalysts for the oxidative coupling of methane, Catalysis Today 137 (2008) 80 – 89. [14] L. Mleczko, M. Baerns, Catalytic oxidative coupling of methane—reaction engineering aspects
CR IP T
and process schemes, Fuel Processing Technology 42 (1995) 217 – 248. [15] E. V. Kondratenko, M. Baerns, Oxidative Coupling of Methane, Wiley–VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2008, pp. 3010–3023. URL: http://dx.doi.org/10.1002/ 9783527610044.hetcat0152. doi:10.1002/9783527610044.hetcat0152.
AN US
[16] A. L. Tonkovich, R. W. Carr, R. Aris, Enhanced C2 yields from methane oxidative coupling by means of a separative chemical reactor, Science 262 (1993) 221–223. [17] Y. Jiang, I. V. Yentekakis, C. G. Vayenas, Methane to ethylene with 85 percent yield in a gas recycle electrocatalytic reactor-separator, Science 264 (1994) 1563–1566.
M
[18] A. V. Kruglov, M. C. Bjorklund, R. W. Curr, Optimization of the simulated countercur-
ED
rent moving-bed chromatographic reactor for the oxidative coupling of methane, Chemical Engineering Science 51 (1996) 2945–2950.
PT
[19] N. Lape˜ na-Rey, P. H. Middleton, The selective oxidation of methane to ethane and ethylene in a solid oxide electrolyte reactor, Applied Catalysis A: General 240 (2003) 207–222.
CE
[20] V. R. Choudhary, S. T. Chaudhari, A. M. Rajput, V. H. Rane, Beneficial effect of oxygen distribution on methane conversion and C2 –selectivity in oxidative coupling of methane to C2 –
AC
hydrocarbons over lanthanum-promoted magnesium oxide, J. Chem. Soc., Chem. Commun. 0 (1989) 1526–1527.
58
ACCEPTED MANUSCRIPT
[21] Y. Lu, A. G. Dixon, W. R. Moser, Y. H. Ma, Analysis and optimization of cross-flow reactors with staged feed policies—isothermal operation with parallel–series, irreversible reaction systems, Chemical Engineering Science 52 (1997) 1349 – 1363.
CR IP T
[22] C. Hamel, S. Thomas, K. Sch¨ adlich, A. Seidel-Morgenstern, Theoretical analysis of reactant dosing concepts to perform parallel–series reactions, Chemical Engineering Science 58 (2003) 4483 – 4492.
[23] Y. Lu, A. G. Dixon, W. R. Moser, Y. H. Ma, Oxidative coupling of methane in a modified
AN US
γ–alumina membrane reactor, Chemical Engineering Science 55 (2000) 4901 – 4912.
[24] S. Bhatia, C. Y. Thien, A. R. Mohamed, Oxidative coupling of methane (OCM) in a catalytic membrane reactor and comparison of its performance with other catalytic reactors, Chemical Engineering Journal 148 (2009) 525 – 532.
M
[25] A. L. Y. Tonkovich, R. W. Carr, A simulated countercurrent moving-bed chromatographic reactor for the oxidative coupling of methane: Experimental results, Chemical Engineering
ED
Science 49 (1994) 4647 – 4656.
PT
[26] H. R. Godini, H. Arellano-Garcia, M. Omidkhah, R. Karimzadeh, G. Wozny, Model-based analysis of reactor feeding policies for methane oxidative coupling, Industrial & Engineering
CE
Chemistry Research 49 (2010) 3544–3552. [27] E. Lenhart, E. Esche, H. Arellano-Garcia, L. Biegler, Oxidative coupling of methane: Optimal
AC
operating policies for a membrane reactor network, Chemie Ingenieur Technik 84 (2012) 1980– 1988.
[28] N. Holst, S. Jaˇso, H. R. Godini, S. Gl¨ oser, H. Arellano-Garcia, G. Wozny, J. Steinbach, Two-
59
ACCEPTED MANUSCRIPT
dimensional model for oxidative coupling of methane in a packed-bed membrane reactor, Chemical Engineering & Technology 35 (2012) 294–301. [29] E. Esche, H. Arellano-Garcia, L. Biegler, G. Wozny, Two-dimensional modeling of a packed-bed
CR IP T
membrane reactor for the oxidative coupling of methane, Chemical Engineering Transactions 29 (2012) 1537 – 1542.
[30] E. Esche, H. Arellano-Garcia, L. T. Biegler, Optimal operation of a membrane reactor network, AIChE Journal 60 (2014) 170–180.
AN US
[31] J. Makibar, A. R. Fernandez-Akarregi, L. D´ıaz, G. Lopez, M. Olazar, Pilot scale conical spouted bed pyrolysis reactor: Draft tube selection and hydrodynamic performance, Powder Technology 219 (2012) 49–58.
[32] M. Amutio, G. Lopez, M. Artetxe, G. Elordi, M. Olazar, J. Bilbao, Influence of temperature
M
on biomass pyrolysis in a conical spouted bed reactor, Resources, Conservation and Recycling
ED
59 (2012) 23–31.
[33] G. Alwan, M. Aldahhan, F. Mehdi, Multiobjective Optimization of a Spouted Bed Reactor:
PT
Improving the Stability and Operability, LAP Lambert Academic Publishing, 2012. [34] H. B. Vuthaluru, D. Zhang, Effect of coal blending on particle agglomeration and defluidisation
CE
during spouted-bed combustion of low-rank coals, Fuel Processing Technology 70 (2001) 41 – 51.
AC
[35] B. H. Song, A. P. Watkinson, Three-stage well-mixed reactor model for a pressurized coal gasifier, The Canadian Journal of Chemical Engineering 78 (2000) 143–155.
60
ACCEPTED MANUSCRIPT
[36] M. Olazar, J. M. Arandes, G. Zabala, A. T. Aguayo, J. Bilbao, Design and operation of a catalytic polymerization reactor in a dilute spouted bed regime, Industrial & Engineering Chemistry Research 36 (1997) 1637–1643.
CR IP T
[37] B. D. Hook, H. Littman, M. H. Morgan, Y. Arkun, A priori modelling of an adiabatic spouted bed catalytic reactor, The Canadian Journal of Chemical Engineering 70 (1992) 966–982.
[38] R. S. Krzywanski, N. Epstein, B. D. Bowen, Multi-dimensional model of a spouted bed, The Canadian Journal of Chemical Engineering 70 (1992) 858–872.
AN US
[39] M. Olazar, G. Lopez, H. Altzibar, A. Barona, J. Bilbao, One-dimensional modelling of conical spouted beds, Chemical Engineering and Processing: Process Intensification 48 (2009) 1264– 1269.
[40] A. Niksiar, M. Sohrabi, A novel hydrodynamic model for conical spouted beds based on
M
streamtube modeling approach, Powder Technology 267 (2014) 371 – 380.
ED
[41] P. N. Kechagiopoulos, S. S. Voutetakis, A. A. Lemonidou, Cold flow experimental study and computer simulations of a compact spouted bed reactor, Chemical Engineering and Processing:
PT
Process Intensification 82 (2014) 137 – 149. [42] I. Sanchez, G. Mazza, G. Flamant, D. Gauthier, A streamtube non-isothermal spouted-bed
CE
reactor mathematical model, Chemical Engineering Science 55 (2000) 193–202. [43] Z. Stansch, L. Mleczko, M. Baerns, Comprehensive kinetics of oxidative coupling of methane
AC
over the La2 O3 /CaO catalyst, Industrial & Engineering Chemistry Research 36 (1997) 2568–
2579.
61
ACCEPTED MANUSCRIPT
[44] K. J. Smith, Y. Arkun, H. Littman, Studies on modeling and control of spouted bed reactors—i: Reactor modeling, Chemical Engineering Science 37 (1982) 567–579. [45] L. Biegler, Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical
CR IP T
Processes, MPS-SIAM Series on Optimization, Society for Industrial and Applied Mathematics, 2010.
[46] A. W¨achter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical Programming 106 (2006) 25–57.
AN US
[47] L. T. Biegler, An overview of simultaneous strategies for dynamic optimization, Chemical Engineering and Processing: Process Intensification 46 (2007) 1043 – 1053. [48] M. Olazar, M. J. San-Jos´e, S. Alvarez, A. Morales, J. Bilbao, Measurement of particle velocities
Research 37 (1998) 4520–4527.
M
in conical spouted beds using an optical fiber probe, Industrial & Engineering Chemistry
ED
[49] M. J. San-Jos´e, M. Olazar, S. Alvarez, J. Bilbao, Local bed voidage in conical spouted beds, Industrial & Engineering Chemistry Research 37 (1998) 2553–2558.
PT
[50] A. Niksiar, M. Sohrabi, A. Rahimi, A correction on a one-dimensional model for conical spouted beds published in chemical engineering and processing, 48 (2009) 1264-1269, Chemical
CE
Engineering and Processing: Process Intensification 70 (2013) 289 – 291. [51] L. Mleczko, K.-J. Marschall, Performance of an internally circulating fluidized–bed reactor for
AC
the catalytic oxidative coupling of methane, The Canadian Journal of Chemical Engineering 75 (1997) 610–619.
62
ACCEPTED MANUSCRIPT
[52] H. R. Godini, S. Jaˇso, S. Xiao, H. Arellano-Garcia, M. Omidkhah, G. Wozny, Methane oxidative coupling: Synthesis of membrane reactor networks, Industrial & Engineering Chemistry Research 51 (2012) 7747–7761.
CR IP T
[53] S. Jaˇso, Modeling and Design of the Fluidized Bed Reactor for the Oxidative Coupling of Methane, Ph.d. thesis, Department of Chemical Engineering, Berlin, Germany, 2012.
[54] D. Schweer, L. Meeczko, M. Baerns, Ocm in a fixed-bed reactor: limits and perspectives, Catalysis Today 21 (1994) 357 – 369.
AN US
[55] L. Mleczko, U. Pannek, M. Rothaemel, M. Baerns, Oxidative coupling of methane over a la2o3/cao catalyst. optimization of reaction conditions in a bubbling fluidized-bed reactor, The Canadian Journal of Chemical Engineering 74 (1996) 279–287.
[56] T. P. Tiemersma, Integrated autothermal reactor concepts for oxidative coupling and reforming
M
of methane, Ph.d. thesis, Chemical Engineering Department, Enschede, The Netherlands, 2010.
ED
[57] C. T. Tye, A. R. Mohamed, S. Bhatia, Modeling of catalytic reactor for oxidative coupling of methane using La2 O3 /CaO catalyst, Chemical Engineering Journal 87 (2002) 49–59.
PT
[58] U. Pannek, L. Mleczko, Comprehensive model of oxidative coupling of methane in a fluidizedbed reactor, Chemical Engineering Science 51 (1996) 3575 – 3590.
CE
[59] C. Y. Thien, A. R. Mohamed, S. Bhatia, Process optimization of oxidative coupling of methane for ethylene production using response surface methodology, Journal of Chemical Technology
AC
& Biotechnology 82 (2007) 81–91.
[60] N. A. S. Amin, S. E. Pheng, Influence of process variables and optimization of ethylene yield
63
ACCEPTED MANUSCRIPT
in oxidative coupling of methane over Li/MgO catalyst, Chemical Engineering Journal 116 (2006) 187–195. [61] U. Pannek, L. Mleczko, Reaction engineering simulations of oxidative coupling of methane in
CR IP T
a circulating fluidized-bed reactor, Chemical Engineering & Technology 21 (1998) 811–821. [62] S. Jaˇso, H. Godini, H. Arellano-Garcia, M. Omidkhah, G. Wozny, Analysis of attainable reactor performance for the oxidative methane coupling process, Chemical Engineering Science 65 (2010) 6341 – 6352.
AN US
[63] E. Lenhart, Optimal Design and Operation of a Membrane Reactor Network, Studienarbeit, Fachgebiet Dynamik und Betrieb technischer Anlagen & Carnegie Mellon University – Chemical Engineering, Berlin, Germany, 2010.
[64] N. Epstein, J. Grace, Spouted and Spout-Fluid Beds: Fundamentals and Applications, Cam-
M
bridge University Press, Cambridge, 2011. URL: http://books.google.com.br/books?id=
ED
yMfe4GTWRNsC.
[65] W. Kiatkittipong, T. Tagawa, S. Goto, S. Assabumrungrat, K. Silpasup, P. Praserthdam,
PT
Comparative study of oxidative coupling of methane modeling in various types of reactor, Chemical Engineering Journal 115 (2005) 63–71.
CE
[66] Dippr Diadem, Design institute for physical properties (Dippr) – information and data evaluation manager, 2000.
AC
[67] R. Sinnott, Chemical Engineering Design, volume 6 of Coulson & Richardson’s Chemical Engineering Series, 4 ed., Elsevier Butterworth–Heinemann, Oxford, United Kingdom, 2005.
64
ACCEPTED MANUSCRIPT
[68] R. Reid, J. Prausnitz, B. Poling, The Properties of Gases and Liquids, Chemical engineering series, 4 ed., McGraw–Hill, New York, USA, 1987. [69] V. Verfahrenstechnik, Chemieingenieurwesen, V. Gesellschaft (Eds.), VDI Heat Atlas, 2 ed.,
AC
CE
PT
ED
M
AN US
CR IP T
Springer, Rodenbach, Germany, 2010.
65