Catalyst simulations based on coupling of 3D CFD tool with effective 1D channel models

Catalysis Today 188 (2012) 87–93

Contents lists available at SciVerse ScienceDirect

Catalysis Today journal homepage: www.elsevier.com/locate/cattod

Catalyst simulations based on coupling of 3D CFD tool with effective 1D channel models ˇ epánek a , P. Koˇcí a,∗ , M. Marek a , M. Kubíˇcek b J. Stˇ a b

Institute of Chemical Technology, Prague, Department of Chemical Engineering, Technická 5, CZ 166 28 Praha, Czech Republic Institute of Chemical Technology, Prague, Department of Mathematics, Technická 5, CZ 166 28 Praha, Czech Republic

a r t i c l e

i n f o

Article history: Received 11 August 2011 Received in revised form 21 December 2011 Accepted 26 January 2012 Available online 30 March 2012 Keywords: 3D CFD Monolith NSRC LNT Simulation Exhaust gas aftertreatment

a b s t r a c t Spatially 1D simulations of exhaust gas catalysts are often used in the automotive industry because of short computation times. In fully 3D CFD-based model, the computational requirements increase by several orders of magnitude. A novel approach to 3D simulations of complete exhaust gas aftertreatment systems is proposed in this paper. It is based on a modular combination of standard CFD tool with a user-defined 1D channel model. The employment of external 1D channel model enables more effective computations and gives maximum freedom in definition of realistic reaction kinetics, including nonlinear rate laws and fully transient solution of adsorbed components. The approach is validated by comparison of the results from the stand-alone 1D model and single channel CFD simulation. Dynamic simulations of CO oxidation light-off and lean/rich operation of NOx storage catalyst (NSRC, LNT) are then conducted in a test geometry with 10 representative channels and skewed inlet/outlet pipes. Finally, a complex system including bent pipes and conical expansion/constriction parts is simulated with 80 representative channels. It is demonstrated that relatively large deviations can exist between the individual channels. Non-uniform inlet flow distribution affects light-off behaviour as well as breakthrough of adsorbed components, and it may also lead to a different extent of catalyst ageing in individual channels. The combined 3D CFD and 1D channels simulation approach represents an effective tool for analysis of such problems and can help with the design of more efficient exhaust gas aftertreatment systems. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The increasing computer power brings significant benefits for design of catalytic reactors. In recent years, spatially 1D simulations of monolith converters have been commonly used in the automotive industry [1]. The spatially 1D models provide usable results within a short computation time, although they have certain disadvantages. A single channel represents the whole monolith converter, therefore the 1D model cannot effectively describe the influence of flow variations through the individual channels of the monolith, resulting from inlet/outlet pipes geometry and monolith housing shape, as has been demonstrated with 2D and 3D monolith models [2–7].

Abbreviations: 1D, spatially 1-dimensional; 3D, spatially 3-dimensional; CFD, computational fluid dynamics; DOC, diesel oxidation catalyst; LNT, lean NOx trap (=NSRC); NH3 -SCR, selective catalytic reduction of NOx by ammonia; NSRC, NOx storage and reduction catalyst (=LNT); TWC, three-way catalyst. ∗ Corresponding author. Tel.: +420 22044 3293; fax: +420 22044 4320. E-mail address: [email protected] (P. Koˇcí). URL: http://www.vscht.cz/monolith (P. Koˇcí). 0920-5861/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cattod.2012.01.038

However, when all hydraulic flow properties, energy exchange, and catalytic reactions inside monolith channels are described in all three spatial dimensions using CFD, a large system of equations needs to be solved and such computations are extremely computerresource consuming [8]. Parallel systems with large amounts of memory are necessary to conduct such studies, and yet the processor time demand is considerable (a single simulation may take several days). Further disadvantage of many CFD solvers is the lack of support for realistic reaction kinetics on the catalytic surface. In the automotive exhaust gas aftertreatment applications, the reaction kinetics are strongly non-linear and often involve highly transient evolution adsorbed species, which is quite difficult to implement directly into a conventional CFD software. In this paper, we present an alternative approach enabling efficient 3D simulations of exhaust gas aftertreatment systems. It is based on a modular combination of standard 3D CFD software with a user-defined 1D channel model that enables more effective solution of the processes inside the catalytic monolith channel. The governing CFD software calculates the flow field in the pipes and heat transfer between individual monolith channels and manages the calls of 1D model for individual representative channels. All processes inside the monolith channels are then simulated by the 1D

J. Sˇ tˇepánek et al. / Catalysis Today 188 (2012) 87–93

88

2. Methodology Nomenclature Latin letters a density of external surface area in monolith, m2 m−3 c concentration in bulk gas, mol m−3 cs concentration in washcoat pores, mol m−3 g cp heat capacity of the gas, J kg−1 K−1 s cp effective heat capacity of the solid, J kg−1 K−1 H enthalpy, J mass transfer coefficient, m s−1 kc kh heat transfer coefficient, J m−2 K−1 s−1 Rj rate of jth reaction, mol m−3 s−1 (catalytic washcoat) sm mass source term in the 3D CFD model, kg m−3 s−1 se energy source term in the 3D CFD model, W m−3 t time, s T temperature, K Ts temperature of the solid phase, K u linear gas velocity in the 3D CFD model, m s−1 v linear gas velocity in the 1D channel model, m s−1 x spatial coordinate in the 3D CFD model yk mol. fraction of gas component k, 1 z spatial coordinate along the monolith channel in the 1D model, m Greek letters Hr standard reaction enthalpy, J mol−1 g ε void fraction of the reactor (fraction of open frontal area), 1 washcoat porosity, 1 εs s volume fraction of catalytic washcoat in the solid phase, 1  heat conductivity, W m−1 K−1  stoichiometric coefficient, mol mol−1 surface site coverage, 1  cap concentration of active surface sites, mol m−3

model. This approach represents a compromise between classical spatially 1D modelling and a fully 3D simulation. The external 1D channel model enables the user to include tailored numerical solution methods leading to faster computation times, virtually arbitrary reaction kinetics including non-linear rate laws with fully transient solution of adsorbed species (this feature is often missing in conventional CFD tools). It also ensures the compatibility between the stand-alone 1D and fully 3D simulations, which is very important for the design of exhaust gas aftertreatment systems. The parallelization methods available in the CFD software can be easily employed to further decrease the computation time for systems with many representative channels.

The simulation uses equidistant time nodes with time steps typically 0.1–0.2 s (generally the time step can be regulated adaptively). Within each time step, the CFD software calculates gas flow, pressure and temperature distribution. The calculated data are passed into the 1D channel model, together with the stored flow, concentration and temperature profiles along the channel at the beginning of the time step. The 1D solver then integrates evolution of concentration, pressure and temperature profiles, and enthalpy flux into the solid phase along each monolith channel. The intra-channel enthalpy flux (consisting of gas–solid heat transfer and reaction heat at each axial position inside the channel) is computed in the 1D model, using an axial solid-phase temperature profile provided by the 3D thermal conduction model of the monolith that is solved by the CFD software. The axial discretization of the monolith in CFD software matches the discretization of the 1D model, and the solid-phase temperature passed to the 1D model is the volume average of the cells in the solid that correspond to the current representative channel at each axial position. In turn, the 1D model passes to the 3D conduction simulation the above mentioned intra-channel heat source profile (the term se at each axial position). This enthalpy source profile is applied in the conduction model over the cells that are associated with the space of the current channel. This process is iterative – within each iteration, the calculated pressure losses and energy balances obtained from the 1D channel model and the CFD are compared for each representative channel and corrected until they match. When the compared values match within a given tolerance, next time step is initiated. The interaction between the CFD software and the 1D model is illustrated in Fig. 1. Calls to the 1D channel computation do not need to be performed with each 3D CFD iteration, typically one call for every four iterations is used until the CFD part is converged (each time step). The relative time spent computing the 3D CFD part and the 1D channel part depends on the relative complexity of these two parts (mainly the number of reactions and components in the reaction kinetics, and cell count in the CFD model) as well as the number of representative channels. In general, the channel part of the computation has never been observed to be more than 50% of the computation time, and most often it is between 5 and 20%. By default, the representative channels are distributed uniformly across the monolith, but the association of non-equal volumes is possible. The number of representative channels affects the solution in the following ways: (i) Higher numbers of representative channels provide higher resolution of the variation in state across the cross section of the device, but increase the computation time. (ii) Too few representative channels may cause problems with the convergence of the CFD problem. This depends on the degree of variation of the flow state over the inlet of the device.

Fig. 1. Scheme of interactions in the coupled CFD and 1D-channel models.

J. Sˇ tˇepánek et al. / Catalysis Today 188 (2012) 87–93

89

Fig. 2. Single channel test system geometry.

2.1. 3D CFD model Recently, a suitable interface enabling the use of external 1D channel models has been developed in the frame of CFD software package StarCD (plugin ES-AFTER [9]). The CFD solution is based on mass (1) and momentum (2) balance describing general fluid flow, along with energy balances (3, 4). In Cartesian coordinates, the equations are written in the following form [9,10]: Fig. 4. Test system geometry with computed linear velocity profile.

∂ ∂ ( uj ) = sm + ∂xj ∂t

(1)

∂( ui ) ∂ ∂p ( uj ui − i,j ) = − + + si ∂t ∂xj ∂xi

(2)

balance of the flowing gas (8) and the enthalpy balance of the solid phase (9):

∂( H) ∂ ∂p ( uj H + Fh,j − ui i,j ) = − + + si ui + sh , ∂t ∂xj ∂t

(3)

∂ck (z, t) ∂(v · ck ) kc a s =− + g c(yk − yk ), ε ∂t ∂z

where H = (1/2)ui ui + h.

∂cks (z, t)

∂( e) ∂Fe,j + se , = ∂t ∂xj

∂t

(4)

k = 1, . . . , K

1 kc a c(yk − yks ) + s k,j Rj , ε εs (1 − εg )ϕs

(5)

J

=

k = 1, . . . , K

j=1

(6)

where the term Fe,j describes thermal conductivity, and se is local enthalpy source (consisting of gas–solid heat transfer and reaction heats inside the channel).

∂

m (z, t)

∂t

=

1

J 

cap

m

m,j Rj ,

m = 1, . . . , M

(7)

j=1

2.2. 1D channel model Heterogeneous, spatially 1D plug-flow model of catalytic monolith channel with surface deposition of components [11,12] has been used for the simulations of the processes inside the representative channels. The model considers the following balances: mass balances of individual components in the flowing gas (5), in the washcoat pores (6) and on the catalyst surface (7), total enthalpy

a)

∂T (z, t) ∂T g g kh a s = −v cp + g (T − T ) ε ∂t ∂z

s cps

 kh a ∂T s (z, t) ∂ Ts + (T − T s ) − ϕs Hr,j Rj = s g 2 1 − ε ∂t ∂z

g

Single channel 1D+CFD 1D model

530 Single channel 1D+CFD 1D model

525 520 (K)

25 515

out

20

T

(ppm) out

(9)

j=1

30

NO

(8) J

2

b)

40 35

g cp

510

15 505

10

500

5 0 250

255

260 Time (sec)

265

270

495 250

260

270

280 290 Time (sec)

300

310

320

Fig. 3. Outlet NO concentration (a) and gas temperature (b) during NSRC lean/rich cycling in stand-alone 1D model and single channel CFD + 1D-channels test system. Space velocity GHSV = 60 000 h−1 , Tin = 498 K.

J. Sˇ tˇepánek et al. / Catalysis Today 188 (2012) 87–93

90

Boundary conditions used at the inlet (z = 0) and at the outlet (z = L) of the monolith are:

640

620

Tout (K)

600

Channel 1 Channel 2 Channel 6 Channel 10 1D model single channel geometry Inlet

T = T in

580

∂T s

560

∂z

a)

520 900

950

1000

1050

1100 Time (sec)

1150

1200

1250

1300

0.006

0.005

0.003

CO

out

(%)

0.004

0.002

0.001

b)

0 900

Channel 1 Channel 2 Channel 6 Channel 10 1D model single channel geometry Inlet 950

1000

(10)

= 0 at z = 0 z = L

ck = ckin ,

540

at z = 0

k = 1...K

(11)

at z = 0

(12)

The values of spatially distributed mass and heat transfer coefficients along the monolith channel (kc (z) and kh (z), respectively) are calculated from the correlations proposed in [13]. The reaction rates Rj for NOx storage and reduction catalyst are taken after [12]. The complete 1D channel model including detailed description of reaction kinetics (Rj ) in standard exhaust gas aftertreatment catalysts (TWC, DOC, NSRC and NH3-SCR) is implemented in Fortran into a software package XMR, developed in our research group [14]. 3. Validation of the approach

1050

1100

1150

1200

1250

1300

Time (sec)

Fig. 5. CO oxidation light-off curve, comparison of the calculated 1D-model and the CFD + 1D-channels test geometry: outlet temperature (a) and CO concentration (b) around the light-off temperature. Space velocity GHSV = 60 000 h−1 .

For the purposes of the validation, a simple device consisting of a single channel has been constructed (cf. Fig. 2). The approach has been validated using simulations of controlled lean-rich cycling with a NOx storage and reduction catalyst (NSRC) [12]. The considered inlet gas composition is described in Table 1. Fig. 3 illustrates the comparison of the stand-alone 1D model results with the

45 Channel 1 Channel 2 Channel 6 Channel 10 1D model

40 35 NOout (ppm)

30 25 20 15 10 5

a)

0 0

50

100

150

200

250

300

Time (sec) 45 Channel 1 Channel 2 Channel 6 Channel 10 1D model

40 35

NO

out

(ppm)

30 25 20 15 10 5

b)

0 250

252

254

256

258 Time (sec)

260

262

264

530 Channel 1 Channel 2 Channel 6 Channel 10 1D model

525

Tout (K)

520 515 510 505 500

c)

495 250

260

270

280 Time (sec)

290

300

310

Fig. 6. Outlet NO concentration evolution during NSRC lean/rich cycling: global view (a), detail of the 4th regeneration phase (b), and detail of gas temperature evolution (c) for various representative channels of the test system geometry, compared with stand-alone 1D-model results. Space velocity GHSV = 60 000 h−1 .

J. Sˇ tˇepánek et al. / Catalysis Today 188 (2012) 87–93

91

Table 1 Inlet gas composition and timing used for NSRC lean/rich cycling. Component

Lean phase

Rich phase

CO (%) O2 (%) H2 (%) C3 H6 (ppm) H2 O (%) CO2 (%) NO (ppm)

(s)

0 9.5 0 0 7 7 150 60

3.3 0.5 1.1 3000 7 7 150 5

Fig. 8. Temperature profile in the converter 1 s after the beginning of the 4th rich phase of the NSRC lean/rich cycling. Space velocity GHSV = 60 000 h−1 , Tin = 498 K.

Fig. 7. Used capacity of NOx storage sites on the catalyst surface at the end of the 1st lean phase (a) and at the end (b) of the 5th lean phase of NSRC regeneration. Space velocity GHSV = 60 000 h−1 , Tin = 498 K.

CFD + 1D-channels approach. It can be seen that the calculated results match well. 4. Results and discussion For purposes of simulating the influence of flow distribution on exhaust gas conversion and temperature, a test system geometry with skewed inlet and outlet pipes has been constructed with 10 representative channels, cf. Fig. 4. The CO oxidation light-off (temperature ramp from 373 K to 773 K at 10 K/min; for composition see Table 2) has been calculated and compared to the stand-alone 1D-model, cf. Fig. 5. The outlet values differ mainly around the light-off temperature, where the gas flow non-uniformities cause different flow rates for particular channels, resulting in different temperatures and conversions. Furthermore, periodic operation of NSRC has been simulated using the test geometry. The inlet mixture composition was the same as in the validation case (see Table 1). In Fig. 4, the calculated Table 2 CO oxidation light-off. Component

Concentration

CO (%) O2 (%)

0.5 1.0

Fig. 9. Complex system geometry including bent inlet pipe and expansion/constriction cones. Different colors illustrate the division of the system into individual regions (representative channels, pipe parts) eligible for parallel computations.

92

J. Sˇ tˇepánek et al. / Catalysis Today 188 (2012) 87–93

Fig. 10. (a) Full diesel oxidation catalyst geometry including pipes with a shown gas velocity profile; (b) gas temperature profile; (c) and (d) profile of gas CO concentration.

profile of linear gas velocities in the inlet and outlet pipes is shown 250 s after the simulation start (inlet gas temperature Tin = 498 K). It can be seen that the used geometry leads to a slightly non-uniform flow distribution. The striped area marks the domain, computed by the 1D model. Fig. 6 illustrates the outlet NO concentrations and temperatures during NSRC lean/rich cycling. Again, due to non-uniform flow distribution, the outlet concentrations and temperatures from the individual channels differ. The outlet temperature peak is much broader than the concentration peak due to substantial heat capacity of the monolith solid phase. The standard 1D model of the converter gives results close the centre of the range for individual channels from the CFD + 1D-channel simulation, which again demonstrates a correct coupling of the models. It can be seen that the differences in NO breakthrough from individual channels increase with time, cf. Fig. 6a. It is caused by the fact that the regeneration of the surface is not complete at this

relatively lower temperature, and the amount of the unreduced NOx remaining adsorbed on the surface is also different for individual channels. Fig. 7 illustrates the amount of NOx stored on the catalyst surface at the end of the 1st and the 5th lean phase. The surface gets more saturated with each consequent lean phase, which results in gradually increasing NOx breakthrough. For further details regarding the reaction kinetics, see [12]. Fig. 8 shows the temperature profile in the converter 1 s after the start of the 3rd rich phase of NSRC regeneration. A temperature increase can be observed at the inlet of the converter, which is caused by reaction heat of the rich mixture with high CO, H2 and HC concentrations. Finally, a dynamic simulation of CO oxidation light-off has been performed with 80 representative channels within a fully 3D converter geometry, including expansion and constriction cone and a relatively sharply bent inlet pipe (cf. Fig. 9).

J. Sˇ tˇepánek et al. / Catalysis Today 188 (2012) 87–93

non-uniform inlet flow distribution. This non-uniformity affects light-off behaviour as well as breakthrough of adsorbed components, and it may also lead to a different extent of catalyst ageing in individual channels. The combined 3D CFD and 1D channels simulation approach represents an effective tool for analysis of such problems and can help with the design of more efficient exhaust gas aftertreatment systems.

0.6

0.5

0.4 COout (%)

93

0.3

0.2

Acknowledgements

0.1

0 2

4

6

8

10

12

14

Time (sec) Fig. 11. CO outlet concentration evolution for various channels of the full monolith.

The considered inlet gas composition is listed in Table 2. The inlet gas temperature was 573 K, while the initial temperature of the solid phase was 293 K. This lead to a fast heating of the monolith front part, and local CO light-off. Fig. 10a depicts the flow field with significant non-uniformities in front of the monolith, involving a swirl. This causes low channel inlet velocities in certain parts of the monolith. The flow distribution affects the heat-up time, and thus also CO conversions (cf. Fig. 10b and c). The CO light-off leads to formation of a hot zone near monolith inlet where the reaction heat is released (cf. Fig. 10d). Fig. 11 shows the evolution of CO outlet concentration for various channels of the full monolith. It confirms that the light-off strongly depends on the channel position in the monolith, due to significant flow non-uniformities (cf. Fig. 10). 5. Conclusions A novel approach for exhaust gas aftertreatment modelling has been presented, based on a combination of standard CFD software with an independent user-defined model of monolith channel. The results of the stand-alone 1D channel model match the results of a single channel system calculated by the CFD + 1D-channels approach, which validates correct coupling of the models. Dynamic simulations of CO oxidation light-off and NSRC lean/rich operation have been then conducted in a test geometry with 10 representative channels and skewed inlet/outlet pipes, and with 80 representative channels in a complex 3D monolith systems including bent inlet and outlet popes as well as conical expansion and constriction parts. It was demonstrated that relatively large deviations can exist between the conditions in individual channels resulting from

This work has been financially supported by the Czech Science Foundation (GACR 104/08/H055), Czech Ministry of Education (MSM 6046137306) and by the specific university research (MSMT 21/2011). The authors thank to Mike Weaver (CD-Adapco) for support and useful discussions. References [1] A. Güthenke, D. Chatterjee, M. Weibel, B. Krutzsch, P. Koˇcí, M. Marek, I. Nova, E. Tronconi, Current status of modelling lean exhaust gas aftertreatment catalysts, Adv. Chem. Eng. 33 (2007) 103–211. [2] J. Braun, T. Hauber, H. Többen, J. Wundmann, P. Zacke, D. Chatterjee, C. Correa, O. Deutschmann, L. Maier, S. Tischner, J. Warnatz, Three-dimensional simulation of the transient behavior of a three-way catalytic converter, SAE Tech. Paper 2002-01-0065, 2002. [3] S.-J. Shuai, J.-X. Wang, Unsteady temperature fields of monoliths in catalytic converters, Chem. Eng. J. 100 (2004) 95–107. [4] D.N. Tsinoglou, G.C. Koltsakis, D.K. Missirlis, K.J. Yakinthos, Transient modelling of flow distribution in automotive catalytic converters, Appl. Math. Model. 28 (2004) 775–794. [5] S.F. Benjamin, Z. Liu, C.A. Roberts, Automotive catalyst design for uniform conversion efficiency, Appl. Math. Model. 28 (2004) 559–572. [6] S.F. Benjamin, C.A. Roberts, Three-dimensional modelling of NOx and particulate traps using CFD: A porous medium approach, Appl. Math. Model. 31 (2007) 2446–2460. [7] W. Guojiang, T. Song, CFD simulation of the effect of upstream flow distribution on the light-off performance of a catalytic converter, Energy Convers. Manage. 46 (2007) 2010–2031. [8] A. Kumar, S. Mazumder, Towards simulation of full-scale monolithic catalytic converters with complex heterogeneous chemistry, Comput. Chem. Eng. 34 (2010) 135–145. [9] M. Weaver, Coupling 1-D aftertreatment models with 3-D flow and solid heat conduction simulations in STAR-CD, in: DOE Crosscut Workshop on Lean Emissions Reduction Simulation, CLEERS, 20–22 April, University of Michigan, 2010. [10] StarCD User’s Guide – Methodology, CD-Adapco, 2009, pp. 15–31. [11] P. Koˇcí, M. Schejbal, T. Gregor, J. Trdliˇcka, M. Kubíˇcek, M. Marek, Transient behaviour of catalytic monolith with NOx storage capacity, Catal. Today 119 (2007) 64. ˇ epánek, S. ˇ Bártová, M. Marek, M. Kubíˇcek, V. Schmeißer, D. [12] P. Koˇcí, F. Plát, J. Stˇ Chatterjee, M. Weibel, Global kinetic model for the regeneration of NOx storage catalyst with CO, H2 and C3 H6 in the presence of CO2 and H2 O, Catal. Today 147S (2009) S257–S264. [13] K. Ramanathan, V. Balakotaiah, D.H. West, Light-off criterion and transient analysis of catalytic monoliths, Chem. Eng. Sci. 58 (2003) 1381–1405. [14] Monolith Research Group, http://www.vscht.cz/monolith, 2011.