CFD simulations of Industrial Steam Cracking Reactors: Turbulence

Subscriber access provided by READING UNIV

Article

CFD simulations of industrial steam cracking reactors: turbulence-chemistry interaction and dynamic zoning Pieter A Reyniers, Carl M. Schietekat, Bo Kong, Alberto Passalacqua, Kevin Marcel Van Geem, and Guy B Marin Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b02492 • Publication Date (Web): 28 Nov 2017 Downloaded from http://pubs.acs.org on December 4, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

CFD simulations of industrial steam cracking reactors: turbulence-chemistry interaction and dynamic zoning Pieter A. Reyniers a, Carl M. Schietekat a, Bo Kong b, Alberto Passalacqua b, Kevin M. Van Geem a,*, Guy B. Marin a

a

Ghent University, Laboratory for Chemical Technology, Technologiepark 914, 9052 Gent,

Belgium. b

Department of Chemical and Biological Engineering, Iowa State University, Ames, IA,

United States *

Corresponding author: Technologiepark 914, 9052 Gent, Belgium;

[email protected]

ACS Paragon Plus Environment

1 of 49

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 49

Abstract: Modeling finite-rate chemistry in turbulent reacting flows is challenging because of the large span in length and time scales. Reynolds-averaged Navier-Stokes equations-based simulations do not resolve turbulent fluctuations and hence neglect their effect on the reaction rates. Turbulence-chemistry interaction was accounted for in RANS simulations via quadraturebased integration of the reaction rates, calculated using a temperature probability density function with a presumed Gaussian shape and transported mean and variance. The effect on light olefin yield was 0.1-0.3 wt% absolute. A dynamic zoning method was implemented to reduce the computational cost by performing chemical rate calculations only once for thermodynamically similar cells. Speedups of 50 to 190 were observed while the relative error on conversion remained below 0.05 %. The advantages of the presented methodology were illustrated for a large-scale butane-cracking U-coil reactor.

Keywords: Steam cracking, Computational fluid dynamics, Turbulence-chemistry interaction, Dynamic zoning method, U-coil reactor, Ethene


2 of 49

Page 3 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Modeling finite-rate chemistry in turbulent reacting flows presents a challenge because of the large span in length and time scales. In large eddy simulation (LES), the largest turbulent time scales are resolved, while in direct numerical simulation (DNS), all turbulent time scales are explicitly resolved, implying that local fluctuations and their effect on reaction rates are explicitly accounted for. Although requiring substantial computational resources, LES turbulence models have been successfully applied to reacting flows, particularly in combustion research

1-6

. The present work focuses on steam cracking of hydrocarbons, i.e.

pyrolysis. Although pyrolysis reactions are inherently a part of kinetic models for combustion, the simulation requirements for combustion and steam cracking are quite different. Since often the main objectives of a CFD simulation of a combustion process are the prediction of heat release, flame temperature and flame stability, reduced kinetic models and/or mixture fraction approaches are typically sufficient. The main exception is the simulation of formation during combustion, which requires a detailed kinetic model. Even then formation can be accurately calculated in a post-processing step, decoupling the flow dynamics from the chemistry calculation 7-9. The alternative is to use tabulation methods such as in situ adaptive tabulation to reduce the computational cost of calculating the species rates of production with a large kinetic network 10, which has also been applied outside of the combustion community

11

. Other techniques precompute the species rates of formation to

minimize redundancy in the calculations

12

. Whereas typically combustion proceeds to full

conversion to water and carbon dioxide, steam cracking reactors are operated at a partial conversion to maximize the production of valuable light olefins. Hence, the minimum size of


3 of 49


Page 4 of 49

the kinetic models, necessary for steam cracking reactor simulations, is larger than that required for typical combustion simulations. There is a trade-off between the level of accuracy that can be attained in the description of the flow phenomena on the one hand and the description of the chemical phenomena on the other hand. The largest available kinetic mechanisms, consisting of hundreds of thousands of reactions amongst thousands species, can only be solved in 1D reactor models, i.e. the ideal plug flow reactor. The paramount objective of simulations using ideal reactor models is the accurate simulation of the product yields. However, not all industrial reactors can be simulated via these ideal reactor models

13

because flow phenomena occur that, by virtue of

being 2D or 3D in nature, cannot be simulated with a 1D model but have a significant influence of the simulated product yields. A full 3D description of these reactors is necessary, but in order to obtain results at a reasonable computational cost and within a reasonable timeframe, the kinetic mechanisms describing the chemical phenomena have to be simplified, which reduces the accuracy of the prediction of the end-of-pipe yields compared to using a full kinetic mechanism. Nevertheless, for some enhanced steam cracking reactors, the 3D simulation with a less detailed kinetic mechanisms can provide a more accurate prediction of the product yields compared to a 1D simulation with a detailed kinetic mechanism due to the large effect of the flow phenomena on the chemistry. The Reynolds number in steam cracking reactors ranges from 100,000 to more than 200,000 14

. The grid requirements for wall-resolved LES scale to the power of 1.8 with the Reynolds

number

15, 16

, currently limiting LES to Reynolds numbers of approximately 100,000 for

reasons of computational cost. The combination of the large kinetic model and the large grid size renders a LES of an industrial steam cracking reactor computationally very demanding. Therefore, Reynolds-averaged Navier-Stokes equations are used in this work, which do not


4 of 49

Page 5 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


take turbulence into account explicitly but model its effect on the flow, in line with previous work

17-20

. However, neglect of the turbulent fluctuations by using the mean temperature and

concentrations for the calculation of reaction rates can induce an error:

̅ ≠ ,

(1)

and this mainly due to the highly nonlinear dependence of the reaction rate coefficient on the temperature in the Arrhenius equation: ,

≠ =

(2)

The effect of turbulent fluctuations on the reaction rates can be accounted for by introducing scalar-variance transport equations for temperature and species concentrations. As the temperature dependence (exponential) is stronger than the concentration dependence (first or second order, in pyrolysis kinetic models), only the effect of temperature fluctuations is accounted for in this work when considering turbulence-chemistry interaction. Scalarvariance transport equations have been used for several applications, one of the most important being the simulation of turbulent premixed flames

21-23

or turbulent non-premixed

flames 24, 25. To reduce the computational load in those turbulent premixed flame simulations, a detailed chemistry model based on elementary reactions can be replaced with a method using the concept of the reaction progress variable

26

, which is defined in such a way that it

spans the interval between unburned and fully burned fuel by increasing monotonically from zero to one. Transport equations for the progress variable mean and variance are solved, and the species concentrations are reconstructed in a post-processing step based on a predefined probability density function (PDF) of the reaction progress variable. Scalar-variance transport equations have also been used to predict formation in turbulent combustion 7. In a post-


5 of 49


Page 6 of 49

processing step, the formation rate of is calculated based on the temperature mean and variance and mean species concentrations. In spite of the large computational cost associated with coupling three-dimensional flow models with detailed chemical mechanisms. A trend of utilizing such a strategy has emerged in combustion research

27, 28

but spread fast to other fields

13, 29, 30

. The improved accuracy

comes at a high computational cost, to the point where the chemical rate calculations account for a significant part of the total solution time. Furthermore, three-dimensional reactor technologies, aiming at enhanced heat transfer, induce flow fluctuations and turbulence via modifications of the wall. Hence, an accurate description of the near-wall region using a lowReynolds formulation of the turbulence model is desirable. This requires a first-cell "

dimensionless wall distance = τ ⁄ # % & below unity, yielding a first-cell thickness in the µm range. Given the large length of steam cracking reactors, i.e. from 10 to 100 m, this results in grids containing tens of millions of cells. The chemical rate calculations are carried out in every computational cell at every time step, so there is a certain degree of redundancy as the resulting species production rates for thermodynamically similar cells, i.e. cell with nearly the same temperature and species concentrations, are nearly identical. The dynamic zoning method as introduced by Liang et al.

31

, in which these thermodynamically similar

cells were grouped together, and the species production rate calculations were performed only once for an averaged cell state. The calculated rates are subsequently mapped to the relevant cells, and the species conservation equation is solved using the mapped production rate, thus realizing a considerable speedup due to the lower number of chemical rate calculations. In this work, a code for the three-dimensional simulation of steam cracking reactors, based on the open source CFD package OpenFOAM® 2.2.0 was developed


32

. The effect of turbulent

6 of 49

Page 7 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


temperature fluctuations on the reaction rates is accounted for by solving a temperature variance transport equation. To reduce the computational cost, the code also incorporates a dynamic zoning method. In the following, the governing equations are summarized and the numerical methods used are introduced. Next, the importance of turbulence-chemistry interaction in steam cracking reactors is assessed, the speedup of the dynamic zoning method is investigaged. Finally, a large-scale steam cracking reactor is simulated to demonstrate the effectiveness of this code.

2. Governing equations 2.1. Conservation equations Only the most important aspects of the simulation framework are discussed here, the reader is referred to the Supplementary Material for additional information. Steam cracking reactors are gas-phase flow tubular reactors made of iron-nickel-chromium alloys. To account for the conductive heat transfer in the reactor metal, both the solid phase, i.e. the reactor wall, and the fluid phase, i.e. the process gas, are modeled. Since the process is non-isothermal and non-isobaric, both the momentum and energy equations need to be solved for the fluid phase 33. The ideal gas law can be used for the fluid phase as the compressibility factor of the mixture at high temperatures is low. In the solid region of the reactor, encompassing the metal reactor wall, only the conductive heat transfer equation is solved. The governing steady-state conservation equations are summarized in Table 1. These equations are valid in general and can be combined with any of the available RANS standard turbulence models of OpenFOAM® or new user-implemented turbulence models to account for the effects of turbulence on the flow.


7 of 49


Page 8 of 49

In the Navier-Stokes equations, the laminar viscosity of the ideal gas mixture is calculated from the species laminar viscosity using an ideal mixing law. The species laminar viscosities '(,) are calculated from kinetic theory using the species Lennard-Jones well depth *+, and -+, . Most values for *+, and -+, are taken from literature

RMG’s TransportDataEstimator

35

34

, with the remaining obtained from

. In the energy equation (16),

./,0// 123

is the thermal

diffusivity containing contributions of laminar and turbulent conduction. The laminar conductivity coefficient is calculated from kinetic theory while the turbulent conductivity coefficient is calculated from the turbulent viscosity. The heat of reaction 45 is calculated from the species heat of formation. The species continuity equations (17) are solved for 6 − 1 species, with the concentration of the last species calculated from the global mass

balance (14). The effective diffusion coefficient 9:;;,) in the species continuity equations also has a laminar and a turbulent contribution. 17. To account for the effect of local turbulent fluctuations of temperature on the reaction rates, a temperature probability density function (PDF) is used as explained in section 3.2. The temperature variance is obtained by solving the scalar variance equation (3). The derivation of this equation relies on the eddy diffusivity concept of Kolmogorov’s theory for the unclosed terms 36. The laminar transport contribution in the first term on the right-hand side of equation (3) is neglected because of the high Reynolds number flows considered in this work.

-> = ∇ ∙ @ ∇ =

'A > 'A ∇ > − * ∇- B + 2.85 -A 4HA

(3)

There are two distinct ways of determining the unclosed temperature variance dissipation * . In fully developed flow at high Reynolds numbers, the algebraic equilibrium model, i.e. equation (4) can be used. The parameter & is typically set to 2.0, based on the experimental work of Béguier et al. 37.


8 of 49

Page 9 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


* * = & - >

(4)

In reality, the relation between the temperature variance dissipation * and the characteristics of the turbulent flow field is more subtle. By solving an additional conservation equation (5) for the dissipation of the temperature variance, its value can be determined by rigorously accounting for the occurring physical mechanisms 38, 39.

* = ∇ ∙ @ ∇ =

'A *> ∇* B − > > − I 4HA -

>

J J * ∇ > + L 'A |4| * + K 'A 4HA

(5)

The terms with the model parameters > and I are the dissipation terms due to scalar and mechanical destruction of fluctuations, respectively. The terms with the model parameters K

and L are the production terms due to scalar gradients and velocity gradients, respectively.

Many sets of values have been proposed for the model parameters > , I , K and L . An

extensive overview is given by Sanders and Gokalp

40

. The parameter J is part of the

turbulence model and equal to 0.09 in most k-ε turbulence models. The term |4| is the norm of the strain rate tensor defined as:

4) =

1 O P O P) N + R 2 O Q) O Q

(6)

2.2. Species rate of formation Thermal cracking of hydrocarbons proceeds mainly through a radical mechanism. Three major reaction families must be considered: carbon-carbon and carbon-hydrogen bond scission and the reverse radical-radical recombination reactions, intra- and intermolecular hydrogen abstraction reactions, and the radical addition to olefins and the reverse β scission of radicals, both intra- and intermolecular. The adopted kinetic models in this work are singleevent microkinetic models that only consider elementary reactions. The kinetic models were


9 of 49


Page 10 of 49

derived from the model of Van de Vijver et al. 28 by removing unimportant species identified via a combination of sensitivity analysis PRO

44

41, 42

and rate of production analysis

43

. CHEMKIN

was used to calculate the sensitivity coefficients and the production rates of the

species. In agreement with Brock et al. 45, reactions were included in the skeletal mechanism if the sensitivity coefficients were higher than 0.1 or if the net rate was higher than 5% of the net rate of the fastest step in the considered time interval. The validity of the reduced network is shown in the Supplementary Material by comparing one-dimensional plug flow simulation results obtained with the reduced and the full models. Thermodynamic data of most molecules and radicals was derived from first principles CBS-QB3 calculations of Sabbe et al.

46, 47

.

Missing data was obtained from RMG’s ThermoDataEstimator 35. The average reaction rate coefficients are calculated by integrating a prescribed temperature PDF to account for turbulence-chemistry interaction: V

= S T U W

(7)

with the mean reaction rate coefficient used for the calculation of the reaction rate ̅ = [\ ∏)]& %,) H) Y,Z . A Gaussian distribution is assumed for the temperature PDF:

T =

1

&

2 ^ -> >

# > _`#

(8)

with the mean temperature obtained from the energy equation (16), and -> the temperature variance. To obtain a value for -> in each grid cell, the temperature variance transport

equation (3) is solved. Based on user input, turbulence-chemistry interaction can be turned on/off in the code, i.e. whether equation (2) or (7) is solved in the calculation of the species rate of formation. If turbulence-chemistry interaction is considered, the algebraic closure


10 of 49

Page 11 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model, i.e. equation (4), or one of the fifteen sets of implemented parameter values for equation (5) can be used.

3. Numerical methods 3.1. Conservation equations OpenFOAM® applies the finite-volume method on co-located grids and adopts Gaussian integration. A segregated solver was used to solve the conservation equations given in Table 1. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE)

48

algorithm

was used to solve the Navier-Stokes equations. The mesh orthogonality was high in all simulations, thus no non-orthogonal correction was applied. The second-order central differencing scheme ‘Gauss linear’ was used for the discretization of gradients while for the Laplacian terms ( ∇ ∙ 9 ∇a ) the ‘Gauss linear uncorrected’ scheme was used. The interpolation scheme was ‘linear’ and the default surface normal gradient scheme was ‘corrected’, although the explicit non-orthogonal corrector included in this scheme was not necessary due to the high mesh orthogonality. For the majority of the divergence terms of the form ∇ ∙ … (excluding Laplacian terms), the ‘bounded Gauss upwind’ scheme was used, a first-order bounded scheme. One major exception is the advection term of species mass fraction, for which the ‘Gauss limitedLinear01’ scheme was used, imposing a lower limit of zero and an upper limit of unity on the solution of the transport equations for the species mass fraction. The ‘Gauss limitedLinear01’ scheme is a linear scheme (second order, unbounded) that limits towards upwind (first order bounded) in the presence of a strongly varying gradient. The 01 implies that the bounding in the upwind scheme is done between 0 and 1 by means of the Sweby limiter 49, 50.


11 of 49


Page 12 of 49

For the turbulence model, the shear-stress transport (SST) k-ω model was used to calculate the eddy viscosity 51 as it proved to give reasonably accurate results for similar flows 52. When using detailed kinetic models, the species conservation equations become stiff because of the significant difference in time scales associated with species and reactions. The typical lifetime of a radical species is several orders of magnitude smaller than the lifetime of a non-radical species which results in a large stiffness ratio, which is the ratio of the largest and smallest eigenvalue of the chemical Jacobian matrix. To tackle this issue, the pseudo-steady state assumption (PSSA) is applied to the radicals on-the-fly using a numerical algebraic solver. Besides reducing the stiffness of the species conservation equations, also the number of species conservation equations is reduced as the concentrations of the PSS species can be calculated in the algebraic solver based on the concentrations of the non-PSS species 13.

3.2. Turbulence-chemistry interaction When turbulence-chemistry is accounted for, the integral in equation (7) has to be calculated in each cell at every iteration. A Gaussian quadrature-based approach is adopted to efficiently perform this integration numerically. The combination of equation (7) and (8) gives expression (9) for the mean reaction rate coefficient: V

= c = 2 ^ -> &/> S exp @− V

− > ch, B exp N− R U i 2 ->

(9)

A linear transformation of to the normalized variable Q is performed. The resulting expression for the mean reaction rate coefficient is given below.

Q=

− 2 -> &/>


(10)

12 of 49

Page 13 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


V

= ^ &/> S exp N− V

ch, > &/> i 2 - Q

+

R exp−Q > UQ

(11)

The right-hand side of this expression can be re-written as follows, which conforms with the template for Gaussian quadrature based integration:

^

&/>

V

S exp−Q V

>

V ch, &/> j exp N− Rk UQ = ^ S lQ mQ UQ i 2 -> &/> Q + V

with lQ = exp −Q > and mQ = exp N−

,

n> _`#

"/#

o

(12)

R

The mean reaction rate coefficient is hence numerically calculated using equation (13).

≅ ^

[ & > q r) ) )]&

exp N−

ch,)

i 2 - > &/> Q) +

R

(13)

A quadrature formula with N = 7 was found to give satisfactory results. The on-the-fly integration of equation (7) thus breaks down to the product of the 7 pre-calculated weights and the reaction rate evaluated at the 7 abscissae, i.e. temperatures. As the abscissae and weights are independent of the adopted reaction network, these values are taken from literature

53

. The abscissae Q) are equal to -2.651961, -1.673552, -0.8162879, 0.0, 0.816279,

1.673552 and 2.651961 with the corresponding weights r) being 0.0009717812, 0.05451558, 0.4256073, 0.8102646, 0.4256073, 0.05451558 and 0.0009717812 respectively.

3.3. Dynamic zoning method Calculation of the species rates of formation accounts for the largest share of CPU time. Hence, significant speed-up can be realized if the chemistry calculations would not need to be carried out on the same fine grid as the flow calculations. To evaluate the potential speed-up and the trade-off with simulation accuracy, the dynamic zoning partitioning scheme based on the work of Liang et al.

31

was tested. This scheme comprises three steps: grouping of


13 of 49


Page 14 of 49

“thermodynamically similar” cells into zones, calculation of the rates of the formation based on the zonal averages and mapping of the zonally averaged solution back to the individual cells. For combustion applications, several advanced algorithms

31, 54

have been proposed for

grouping of cells because of the strongly varying conditions during combustion, e.g. from highly stratified to near-homogeneous. This implies that, for certain outliers, e.g. cells where ignition starts, the algorithm must assure that the accuracy is maintained which results in a significant CPU overhead. In the case of pyrolysis and steam cracking, the thermodynamic state is more uniform over the grid compared to combustion applications because of the endothermic nature of the pyrolysis reactions. Therefore, a simple, fast, uniform zoning method can be applied. When disregarding turbulence-chemistry interaction, the v w . Hence, thermodynamic state that determines the rates of formation in cell is s t , u strictly speaking, grouping of cells should be performed by comparing all species concentrations and temperature. This is inapplicable for two reasons: the CPU time associated with comparing of all these variables would be prohibitively high and it is far from straightforward to define an appropriate criterion for similarity for each species individually as one should know the expected concentration span of each species in the reactor

55

.

However, as an approximation, the thermodynamic state can be simply thought of as being uniquely represented by a limited number of so-called features U 31, e.g. the temperature t and the yield of one characteristic species. This agrees with the conclusion of Van Geem et al. who stated that two so-called severity indices are sufficient to unambiguously characterize the product yield for a given feedstock

56

. Figure 1 shows the rate of formation for ethene as a

function of mean temperature and dry methane yield for all cells in the simulation of a typical butane-cracking reactor. Indeed, the rate of formation for ethene is a smooth function of mean temperature and dry methane yield and grouping of cells can be performed based on these two ACS Paragon Plus Environment

14 of 49

Page 15 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


features. In each iteration step, the minimum and maximum values of the selected features U are determined. Next, based on user-supplied thresholds * x on each feature U, the cells are uniformly grouped into zones as visualized in Figure 1. Depending on the values for the usersupplied thresholds, the number of zones in representative cases varies between 102 and 104. Based on the zonal averages, the species rates of formation and the heat of reaction are calculated for each zone 31. Finally, the species rates of formation and heat of reaction in each cell, necessary for solving equations (16) and (17), are set equal to the values of their associated zone.

4. Results and discussion 4.1. Impact of turbulence-chemistry interaction 4.1.1. Validation with DNS data In order to validate the usage of the temperature variance conservation equation, a comparison is made with the DNS data of non-isothermal turbulent pipe flow calculated by Redjem-Saad et al.

57

. As in the work of Redjem-Saad, the flow is modeled incompressible, and the

temperature is treated as a passive scalar. Following the methodology of Patankar et al.

58

,

streamwise periodic boundary conditions are applied. By applying a uniform heat flux to the tube wall, the fluctuating components of pressure and temperature can be isolated from the mean streamwise gradient. More details on the applied methodology can be found in Van Cauwenberge et al. 52. Figure 2 compares the results of the mean velocity (A), temperature (B), the turbulent kinetic energy (C) and the root-mean-square temperature (D) as a function of the dimensionless wall


15 of 49


Page 16 of 49

distance. Good agreement is seen for the mean variables, although both mean velocity and temperature are slightly overpredicted near the tube center. The turbulent kinetic energy is underpredicted near the wall and overpredicted in the center of the tube, which leads to considerable deviations in the simulated dimensionless root-mean-square temperature. The algebraic closure model for the temperature variance dissipation * calculates the correct magnitude for the root-mean-square temperature fluctuations, but the maximum is shifted towards the tube center. The simulations in which the temperature variance dissipation is calculated via a dedicated transport equation significantly overpredict the root-mean-square temperature. In other words, the dissipation rate is underpredicted. Using different sets of parameters for the transport equation for the dissipation of the temperature variance shows considerable differences, as exemplified in Figure 2 D, when using the parameters of Dibble et al. 59 and Elghobashi and LaRue 60. Since solving a dedicated conservation equation for the temperature variance dissipation does not yield realistic results, the algebraic closure model was used in the remainder of this work.

4.1.2. Turbulence-chemistry interaction To assess the importance of turbulence-chemistry interaction, a validation case is used that considers an n-butane-cracking reactor of 10.55 m long with an internal diameter of 0.0302 m and a wall thickness of 0.00675 m. An entrance zone of 2.0 m was simulated to ensure fully developed flow at the reactor inlet. The feedstock was pure n-butane at a mass flow rate of 140 kg h-1 and steam was added at a dilution of 0.30 kg kg-1. The coil outlet pressure was set to 233.15 kPa abs and the coil inlet temperature was set to 893 K. The adopted single-event microkinetic network has 20 molecules and 24 radicals. Both the fluid and the metal tube are modeled. The heat flux profile shown in Figure 3 was applied to the reactor outer wall.


16 of 49

Page 17 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4 (B) show the temperature standard deviation - as a function of the axial coordinate for the considered validation case. As the mean temperature profile in Figure 4 (A) is quite smooth, the source term 2

Jy

z2y

∇ > in equation (3), which depends on the mean temperature

gradient, does not generate a large uncertainty on the temperature, keeping the temperature standard deviation relatively low. The maximum temperature standard deviation is around 23 K. Towards the end of the reactor, the temperature standard deviation decreases as the temperature difference between the tube center and the inner wall decreases, and as the temperature variance dissipation * increases. The importance of turbulence-chemistry interaction depends largely on the existing radial temperature gradient in the reactor, which, in turn, is determined by the reactor diameter, convective heat transfer coefficient and heat flux. The effect of turbulence-chemistry interaction on the radial temperature profile is negligible as seen from Figure 5. The species rates of formation are slightly affected, having a small effect on the product yields as discussed further. However, the change in rates of formation is too small to influence the thermal power from reactions to the extent of having a significant influence on the mean temperature. Table 2 summarizes the results of the cases with and without accounting for turbulence-chemistry interaction. Note that only the results obtained with the algebraic temperature variance dissipation model are shown. The conversion is around to 90% to be representative for the industrial situation of n-butane cracking at moderate severity in a Millisecond reactor. Note that substantial conversion takes place in the adiabatic section and even in the beginning of the transfer line heat exchanger in the case of Millisecond furnaces. Considering turbulence-chemistry interaction results in a small change for the ethene and propene yield, i.e. the ethene yield increases by 0.15 wt% while the propene yield decreases


17 of 49


Page 18 of 49

by 0.07 wt%. The effect of turbulence-chemistry interaction on conversion and product yields is small, and results in a change in ethene and propene yields that is smaller than the effect of a three-dimensional reactor technology on these yields as shown in simulations of finned reactors and ribbed reactors 13, 17, 61.

4.2. Dynamic zoning method The dynamic zoning method was evaluated with the simulation of a butane-cracking reactor. The goal of these simulations is to select appropriate features U for the zoning algorithm and

to find user-specified thresholds * x for the features that balance accuracy against

performance in terms of CPU time. All simulations were performed on an Intel Xeon E5-2670 octo-core processor with 16 logical cores via Intel® Hyper-Threading Technology, i.e. a maximum of 16 threads was run simultaneously. As stated in section 3.3, the rate of formation of a species in a cell is determined by the

v w in that cell. Hence, an obvious choice for a first feature is the thermochemical state s t , u

mean temperature . As a second feature, the n-butane conversion { seems the most v . However, when conversion reaches straightforward to represent the concentrations vector u

100 % in the reactor, all cells with a conversion of 100 % will be grouped into a single zone, although the cracking severity in these cells can strongly vary. Indeed, secondary reactions, such as converting light olefins to aromatics, continue even after the feed molecule is fully consumed. Hence, the cell concentration vectors and the correct species rates of formation will differ significantly within this zone and will hence not be calculated accurately via the dynamic zone method with conversion as one of the features. The simulation results with nbutane conversion as second feature are included in the Supplementary Material but are not further discussed.


18 of 49

Page 19 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Rather than conversion, the methane weight fraction was evaluated as second zoning feature. To remove the effect of dilution on the zoning feature threshold, the dry methane weight fraction is used, i.e. normalized using the total hydrocarbon weight fraction. The dry methane weight fraction increases monotonically with cracking severity and has been proposed previously as a cracking severity index

56

. Figure 6 (A) and (B) show the simulated outlet

conversion and ethene yield respectively as a function of * |}~ . The error on conversion and consequently on the ethene yield compared to the fully resolved simulation is seen to be largely determined by * |}~ while the effect of * is rather small, except for the ethene yield at * |}~ = 100 r% and * = 1. Appropriate values for * and * |}~ are 10 K and 5 wt% respectively. Figure 7 shows the CPU time consumed by calculation of the reaction rates and rates of formation (chemistry time) and the overhead of the zoning algorithm (zoning time). The zoning overhead includes the gathering of the cell data from all CPUs, calculation of the number of zones, distribution of the cells over the different zones and calculation of the zone averages. The CPU time used for zoning and chemistry calculations in the simulation with * |}~ = 1 % and * = 1 K is comparable to the CPU time in the fully resolved case. For all other simulations, a speedup factor ranging from 3.7 to 125 is obtained. The simulation with * |}~ = 5 % and * = 5 which shows reasonable agreement to the fully resolved simulation as shown in Figure 6 has a speedup factor of 18. A significant overhead of the zoning algorithm is seen for all simulations, which is largely associated with the collection of the data from the different logical cores on one logical core to distribute the cells over their respective zones and calculation of the zonal averages. As v w are also geometrically close to one another, cells with a similar thermochemical state s t , u zoning can alternatively be performed logical core-wise, i.e. on each logical core the cells are


19 of 49


Page 20 of 49

grouped into zones based on the user-specified feature thresholds without inter-logical core communication. This methodology was tested for the parametric study using the methane weight fraction and mean temperature as features. Figure 8 (A) and (B) show the simulated outlet conversion and ethene yield respectively as a function of * |}~ . As zoning is performed per logical core, the number of zones is higher compared to the previous parametric studies. Hence, the agreement to the fully resolved simulation data is much better than in the previous parametric studies. For all simulations, the absolute difference to the outlet conversion in the fully resolved simulation is below 0.04 %. As before, the error on conversion and consequently on the ethene yield compared to the fully resolved simulation is seen to increase with increasing * |}~ while the effect of * is rather small. Figure 9 shows the CPU times consumed by chemistry calculation and by zoning per logical core for the 12 cases. Note that the bar for fully resolved simulation does not fit on this figure due to the choice of the y-axis range: the scale is different than in Figure 7 to be able to see the difference between the cases in the parametric study. Compared to the fully resolved simulation, a speedup factor from 50 to 190 is obtained. A much higher speedup factor is obtained compared to zoning over the entire grid. The reason is two-fold. First, the timeconsuming gathering of the data spread over the different logical cores is avoided. Second, the calculation of the zonal averages per core is faster as each zone has fewer cells. Approximatively the same total number of floating point operations is still required to calculate the zonal averages as when zoning over the entire grid is applied, however, in this case, the floating-point operations can be carried out in parallel as the information is processed on each core individually, so the total wall time for averaging is substantially lower. Little CPU time reduction is seen by increasing * |}~ above 5 wt% for a fixed * .


20 of 49

Page 21 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This is a result of the zoning being performed per logical core. As 16 logical cores are used and the maximum methane weight fraction is around 20 wt%, the average methane weight fraction range per logical core is around 1.25 wt%. Hence, zoning with a threshold * |}~ above 5 wt% results in almost no increase in the number of zones. To summarize, zoning per logical core has a higher accuracy and a better performance compared to zoning over the entire grid and is therefore recommended. Finally, it should be noted that, as in this case steady-state simulations are performed, the feature threshold values can be changed along the run time of the simulation. At the start of a simulation, high feature threshold values can be applied, i.e. * = 10 K and * |}~ = 100 wt% . This results in fast iterations that still guarantee a solution close to the fully resolved solution. Afterwards, the feature threshold values can be refined to * = 1 K and * |}~ = 1 wt% to move closer to the fully resolved solution. Ultimately, zoning can be disabled to calculate the rates of production in each cell individually and converge to the final solution.

4.3. Simulation of a large scale steam cracking reactor The advantages of the presented methodology are illustrated for a large-scale n-butanecracking U-coil reactor. A U-coil reactor has two passes, referred to as the inlet and outlet leg, respectively. Both legs are connected by a return bend and a joint where the diameter gradually expands from the inlet leg diameter to the outlet leg diameter. Figure 10 (A), (B) and (C) show a front, side and top view of the lower part of the reactor, i.e. focusing on the return bend. The first part of the return bend is an S-bend moving the reactor outside the plane defined by the inlet and outlet leg. Downstream, a U-bend connects the S-bend to the outlet leg. The most upstream part of the outlet leg is a joint where the diameter gradually expands. The reactor geometry details and process conditions are summarized in Table 3. To ensure ACS Paragon Plus Environment

21 of 49


Page 22 of 49

fully developed flow at the reactor inlet, an entrance zone, consisting of a 1.5 m long straight tube, is simulated upstream the inlet leg. The diameter of the simulated reactor is double than that of a typical reactor applied in industrial furnaces, which highlights the hydrodynamic effects. However the observations made based on the present simulation can be straightforwardly extrapolated to the industrially applied U-coil reactor. In this case, the adopted single-event microkinetic model contains 149 reactions between 11 molecules and 9 radicals and is obtained by reducing the butane pyrolysis network of Van de Vijver et al.

28

to its main species. The validity of the reduced network was confirmed by

comparison to pilot plant data, as shown in the Supporting Information. The adopted mesh is a structured butterfly grid created using Pointwise® V17.1 R4. The number of cells in the process gas fluid and reactor metal solid region is 1.976 ⋅ 10 and 2.040 ⋅ 10 respectively. An axial heat flux profile was applied to the reactor outer wall, based on a coupled furnacereactor simulations of a relevant section of the furnace as described in the work of Zhang et al. 62. The product yields are given in Table 3. Figure 11 shows the velocity magnitude and the in-plane velocity vectors for the six cross sections depicted on the right-hand side of the figure. The cross sections are labeled from A to F. Note that the in-plane velocity vectors are not scaled to their magnitude and hence only give the direction of the occurring secondary flow. Cross section A is located in the straight inlet leg and a fully developed turbulent flow field is simulated consequently. As the flow enters the S-bend, a low velocity zone is created near the inner curve of the bend and two counter-rotating vortices prevail. By entering the U-bend, the low-velocity region shifts to the inner curve of the U-coil bend from C to D. Downstream the U-bend, the two counter-rotating vortices and the dead zone remain for some length as shown in E and F. The velocity in E and F is also lower, compared to the other cross sections, as the inner diameter has enlarged and,


22 of 49

Page 23 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


hence, the cross-sectional area is larger. The effect of these flow patterns on temperature and conversion is apparent: low near-wall velocity results in low heat transfer and subsequently in high temperature and high conversion. The near-wall low velocity in a U-bend prevails at the inner curve 63-65. As the U-bend is not in the plane defined by the inner and outer leg, denoted ‘p1’ in Figure 10 C, but rather on plane ‘p2’, the near-wall low velocity zone in the outlet leg is not directed towards the inlet leg but is directed towards the start of the U-bend. Hence, it is in a location receiving a high radiative heat flux from the furnace walls and burners. Having this zone directed towards the inlet leg would be beneficial to counteract the high wall temperature caused by the low nearwall velocity as this is the so-called ‘shadow zone’ receiving less heat flux

66, 67

. Figure 12

(A), (B) and (C) show the velocity, temperature and n-butane yield in plane ‘p2’ shown in Figure 10 C in the outlet leg. As shown previously the velocity in the U-bend is higher near the outer curve of the bend. This results in a low velocity zone near the inner curve of the outlet tube that remains up to about 3 m in the outer leg. This in turn causes a high temperature near the inner curve of the outlet leg as shown in Figure 12 (B). Consequently, the n-butane weight fraction in a cross-section of the outlet leg is highly non-uniform. This spatial non-uniformity has its implications on coke formation as well. Due to the high temperature, low velocity and high concentration of precursors, the rate of coke formation near the inner curve of the outlet leg is expected to be high. The coke layer thickness will increase substantially faster than near the outer curve. Eventually, the coke layer will encompass the entire low-velocity zone, narrowing the cross-sectional area for flow. The thick coke layer on the inner curve of the outlet leg is susceptible to coke spalling, that is large pieces of coke breaking off from the wall, increasing the risk of coke accumulation in the


23 of 49


Page 24 of 49

bottom of the bend or the entrainment of fine coke particles to the separation section downstream of the furnaces.


24 of 49

Page 25 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5. Conclusions The importance of turbulence-chemistry interaction during steam cracking was assessed. To this end, a code for the three-dimensional simulation of steam cracking reactors, based on the free and open source CFD package OpenFOAM®, was developed. The code incorporates an algorithm based on a dynamic zoning method that speeds up the chemistry calculations by more than a factor ten, while maintaining a similar accuracy as calculating the chemistry in every cell individually. The effect of turbulent temperature fluctuations is accounted for by solving an extra continuity equation for temperature variance. Comparison of this approach to DNS data for non-reactive, turbulent pipe flow, showed that the root-mean-square temperature is underpredicted near the wall and overpredicted in the tube center. Simulation of reactive flow indicated that including the turbulence-chemistry interaction results in a correction on the ethene and propene yields of around 0.1-0.4 wt% absolute. The capabilities of the approach were demonstrated for a large-scale U-coil reactor fed with butane, which enabled us to visualize secondary flow patterns in the U-bend and an associated high temperature/high conversion zone in the outlet leg near the inner curve of the U-bend, undetectable with one-dimensional simulations. These detailed simulation results show that this code can be used for 3D optimization of reactor geometries.


25 of 49


Page 26 of 49

Acknowledgements Pieter A. Reyniers acknowledges financial support from a doctoral fellowship from the Research Foundation - Flanders (FWO). The authors also acknowledge the financial support from the Long Term Structural Methusalem Funding by the Flemish Government – grant number BOF09/01M00409. The computational work was carried out using the STEVIN Supercomputer Infrastructure at Ghent University, funded by Ghent University, the Flemish Supercomputer Center (VSC), the Research Foundation - Flanders (FWO) and the Flemish Government – Department EWI. The authors are grateful for the fruitful discussions with professor R.O. Fox from Iowa State University.

Supporting information -

Extended discussion on the model equations

-

Validation of the reduced kinetic networks

-

Dynamic zoning based on the mean temperature and the n-butane conversion as features.


26 of 49

Page 27 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nomenclature Abbreviations 1D 3D CFD CPU DNS LES P/E PDF PSS PSSA RANS SIMPLE SST

One-dimensional Three-dimensional Computational fluid dynamics Central processing unit Direct numerical simulation Large eddy simulation Propene/ethene Probability density function Pseudo-steady state Pseudo-steady state assumption Reynolds-averaged Navier Stokes Semi-Implicit Method for Pressure-Linked Equations Shear-stress transport

Roman symbols

Δ;,) ℎ H)

H

u & , > , I , K U 9 ch, mQ ℎ H

Reaction enthalpy [J kg-1] Pre-exponential factor of reaction i [mol m-3 s-1] Concentration of species j [mol m-3] Thermal heat capacity [J kg-1 K-1] Concentration vector Scalar variance model closure terms [-] Feature of the dynamic zoning method Diffusion coefficient [m s-1] Activation energy for reaction i [J mol-1] Gaussian quadrature abscissae Enthalpy [J] Turbulence kinetic energy [m2 s-2] Boltzmann constant [1.38064852 10-23 J K-1] Reaction rate coefficient of reaction i [mol m-3 s-1] Molar mass [kg mol-1] Number of species Number of Gaussian quadrature integration points Pressure [Pa]


27 of 49


T T i i) |4| 4H 45 ∗ = P r)

lQ Q {

Page 28 of 49

Probability density function as function of Prandtl number [-] Reaction rate of reaction i [mol m-3 s-1] Universal gas constant [8.3144598 J mol-1 K-1] Rate of formation of species j [kg m-3 s-1] Strain rate tensor Schmidt number [-] Volumetric thermal power [J m-3 s-1] Temperature [K] Dimensionless temperature [-] Dimensionless temperature in kinetic theory Velocity vector Dimensionless velocity [-] Gaussian quadrature weight point j Gaussian quadrature weights Normalized temperature variable in quadrature based integration Conversion [-] Distance perpendicular to the boundary surface [m] Dimensionless wall distance [-] Yield [-]

Greek symbols

,) *+, *

* x

* |}~ * ' % %,) -+, ->

Reaction order of species i in reaction j [-] Lennard-Jones well depth [J] Temperature scalar dissipation rate [K2 s-1] User-supplied threshold of feature U in the dynamic zoning method User-supplied threshold of methane yield in the dynamic zoning method User-supplied threshold of mean temperature in the dynamic zoning method Thermal conductivity [W m-1 K-1] Viscosity [kg m-1 s-1] Kinematic viscosity [m2 s-1] Stochiometric coefficient of species i in reaction j [-] Density [kg m-3] distance at which the intermolecular potential between two particles is zero [m] Temperature scalar variance [K2]


28 of 49

Page 29 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


̅̅ Ω ∗

Stress tensor [Pa] Wall shear stress [Pa] Temperature dependent intermediate value in kinetic theory

Subscripts and superscripts mm m

Effective Fluid Reaction index Species index Cell index Root mean square Mean


29 of 49


Page 30 of 49

References (1) (2) (3)

(4)

(5) (6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14) (15) (16)

Boileau, M.; Staffelbach, G.; Cuenot, B.; Poinsot, T.; Berat, C., LES of an ignition sequence in a gas turbine engine. Combust. Flame 2008, 154, 2-22. Esclapez, L.; Riber, E.; Cuenot, B., Ignition probability of a partially premixed burner using LES. Proc. Combust. Inst. 2015, 35, 3133-3141. Hakim, L.; Ruiz, A.; Schmitt, T.; Boileau, M.; Staffelbach, G.; Ducruix, S.; Cuenot, B.; Candel, S., Large eddy simulations of multiple transcritical coaxial flames submitted to a high-frequency transverse acoustic modulation. Proc. Combust. Inst. 2015, 35, 14611468. Hannebique, G.; Sierra, P.; Riber, E.; Cuenot, B., Large Eddy Simulation of Reactive Two-Phase Flow in an Aeronautical Multipoint Burner. Flow, Turbul. Combust. 2013, 90, 449-469. Irannejad, A.; Banaeizadeh, A.; Jaberi, F., Large eddy simulation of turbulent spray combustion. Combust. Flame 2015, 162, 431-450. Jones, W. P.; Marquis, A. J.; Wang, F., Large eddy simulation of a premixed propane turbulent bluff body flame using the Eulerian stochastic field method. Fuel 2015, 140, 514-525. Cuoci, A.; Frassoldati, A.; Stagni, A.; Faravelli, T.; Ranzi, E.; Buzzi-Ferraris, G., Numerical Modeling of NOx Formation in Turbulent Flames Using a Kinetic Postprocessing Technique. Energy Fuels 2013, 27, 1104-1122. Stagni, A.; Cuoci, A.; Frassoldati, A.; Faravelli, T.; Ranzi, E., A fully coupled, parallel approach for the post-processing of CFD data through reactor network analysis. Comput. Chem. Eng. 2014, 60, 197-212. Hassan, G.; Pourkashanian, M.; Ingham, D.; Ma, L.; Newman, P.; Odedra, A., Predictions of CO and NOx emissions from steam cracking furnaces using GRI2.11 detailed reaction mechanism – A CFD investigation. Comput. Chem. Eng. 2013, 58, 6883. Adhikari, S.; Sayre, A.; Chandy, A. J., In situ adaptive tabulation (ISAT) for combustion chemistry in a network of perfectly stirred reactors (PSRs). Comput. Chem. Eng. 2017, 97, 124-134. Blasi, J. M.; Kee, R. J., In situ adaptive tabulation (ISAT) to accelerate transient computational fluid dynamics with complex heterogeneous chemical kinetics. Comput. Chem. Eng. 2016, 84, 36-42. Partopour, B.; Dixon, A. G., Computationally efficient incorporation of microkinetics into resolved-particle CFD simulations of fixed-bed reactors. Comput. Chem. Eng. 2016, 88, 126-134. Reyniers, P. A.; Schietekat, C. M.; Van Cauwenberge, D. J.; Vandewalle, L. A.; Van Geem, K. M.; Marin, G. B., Necessity and Feasibility of 3D Simulations of Steam Cracking Reactors. Ind. Eng. Chem. Res. 2015, 54, 12270-12282. Zimmermann, H.; Walzl, R., Ethylene. In Ullmann's Encyclopedia of Industrial Chemistry, Wiley-VCH Verlag GmbH & Co. KGaA: 2000. Chapman, D. R., Computational Aerodynamics Development and Outlook. AIAA J. 1979, 17, 1293-1313. Choi, H.; Moin, P., Grid-point requirements for large eddy simulation: Chapman's estimates revisited. Phys. Fluids 2012, 24.


30 of 49

Page 31 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(17) Schietekat, C. M.; Van Cauwenberge, D. J.; Van Geem, K. M.; Marin, G. B., Computational Fluid Dynamics-Based Design of Finned Steam Cracking Reactors. AlChE J. 2014, 60, 794-808. (18) Schietekat, C. M.; van Goethem, M. W. M.; Van Geem, K. M.; Marin, G. B., Swirl flow tube reactor technology: An experimental and computational fluid dynamics study. Chem. Eng. J. 2014, 238, 56-65. (19) Stefanidis, G. D.; Merci, B.; Heynderickx, G. J.; Marin, G. B., CFD simulations of steam cracking furnaces using detailed combustion mechanisms. Comput. Chem. Eng. 2006, 30, 635-649. (20) Liu, J.; Chen, S.; Liu, Z.; Peng, K.; Zhou, N.; Huang, X.; Zhang, T.; Zheng, C., Mathematical Modeling of Air– and Oxy–Coal Confined Swirling Flames on Two Extended Eddy-Dissipation Models. Ind. Eng. Chem. Res. 2012, 51, 691-703. (21) Rogerson, J. W.; Swaminathan, N.; Tanahashi, M.; Shiwaku, N., Analysis of progress variable variance equations using DNS data. In Third European Combustion Meeting, Chania, Crete, 2007. (22) Chakraborty, N.; Swaminathan, N., Effects of Lewis Number on Scalar Variance Transport in Premixed Flames. Flow, Turbul. Combust. 2011, 87, 261-292. (23) Pei, Y.; Hawkes, E.; Kook, S., A Comprehensive Study of Effects of Mixing and Chemical Kinetic Models on Predictions of n-heptane Jet Ignitions with the PDF Method. Flow, Turbul. Combust. 2013, 91, 249-280. (24) Bazooyar, B.; Shariati, A.; Hashemabadi, S. H., Turbulent Non-premixed Combustion of Rapeseed Methyl Ester in a Free Shear Swirl Air Flow. Ind. Eng. Chem. Res. 2016, 55, 11645-11663. (25) Singh, K. D.; Dabade, T.; Vaid, H.; Gangadharan, P.; Chen, D.; Lou, H. H.; Li, X.; Li, K.; Martin, C. B., Computational Fluid Dynamics Modeling of Industrial Flares Operated in Stand-By Mode. Ind. Eng. Chem. Res. 2012, 51, 12611-12620. (26) Poinsot, T.; Veynante, D., Theoretical and Numerical Combustion. R.T. Edwards, Inc.: Philadelphia, PA, USA, 2005. (27) Lu, T.; Law, C. K., Toward accommodating realistic fuel chemistry in large-scale computations. Prog. Energy Combust. Sci. 2009, 35, 192-215. (28) Van de Vijver, R.; Vandewiele, N. M.; Bhoorasingh, P. L.; Slakman, B. L.; Seyedzadeh Khanshan, F.; Carstensen, H.-H.; Reyniers, M.-F.; Marin, G. B.; West, R. H.; Van Geem, K. M., Automatic Mechanism and Kinetic Model Generation for Gas- and Solution-Phase Processes: A Perspective on Best Practices, Recent Advances, and Future Challenges. Int. J. Chem. Kinet. 2015, 47, 199-231. (29) Dijkmans, T.; Schietekat, C. M.; Van Geem, K. M.; Marin, G. B., GPU based simulation of reactive mixtures with detailed chemistry in combination with tabulation and an analytical Jacobian. Comput. Chem. Eng. 2014, 71, 521-531. (30) Fang, Z.; Qiu, T.; Chen, B., Improvement of ethylene cracking reaction network with network flow analysis algorithm. Comput. Chem. Eng. 2016, 91, 182-194. (31) Liang, L.; Stevens, J. G.; Farrell, J. T., A Dynamic Multi-Zone Partitioning Scheme for Solving Detailed Chemical Kinetics in Reactive Flow Computations. Combust. Sci. Technol. 2009, 181, 1345-1371. (32) Weller, H. G.; Tabor, G.; Jasak, H.; Fureby, C., A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 1998, 12, 620631.


31 of 49


Page 32 of 49

(33) Van Geem, K. M.; Žajdlík, R.; Reyniers, M.-F.; Marin, G. B., Dimensional analysis for scaling up and down steam cracking coils. Chem. Eng. J. 2007, 134, 3-10. (34) Poling, B. E.; Prausnitz, J. M.; O'Connell, J. P., The properties of gases and liquids. McGraw-Hill: 2001. (35) Green, W. H.; Allen, J. W.; Ashcraft, R. W.; Beran, G. J.; Goldsmith, C. F.; Harper, M. R.; Jalan, A.; Magoon, G. R.; Matheu, D. M.; Petway, S.; Raman, S.; Sharma, S.; Van Geem, K. M.; Song, J.; Wen, J.; West, R. H.; Wong, A.; Wong, H.-W.; Yelvington, P. E.; Yu, J. RMG - Reaction Mechanism Generator v3.3. http://rmg.sourceforge.net/ (36) Fox, R. O., Computational Models for Turbulent Reacting Flows. Cambridge University Press: 2003. (37) Béguier, C.; Dekeyser, I.; Launder, B. E., Ratio of scalar and velocity dissipation time scales in shear flow turbulence. Phys. Fluids (1958-1988) 1978, 21, 307-310. (38) Pope, S. B., Consistent modeling of scalars in turbulent flows. The Physics of Fluids 1983, 26, 404-408. (39) Jones, W. P.; Musonge, P., Closure of the Reynolds stress and scalar flux equations. Phys. Fluids 1988, 31, 3589-3604. (40) Sanders, J. P. H.; Gokalp, I., Scalar dissipation rate modelling in variable density turbulent axisymmetric jets and diffusion flames. Phys. Fluids 1998, 10, 938-948. (41) Caracotsios, M.; Stewart, W. E., Sensitivity analysis of initial value problems with mixed odes and algebraic equations. Comput. Chem. Eng. 1985, 9, 359-365. (42) Turányi, T., Sensitivity analysis of complex kinetic systems. Tools and applications. J. Math. Chem. 1990, 5, 203-248. (43) Turányi, T.; Bérces, T.; Vajda, S., Reaction rate analysis of complex kinetic systems. Int. J. Chem. Kinet. 1989, 21, 83-99. (44) Kee, R.; Rupley, F.; Miller, J.; Coltrin, M.; Grcar, J.; Meeks, E.; Moffat, H.; Lutz, G.; Dixon-Lewis, A.; Smooke, M.; Warnatz, J.; Evans, G.; Larson, R.; Mitchell, R.; Petzold, L.; Reynolds, W.; Caracotsios, M.; Stewart, W.; Glarborg, P.; Wang, C.; Adigun, O.; Houf, W.; Chou, C.; Miller, S.; Ho, P.; Young, D. CHEMKIN Release 401. Reaction Design Inc. . CA: San Diego, 2004. (45) Brock, E. E.; Savage, P. E.; Barker, J. R., A reduced mechanism for methanol oxidation in supercritical water. Chem. Eng. Sci. 1998, 53, 857-867. (46) Sabbe, M. K.; De Vleeschouwer, F.; Reyniers, M. F.; Waroquier, M.; Marin, G. B., First Principles Based Group Additive Values for the Gas Phase Standard Entropy and Heat Capacity of Hydrocarbons and Hydrocarbon Radicals. J. Phys. Chem. A 2008, 112, 12235-12251. (47) Sabbe, M. K.; Saeys, M.; Reyniers, M. F.; Marin, G. B.; Van Speybroeck, V.; Waroquier, M., Group additive values for the gas phase standard enthalpy of formation of hydrocarbons and hydrocarbon radicals. J. Phys. Chem. A 2005, 109, 7466-7480. (48) Patankar, S. V., Numerical Heat Transfer and Fluid Flow. Taylor & Francis: 1980. (49) Sweby, P. K., High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws. SIAM Journal on Numerical Analysis 1984, 21, 995-1011. (50) Jasak, H.; Weller, H. G.; Gosman, A. D., High resolution NVD differencing scheme for arbitrarily unstructured meshes. Int. J. Numer. Methods Fluids 1999, 31, 431-449. (51) Menter, F. R., Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994, 32, 1598-1605.


32 of 49

Page 33 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(52) Van Cauwenberge, D. J.; Schietekat, C. M.; Floré, J.; Van Geem, K. M.; Marin, G. B., CFD-based design of 3D pyrolysis reactors: RANS vs. LES. Chem. Eng. J. 2015, 282, 66-76. (53) Love, C. H., Abscissas and weights for gaussian quadrature. U.S. Department of Commerce, National Bureau of Standards: 1966. (54) Perini, F., High-dimensional, unsupervised cell clustering for computationally efficient engine simulations with detailed combustion chemistry. Fuel 2013, 106, 344-356. (55) Shi, Y.; Hessel, R. P.; Reitz, R. D., An adaptive multi-grid chemistry (AMC) model for efficient simulation of HCCI and DI engine combustion. Combust. Theory Modell. 2009, 13, 83-104. (56) Van Geem, K. M.; Reyniers, M.-F.; Marin, G. B., Two Severity Indices for Scale-Up of Steam Cracking Coils. Ind. Eng. Chem. Res. 2005, 44, 3402-3411. (57) Redjem-Saad, L.; Ould-Rouiss, M.; Lauriat, G., Direct numerical simulation of turbulent heat transfer in pipe flows: Effect of Prandtl number. Int. J. Heat Fluid Flow 2007, 28, 847-861. (58) Patankar, S. V.; Liu, C. H.; Sparrow, E. M., Fully Developed Flow and Heat Transfer in Ducts Having Streamwise-Periodic Variations of Cross-Sectional Area. J. Heat Transfer 1977, 99, 180-186. (59) Dibble, R. W.; kollman, W.; Farshchi, M.; Schefer, R. W., Second-order closure for turbulent nonpremixed flames: scalar dissipation and heat release effects. Symp. (Int.) Combust., [Proc.] 1988, 21, 1329-1340. (60) Elghobashi, S. E.; LaRue, J. C., The effect of mechanical stream on the dissipation rate of a scalar variance. In Proceedings of the 4th Symposium on Turbulent Shear Flows, Karlsruhe, 1983. (61) Van Cauwenberge, D. J.; Vandewalle, L. A.; Reyniers, P. A.; Floré, J.; Van Geem, K. M.; Marin, G. B., Periodic reactive flow simulation: Proof of concept for steam cracking coils. AlChE J. 2017. (62) Zhang, Y.; Qian, F.; Schietekat, C. M.; Van Geem, K. M.; Marin, G. B., Impact of flue gas radiative properties and burner geometry in furnace simulations. AlChE J. 2015, 61, 936-954. (63) Patankar, S. V.; Pratap, V. S.; Spalding, D. B., Prediction of Turbulent-Flow in Curved Pipes. J. Fluid Mech. 1975, 67, 583-595. (64) Sudo, K.; Sumida, M.; Hibara, H., Experimental investigation on turbulent flow through a circular-sectioned 180 degrees bend. Exp. Fluids 2000, 28, 51-57. (65) Sugiyama, H.; Hitomi, D., Numerical analysis of developing turbulent flow in a 180° bend tube by an algebraic Reynolds stress model. Int. J. Numer. Methods Fluids 2005, 47, 1431-1449. (66) Heynderickx, G. J.; Cornelis, G. G.; Froment, G. F., Circumferential tube skin temperature profiles in thermal cracking coils. AlChE J. 1992, 38, 1905-1912. (67) Heynderickx, G. J.; Froment, G. F., Simulation and Comparison of the Run Length of an Ethane Cracking Furnace with Reactor Tubes of Circular and Elliptical Cross Sections. Ind. Eng. Chem. Res. 1998, 37, 914-922.


33 of 49


Page 34 of 49

Figures

Figure 1: Rate of formation of ethene as a function of methane yield and mean temperature in a butane cracking reactor.


34 of 49

Page 35 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Validation of turbulent-chemistry interaction; Dimensionless mean velocity u+ (A); dimensionless mean temperature T+ (B), dimensionless turbulent kinetic energy k+ (C) and dimensionless root-mean-square temperature (D) as a function of dimensionless wall distance y+:

- DNS Redjem-Saad et al. 57;

from equation (4); LaRue 60;

- u+ = y+;

- law of the wall;

- RANS *

- RANS * from equation (5), parameters from Elghobashi and

- RANS * from equation (5), parameters from Dibble et al. 59. Reprinted

(Adapted or Reprinted in part) with permission from ref. 57. Copyright 2007 Elsevier Inc.


35 of 49


Page 36 of 49

Figure 3: Benchmarking of turbulent-chemistry interaction; heat flux profile on the reactor outer wall as function of the axial coordinate.


36 of 49

Page 37 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: Benchmarking of turbulent-chemistry interaction; (A) mean temperature [K] as function of axial coordinate [m]; symbols: without turbulence-chemistry interaction, line: with turbulence-chemistry interaction; (B) temperature standard deviation [K] as a function of axial coordinate [m].


37 of 49


Page 38 of 49

Figure 5: Benchmarking of turbulent-chemistry interaction; mean temperature [K] as a function of radial coordinate [10-3 m] at an axial position of 2 m (red, rectangles), 5 m (orange, diamonds) and 10 m (green, circles); symbols: without turbulence-chemistry interaction, line: with turbulence-chemistry interaction.


38 of 49

Page 39 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: Benchmarking of the dynamic zoning method with cell clustering over the entire grid based on mean temperature and methane yield; Outlet conversion (A) and ethene yield (B) as a function of * |}~ :

- fully resolved,

- * = 1 ,

- * = 10 .


- * = 5 ,

39 of 49


Page 40 of 49

Figure 7: Benchmarking of the dynamic zoning method with cell clustering over the entire grid based on mean temperature and methane yield; CPU time per iteration [s] for all

simulations in the parametric study using and ¡ :

- chemistry time and

-

zoning time.


40 of 49

Page 41 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8: Benchmarking of the dynamic zoning method with cell clustering per logical core based on mean temperature and methane yield; Outlet conversion (A) and ethene yield (B) as a function of * |}~ :

- fully resolved,

- * = 1 ,

* = 10 .


- * = 5 ,

-

41 of 49


Page 42 of 49

Figure 9: Benchmarking of the dynamic zoning method with cell clustering per logical core based on mean temperature and methane yield; CPU time per iteration [s] for all simulations

in the parametric study using and ¡ :

- chemistry time and

- zoning time. The

time for the fully resolved simulation continues beyond the top of the figure.


42 of 49

Page 43 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10: Simulation of a large-scale U-coil reactor for steam cracking of butane; Front view (A), side view (B) and top view (C).


43 of 49


Page 44 of 49

Figure 11: Simulation of a large-scale U-coil reactor for steam cracking of n-butane; Velocity magnitude [m s-1] in the six cross-sections in the return bend. Vectors show the in-plane velocity, i.e. the velocity component perpendicular to the plane is not considered when drawing the vectors.


44 of 49

Page 45 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12: Simulation of a large-scale U-coil reactor for steam cracking of n-butane; Velocity magnitude [m s-1] (A), mean temperature [K] (B) and n-butane [wt%] (C) in plane ‘p2’ in the outlet leg, viz. Figure 10 (C).


45 of 49


Page 46 of 49

Tables Table 1: Steady-state conservation equations for the fluid phase and the solid phase. Fluid phase Global mass Momentum Energy Species

= 0 ∇ =

= = −∇ + ∇ ∙ ¢£ ∇∙ =

;,:;; 1 > @ℎ + |= | BR = ∇ ∙ N ∇∙N = ∇ℎR + 45 2 H ) = −∇ ∙ 9:;;,) ∇ ) + i) ∇ ∙ =

(14) (15) (16) (17)

Solid phase Energy

∇ ∙ 6 ∇ = 0


(18)

46 of 49

Page 47 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Benchmarking of turbulence-chemistry interaction; Coil outlet temperature, P/E ratio, conversion and product yields. Turbulence-chemistry Coil outlet temperature [K] P/E ratio [wt%/wt%] n-butane conversion [%] Product Yields [wt%] H2 CH4 C2H6 C3H8 C2H4 C3H6 Product Selectivities [wt%] H2 CH4 C2H6 C3H8 C2H4 C3H6

Without 1148.90 0.479 90.20

With 1148.03 0.475 90.43

0.83 19.30 4.50 0.39 37.00 17.71

0.82 19.40 4.53 0.40 37.15 17.64

0.92 21.39 4.98 0.44 41.02 19.64

0.91 21.45 5.01 0.44 41.08 19.51


47 of 49


Page 48 of 49

Table 3: Simulation of a large-scale U-coil reactor for steam cracking of n-butane; Geometry details and process conditions. Reactor geometry Reactor type Inlet leg Length [m] Inner diameter [m] Wall thickness [m] S-bend Length [m] Inner diameter [m] Wall thickness [m] U-bend Length [m] Inner diameter [m] Wall thickness [m] Outlet leg Length [m] Inner diameter [m] Wall thickness [m] Process conditions Butane mass flow rate [kg s-1 reactor-1] Dilution [kg kg-1] Coil inlet temperature [K] Coil outlet pressure [103 Pa abs] Results Coil outlet temperature [K] Product yields [wt%] H2 CH4 C2H6 C3H8 n-C4H10 C2H4 C3H6 1,3-C4H6 1-C4H8 2-C4H8

U-coil 9.15 0.09 0.0116 1.0177 0.09 0.0116 2.32 0.09 0.0116 10.2 0.102 0.0156 0.4215 0.30 853 243.18 1129 0.50 16.55 8.18 1.01 11.00 38.78 20.61 2.47 0.61 0.29


48 of 49

Page 49 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For Table of Contents Only


49 of 49

CFD simulations of Industrial Steam Cracking Reactors: Turbulence

Recommend Documents