Accelerating reactor development with accessible simulation and automated optimization tools

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

Contents lists available at ScienceDirect

Chemical Engineering & Processing: Process Intensification journal homepage: www.elsevier.com/locate/cep

Accelerating reactor development with accessible simulation and automated optimization tools

T

⁎

Eric Daymoa, , Anna Lee Tonkovicha, Matthias Hettelb, Joel Guerreroc a

Tonkomo, LLC, Gilbert, AZ, 85297, United States Institute for Chemical Technology and Polymer Chemistry, Karlsruhe Institute of Technology (KIT), Kaiserstr. 12, 76128, Karlsruhe, Germany c Wolf Dynamics SRL, Genoa, Italy b

A R T I C LE I N FO

A B S T R A C T

Keywords: CFD DAKOTA DUO OpenFOAM DETCHEM Methane CPOX Optimization Monolith catalyst

Methane CPOX is a common route to syngas, while monolith reactors are a common type of process intensified reactor. Given the general need in process intensified systems to extract improved performance, a proof-ofconcept optimization study was conducted to maximize hydrogen productivity and minimize peak reactor temperature for methane CPOX on a coated metallic monolith reactor. In this study, DAKOTA (optimization toolkit), DETCHEM™ (surface reaction and thermodynamics toolkit), and OpenFOAM (multiphysics simulation software) are newly integrated to optimize non-linear reactive configurations. Using an automated derivativefree multi-objective optimization method, the hydrogen productivity per volume (kg H2/m3/s) can be substantially increased while concurrently reducing the peak reaction temperature. Specifically, by allowing the coupled simulation framework to change the distribution and amount of catalyst, along with channel diameter, feed flow rate and inlet temperature, hydrogen productivity can be increased up to 1.8 times over previous literature values with a reduction in peak metal temperature of about 40 °C, or increased by about 1.6 times with a reduction in peak reactor metal temperature of about 120 °C, from a previously reported 818 °C down to below 700 °C. The presented results demonstrate the potential for design optimization in reactive systems in general, and in particular for process intensified components.

1. Introduction Catalytic processes are conducted in flow reactors, characterized by the spatial variation of velocity, temperature, gas-phase species concentrations, and surface coverages. The understanding and optimization of heterogeneous catalytic reactors requires detailed knowledge of reaction mechanisms and mass transport effects. Consequently, the heterogeneously catalyzed surface reactions must be analyzed together with potential homogeneous gas-phase reactions, mass transport between the active surface and the surrounding reactive flow including inside porous bodies, as well as heat transport in both the gas-phase and solid structures. Once the appropriate computational tools are developed, and a representative model of a system is constructed, then numerical simulations can be a powerful resource for analyzing the impact of various design variables on this system’s Quantities of Interest (QoI). Meanwhile, computational constraints and the availability of costeffective analysis tools can inhibit or limit the number of design options that are considered when designing reactive systems. However, modern computer equipment and the availability of low-cost or open source

⁎

automated optimization tools can lead to the discovery of new designs and/or operating modes that can reduce capital and/or operating costs, including energy, side products, emissions, and so on. Particularly for process intensification applications (e.g., a microreactor), where common fabrication techniques like etching allow for a wide range of unique design features, the opportunities for optimization are wideranging. To date, most Computational Fluid Dynamics (CFD) optimization studies have focused on external aerodynamic shape optimization problems and flows without chemistry, particularly when using the multiphysics solver OpenFOAM [1]. For example, OpenFOAM is distributed with an adjoint solver for shape optimization problems, and there is extensive published material dealing with shape optimization and design space exploration studies using OpenFOAM (e.g., [2–6]), with the objective metrics being aerodynamic quantities (i.e., lift, drag, moment, and so on). Meanwhile, chemical engineering applications have benefited from optimization processes involving CFD. For example, Kundu et al. [7] used an artificial neural network (ANN) process to quantify the stability of oil/water emulsions; Eppinger et al. [8]

Corresponding author. E-mail address: [email protected] (E. Daymo).

https://doi.org/10.1016/j.cep.2019.107582 Received 25 December 2018; Received in revised form 27 May 2019; Accepted 24 June 2019 Available online 25 June 2019 0255-2701/ © 2019 Elsevier B.V. All rights reserved.

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

applied response surface methodology (RSM) with CFD to evaluate feed composition, flow rate, temperature, catalyst mass and inert mass for optimal C2 yield for the oxidative coupling of methane (OCM) in a single-pass fixed bed reactor; Eppinger et al. [9] evaluated optimal heating temperature, feed mass flow rate, and feed composition to minimize carbon deposition on the catalyst and maximize hydrogen selectivity. Automated CFD-based catalytic reactor optimization studies (i.e., evaluation of the optimal reactor geometry, along with catalyst placement and loading, and feed composition and conditions) remain almost unexplored. Just as optimization tools are used to reduce drag, improve lift, etc. in aerodynamic applications, the same optimization methods can be applied to catalytic reactors to lower peak temperatures or improve yields, flow distribution, mass transfer, and so on. Such optimization is of significant interest to the process intensification community, particularly since the means to substantially smaller, less energy intensive, higher productivity processes requires creative design solutions which often require optimization to be commercially successful. Thus, in this work, the scale-up and optimization of a methane Catalytic Partial Oxidation (CPOX) monolith reactor is explored, with the intent to minimize peak temperature and maximize hydrogen productivity. The remainder of the manuscript is organized as follows. Section 2 describes the background information on the baseline experimental CPOX monolithic reactor design, along with a summary of prior computational work on this laboratory reactor configuration. This section also outlines the computational details for the present optimization study, followed by the test matrices for the optimization studies intended to identify design conditions that reduce peak solid temperature and increase hydrogen productivity from CPOX in a monolith reactor. Then, in Section 3, the results obtained from the optimization study are discussed. Finally, conclusions and future perspectives are presented in Section 4.

Fig. 1. Sketch of the experimental system (not drawn to scale), from Hettel et al. [11]. Reprinted from Computers & Chemical Engineering, 109, Hettel, M., Daymo, E., & Deutschmann, O., 3D modeling of a CPOX-reformer including detailed chemistry and radiation effects with DUO, 166–178, Copyright 2017, with permission from Elsevier.

flow disturbances at the catalyst inlet, the FHS was positioned 5 mm upstream of the catalyst (CAT). The back heat shield (BHS) was placed directly after the catalyst. The FHS was included in the simulations by Hettel et al. [10,11], but as will be mentioned later, this part was eliminated in the present optimization study (since the need for a heat shield diminishes as the peak temperature is lowered). For an experimental investigation of the radial heat transfer of the CAT, three single channels at different radial positions were examined. The first channel was located in the center of the catalyst (“central channel”), the second in the middle (“middle channel”) between the center and the periphery of the monolith, and the third one close to the periphery of the monolith (“outer channel”) (see Fig. 2). Measured fluid and solid temperature profiles at the central, middle, and outer positions of the FHS, CAT, and BHS are shown in Fig. 3. It is noted that the measured peak solid temperature is around 1000 K, and simulation temperatures were even higher (deviations between the experimental and simulated solid temperatures were discussed in Hettel et al. [11]). For the present study, peak simulated solid temperatures reached

2. Methods 2.1. Design basis: an unoptimized experimental CPOX process In the CPOX process, synthesis gas (a mixture of H2 and CO) is produced from a hydrocarbon fuel and oxygen in a heterogeneously catalyzed reaction. Many different processes occur simultaneously and successively on the catalytic surface. Adsorption and desorption of reactants, intermediates, and products, as well as surface reactions have a major influence on the CPOX process. Thus, the reaction network inside a CPOX reformer is composed of dozens to hundreds of elementary-step reactions. Although a detailed mechanism was used in this work, the global reaction scheme for the catalytic partial oxidation of methane is shown here,

CH 4 +

1 O2 = CO+ 2 H2 . 2

(1)

Hettel et al. [10,11] reported on 3D external and internal coupled models of the experimental reactor setup shown in Fig. 1. The rhodiumcoated monolithic honeycomb catalyst (Delphi Catalyst Inc.) was placed inside a quartz glass tube (Ø = 19 mm) with uncoated monoliths at both sides. The geometrical data of the monoliths were: Ø = 19 mm, length = 10 mm, 600 cpsi (channels per square inch), and wall thickness =80 μm. The thickness of the γ-Al2O3 washcoat on the catalyst was 40 μm and the loading was 4.23 mg/cm3. The experiments were carried out at 4 SLPM (Standard Liters Per Minute) and with a nitrogen dilution of 80 vol.%. The gases were preheated to 463 K, and the reactor was placed inside a furnace to ignite the reaction and to control the temperature (Tfurnace = 523 K). The front heat shield (FHS) and the back heat shield (BHS) (shown in Fig. 1) act as heat shields to diminish the heat loss in the up- and downstream directions due to radiation, and serve additionally as guiding devices for the suction probe and the thermocouple. To avoid

Fig. 2. Radial positions of the three channels investigated in the experiment described by Hettel et al. [11]. Reprinted from Computers & Chemical Engineering, 109, Hettel, M., Daymo, E., & Deutschmann, O., 3D modeling of a CPOX-reformer including detailed chemistry and radiation effects with DUO, 166–178, Copyright 2017, with permission from Elsevier. 2

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

are obtained with the OpenFOAM solver DUO [14]. DUO is the OpenFOAM solver which couples the two computer codes, DETCHEM and OpenFOAM, and is a synonym for the joint utilization of these two programs. OpenFOAM is an open source multi-physics solver with extensive simulation capabilities, such as incompressible flows, heat transfer, structural mechanics, chemical reactions, turbulence modeling, and so on. Being open-source, OpenFOAM offers users freedom to customize and extend its existing functionality. DETCHEM™ [14] is a commercially available package of tools specifically designed for the modeling and simulation of reacting flows based on elementary step mechanisms, particularly for heterogeneous systems such as catalysis, materials synthesis, and fuel cells. In prior work, DUO was evaluated for CPOX applications (for fixed conditions and geometry) by Hettel et al. [10,11]. As described in the work by Hettel et al. [10], a precompiled shared library (programmed in FORTRAN 90), which was made out of several parts of the DETCHEM toolbox, is coupled (linked) with the CFD software OpenFOAM. For every cell in the domain, the values of pressure, gas temperature, wall temperature and concentrations of the gas phase species are provided to the DETCHEM library by the CFD solver. The species flux terms (units of mol/m2/s) for the surface reactions occurring at walls representing catalytic coated surfaces are then calculated using the DETCHEM library at the provided conditions. These species flux terms are returned to the CFD solver and applied to the species balance equation. The overall schematic for a generic coupling between DAKOTA and OpenFOAM is shown in Fig. 4, based on guidelines for setting up an automated optimization study using DAKOTA as described in the DAKOTA User’s Manual [13] and various online resources, such as those offered by Sandia National Laboratory [15] and Wolf Dynamics [16]. In this study, DAKOTA version 6.7, OpenFOAM version 5.0, and an enhanced form of DUO found in DETCHEM version 2.5 are utilized. It is also noted that due to the rectangular nature of the monolith channel, CAD and mesh generation steps described in Fig. 4 are combined in a single step using OpenFOAM’s block-structured mesher. Additional computational details, methods, and assumptions required to complete the present CPOX optimization study are summarized below. With respect to the CFD calculations, Hettel et al. [11] studied the application of DUO for the evaluation of this CPOX system, evaluating the effect of radiation and conduction through the quartz tube surrounding the catalytic system. A multi-step surface reaction mechanism was used, which was developed for the partial oxidation and autothermal reforming of methane over rhodium [17]. The mechanism includes seven gas phase species (CH4, O2, H2, H2O, CO, CO2, N2), 12 surface species (RH, H2O, H, O, OH, CO, C, CH3, CH2, CH, CO2, CH4) and 38 reactions. The same surface mechanism is used for the present optimization work. As described in Hettel et al. [11], the reaction rates are multiplied by a factor P representing the product of the surface ratio Fcat/geo and the effectiveness factor, η :

Fig. 3. Comparison of experimentally measured solid temperature (lines) and fluid temperature (symbols) for three radial channel positions inside the FHS, CAT, and BHS from Hettel et al. [11]. Reprinted from Computers & Chemical Engineering, 109, Hettel, M., Daymo, E., & Deutschmann, O., 3D modeling of a CPOX-reformer including detailed chemistry and radiation effects with DUO, 166–178, Copyright 2017, with permission from Elsevier.

818 °C (1091 K) in the baseline simulation. In summary, the work of Hettel et al. [11] describes a reasonable experimental CPOX system. However, a significant challenge in engineering is scale-up, particularly when there are a wide range of design constraints. A design might consider a metallic monolith, prepared by roll formed substrates (e.g., such as those manufactured by Emitec [12]). A rolled FeCrAlY metallic monolith improves conduction along the reactor length. Thus, even though FeCrAlY has a higher thermal conductivity (e.g., 10.5 W/m/K) than cordierite (e.g., 2 W/m/K), it may perform more adiabatically (with less thermal spread between the outer and inner channels because of heat transfer resistance). However, peak temperatures need to be lower (e.g., possibly below 800 °C) with a metallic monolith than with a ceramic monolith. It is also desirable to lower peak temperature to increase catalyst lifetime by reducing sintering, volatility, and thermal stresses that promote catalyst delamination. To maintain or even improve hydrogen productivity (e.g., kg H2/s or kg H2/m3 reactor volume/s) while lowering the peak solid temperature, design parameters that may need to be changed include channel length, channel diameter, catalyst loading (i.e., the factor Fcat/geo , as discussed below), and inlet flow rate. As well, in the event incorporating the FHS and BHS is not feasible in a design (e.g., due to space limitations), one may wish to leave a small front segment of the catalyst uncoated, to serve as an integrated recuperator, thereby lowering peak metal and catalyst temperatures. These design parameters, and others, can be readily explored using design space exploration and design optimization, as presented hereafter.

P = Fcat/geo⋅η

(2)

Where Fcat/geo is the ratio of catalytically active surface area to geometric surface area. Mass transfer limitations inside the washcoat are taken into account by means of the effectiveness factor ηi . The effectiveness factor is the ratio of the observed reaction rate to that which would occur in case of no internal diffusion limitations. It is based on the Thiele modulus which considers washcoat thickness, concentration of the species at the boundary between fluid and washcoat, and an effective diffusion coefficient [18]. In Hettel et al. [11] the following values have been estimated based on calculations and measurements: ηi = 0.027, Fcat/geo = 185. Accordingly, the product P = Fcat/geo ∙ηi = 5 is used by the CFD solver to simultaneously account for active catalytic surface area and internal diffusion limitation in the washcoat. This product (Fcat/geo) is adjusted in the present study to represent higher and lower values of catalyst loading. No gas-phase reaction mechanism was employed in Hettel et al.

2.2. Computational setup For purposes of this study, optimization is accomplished with the Design Analysis Kit for Optimization and Terascale Applications (DAKOTA) developed by Sandia National Laboratories [13]. DAKOTA is a widely used open source optimization tool with several user-selectable optimization methods, including gradient-based methods and derivative-free methods (e.g., genetic algorithms). DAKOTA interfaces with external “black-box” software in order to compute the sensitivity of the QoI (e.g., peak solid temperature and hydrogen productivity) as a function of the vector of design parameters under evaluation (e.g., channel hydraulic diameter, catalyst loading, etc.). In this case, the QoI 3

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

Fig. 4. A typical (generic) DAKOTA-OpenFOAM optimization loop.

of the DAKOTA-CFD coupling is that pre-processing, case setup, postprocessing, and data analysis steps are implemented into one single scripting file (the “control script”) so that the optimization process is completely automated. Therefore, no human interaction is needed to recover and format the results from the CFD calculations until the final Data Analytics step (Fig. 4), and if errors should occur during any step in the optimization workflow (including the CFD calculations), the study can be stopped and restarted from a saved state. The workflow specific to this data exchange between DAKOTA and DUO is depicted in Fig. 6 and discussed below. Firstly, the DAKOTA input file is setup to reflect the number and range of design variables (e.g., inlet temperature range, inlet velocity range, etc.), the number of QoI and objective of the optimization study (minimize or maximize). The variables are given the names “x1” to “xN”, where N is the total number of design variables. The derivativefree multi-objective optimization method (MOGA, or Multi-Objective Genetic Algorithm) optimization method [13,23] is selected, along with the maximum number of CFD trials (here, at least 100). Asynchronous function evaluations are chosen so that four individual cases can run concurrently. Then, as depicted in the workflow given in Fig. 6, a “template directory” is created to store the solver input files that are parametrical, that is, subject to change as a result of the optimization process (e.g., files containing the boundary conditions for inlet velocity, inlet temperature, etc.). The automatic update of the parametrical files located in the “template directory” is achieved by using a DAKOTA supplied utility. This utility skims all files located in the “template directory” and automatically inserts the values generated by DAKOTA during design optimization or design exploration study into the predefined locations in the template files. In this workflow, a “base case directory” is also created, where all the files needed to run a DUO simulation are stored. The simulation control script file (denoted by the “Control Script” box in Fig. 6) merges the automatically edited files in the template directory with the base case directory to create a working directory for a specific set of design parameters, and then runs the OpenFOAM mesher and solver (DUO). Following the completion of the CFD simulation, standard OpenFOAM utilities are used to extract the values of the QoI (i.e., peak solid temperature and hydrogen flow rate). Since the post-processing utility generates a lot of extraneous information and further calculations are needed to evaluate hydrogen productivity, Linux shell commands and Python programs are used to isolate the necessary parameters, perform the remaining calculations to convert hydrogen flow rate (kg/s) to hydrogen productivity (kg/m3/s),

[10,11], and the present work, as it was already shown in literature that reactions in the gas phase can be neglected for the given operational conditions [19–22]. Unlike with the previous CPOX simulations described in Hettel et al. [11], for this work the radiation model is disabled, which was found in the referenced paper to have a secondary effect on the peak solid temperature at the center of the monolith, and importantly requires considerable additional computational time. Further, since an objective of the study is to lower the peak solid temperature, radiation will have a correspondingly lessened effect in the optimized system. In order to complete the optimization calculations in a timely manner, 2-D grids instead of 3-D grids are used. The 3-D quarter monolith calculations described in Hettel et al. [11] were completed using an Intel Xeon E5-2699 v3 2.3 GHz dual-CPU (36 core) server. Typical calculation times were around one week on 16 cores yielding a sum of about 2700 core hours per quarter monolith simulation. About 100,000 iterations were needed to get the final solution for each quarter monolith case, with reaction rates updated every 100-time steps. Although it is expected that computational capabilities to improve in the future, allowing more realistic simulations even during the optimization phase of the design cycle, for now, the only means to screen hundreds or thousands of potential design options with the same aforementioned hardware is to reduce the size and scope of the computational domain. The required speedup is accomplished not only by reducing the computational domain to a 2-D (symmetry) geometry, but also by eliminating the FHS or BHS, found in the experimental setup (Fig. 1), from the computational domain. For reference, the difference between a 3-D quarter monolith, 3-D quarter channel, and a 2-D channel are depicted in Fig. 5. Some loss of accuracy is expected in a 2D simulation of a rectangular channel. However, as will be discussed in Section 3, the use of a 2-D simulation for this particular case did not negatively impact the usefulness of the results. For the optimization cases, some variations to the geometry are introduced to allow for variable catalyst coating along the length of the channel (the specific geometry used for each round of optimization is discussed in Section 2.3). Also, the optimization cases presented in Section 2.3 are based on a 0.1 m-long monolith, representative of a 10x scale-up over the 0.01 m long monolith used in the experimental work. The factor of 10x is an arbitrary value but is a reasonable factor to demonstrate the present optimization workflow for scale-up analyses. With these changes, the run time for a 0.1 m long channel is about 2 h. Meanwhile, DAKOTA setup is fairly straightforward. A key feature

4

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

Fig. 5. Computational domains for a Quarter Monolith [10], Quarter Channel [10], and the 2-D channel from the present work. The solid material is gray, and the fluid region is dark (blue). The 3-D quarter channel domain is the center of the 3D quarter monolith, denoted by the red circle in the left figure. The 2-D channel is comprised of the center plane of the quarter channel, denoted by the dashed black line in the right figure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

and format the necessary results into a flat text file suitable for DAKOTA. If more advanced post-processing is needed to process CFD output, ParaView’s [24] Python interface can be utilized to perform advanced calculations and manipulations (e.g., custom “slices” through the calculation domain). Within DAKOTA, the MOGA method reads the aforementioned results file and ranks the degree of fit – the closeness to the “best solution” – according to the optimization objectives (i.e., lower peak solid temperature and higher hydrogen productivity). The genetic algorithm selects new conditions by performing a crossover step where traits of the best design parameters (e.g., catalyst loading, inlet temperature, channel diameter, etc., according to the test matrices shown in Table 1) are used as the starting point for new test conditions. In order to promote global optimization and prevent settling on a local optimum, “mutations” are introduced whereby the best design parameters found in prior runs can be adjusted to significantly different values, in order to ensure the design space is satisfactorily explored. With the experimental basis and the computational method established, it is now possible to conduct the optimization study with the goals of maximizing the hydrogen productivity and minimizing the peak solid temperature of a scaled-up methane CPOX monolithic reactor.

account for the difference in channel volume as the channel length is extended from 0.01 m to 0.1 m. For purposes of this study, in all cases the monolith length is fixed to 0.1 m, and the inlet feed gas is 80% N2 with a carbon/oxygen (C/O) ratio of 1.0, as in the system studied by Hettel et al. [11]. It is noted that in industrial systems that the diluted O2/N2 mixture will be replaced by pure O2 or air, but in order to improve comparisons with the unoptimized experimental system, for purposes of this comparison study the O2/N2 ratio is unchanged from the experimental system. In the first set of optimization tests, the channel is divided into two regions, one with catalyst and the other without catalyst (as described in Section 2.3.1). In the second set of optimization tests, the channel is divided into four regions, one without catalyst, and three regions with independently varying catalyst loadings (as described in Section 2.3.2). 2.3.1. One catalyst coated region test matrix (optimization rounds 1 and 2) In the first set of design optimization studies (Table 1), multiple variables were adjusted by DAKOTA, as listed below:

• inlet velocity (for Rounds 2A to 2D, but not Round 1), • inlet temperature, • channel diameter (for Rounds 2A to 2D, but not Round 1). Note that

2.3. Optimization study

•

Several rounds of optimization cases were conducted to identify conditions that can improve the hydrogen productivity (either in terms of kg H2/s or kg H2/m3/s while minimizing peak solid temperature. As will be discussed further in Section 3, even if optimization is performed on a kg H2/s basis, for results comparison to the laboratory-scale system, ultimately kg H2/m3/s is the preferred metric for comparison between the optimized and baseline laboratory-scale system, in order to

•

only the diameter of the fluid region is adjusted; the solid wall width is constant, percent of the 0.1 m channel wall that remains uncoated (no catalyst activity in the uncoated region, so the catalyst wall acts as a built-in recuperative heat exchanger), and Fcat/geo (the ratio of active catalyst area to geometric area, representing catalyst loading).

Note that for clarity, a change in the round number (1 vs. 2 vs. 3) denotes a change to the design parameter matrix. Lettering denotes a

Fig. 6. Workflow for data exchange between DAKOTA and DUO. White denotes process blocks, dark shaded (blue) document symbols denote unchanging sets of files, and light shaded (green) document symbols denote files that change with each set of design parameters specified by DAKOTA. The grey shaded area denotes the domain of the control script that automatically prepares the OpenFOAM case, runs the DUO solver, post-processes the result, and prepares the results file for DAKOTA. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 5

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

Table 1 Test matrix for Rounds 1, 2A, 2B, 2C, and 2D. Design parameters are the inlet velocity, inlet temperature, channel diameter, % uncoated, and catalyst loading (Fcat/geo ). Round

U inlet Range (m/s)

Inlet T Range (°C)

D channel (mm)

% uncoated

Fcat/geo

Obj Fn Notes

1 2A 2B 2C 2D

3.29 fixed 3.4–5.1 3.4–5.1 3.4–5.1 3.4–5.1

100–500 100–400 100–400 100–400 100–400

0.8765 fixed 0.8–3.2 0.8–3.2 0.8–3.2 0.8–3.2

0–50 0–50 0–50 0–50 0–50

1–10 3–30 3–30 3–30 3–30

Equal (kg Equal (kg Equal (kg 80%-20% 20%-80%

H2/s) H2/s) H2/m3/s) (kg H2/m3/s) (kg H2/m3/s)

# Runs 100 100 100 100 100

approach where spray coating may be used to apply the requisite catalyst loadings to different segments of a reactor (e.g., a laminate construction method, as seen in some microchannel reactors, to ensure the catalyst can be applied before the laminates are joined to form the completed reactor structure). In addition, in Round 3A the productivity objective function is the same as that for 2B (i.e., DAKOTA’s multi-objective optimization method minimizes peak solid temperature and maximizes kg H2/m3/s, with equal weighting for the hydrogen productivity and thermal QoI). The results from round 3A showed no improvement over Round 2D, so for Round 3B the CFD solver productivity QoI is modified. Instead of optimizing around the maximum H2 productivity, the QoI is now based on the “closeness” to the curve of productivity vs. peak solid temperature according to the curve fit (solid line) shown in Fig. 9, as described by the pseudocode in Fig. 10. To reiterate, the optimization objective to minimize peak solid temperature is maintained for Round 3B, but productivity is more heavily weighted the further it exceeds a function that describes the best results from Round 2D. The idea of these response metrics (QoI) is to more heavily weight responses which are further above the fit line (i.e., better than Round 2D), and more negatively weight responses based on their distance below the fit line (same or worse than Round 2D).

change to the objective function weighting for a given set of design parameters (e.g., 2A, 2B, 2C and 2D have the same design parameters delineated in Table 1, and the lettering notes changes to the objective function). The QoI are hydrogen productivity, either expressed as kg H2/s (Rounds 1 and 2A), or kg H2/m3/s (Rounds 2B, 2C and 2D), and maximum solid temperature in degrees Celsius. In Rounds 1, 2A and 2B the QoI are given equal weight. In Round 2C, kg H2/m3/s is weighted 80% and the maximum solid temperature is weighted 20%. In Round 2D, the weightings are reversed, so the kg H2/m3/s is weighted 20% and the maximum solid temperature is weighted 80%. The geometry simulated in Rounds 1 and 2 is depicted in Fig. 7. 2.3.2. Multiple catalyst coated region test matrix (optimization round 3) Given the goal of pushing the limits of improved productivity at minimal peak solid temperature, Round 3 cases (Table 2) were developed to have greater flexibility in the catalyst loading distribution (note that in Round 3 cases, the fluid channel diameter is fixed at 0.81 mm, the approximate optimal value from Rounds 1 through 2D). Specifically, the solid catalytic wall is split into 4 zones. Zone 1 is uncoated and can range in length between 1 and 49 mm. Zone 2 is coated with catalyst and has a length of 50 mm minus the length of Zone 1. Zones 3 and Zones 4 are both 25 mm in length (the total length of Zones 1 to 4 = 100 mm for same total length in the previous rounds of optimization). The catalyst loading in Zones 2, 3, and 4 can be adjusted separately by DAKOTA’s multi-objective optimization method, according to the test matrix outlined in Table 2. The geometry considered for Round 3 of optimization testing is graphically shown in Fig. 8. Increasing discretization and complexity in catalyst location was hypothesized to further improve the productivity and minimize the peak catalyst temperature. It is recognized that multiple catalyst loadings are often not applied to a monolith catalyst, but so long as the loading increases along the length of the catalyst (i.e., Fcat/geo zone 1 < Fcat/geo zone 2 < Fcat/geo zone 3 < Fcat/geo zone 4) it should be feasible to dip coat additional catalyst as a function of channel length. In this study, a constraint was not added so that Fcat/geo increases in each catalyst zone, so it is possible to realize solutions which would be difficult to dip coat (and this is seen in some of the best results in Table 4). In the event there truly were a significantly better solution with catalyst loadings following a trend that is not easily applied with standard dip coating methods, the product may benefit from a redesign to an

3. Results and discussion The consolidated results of the study are shown in Fig. 11, with representative best results for Rounds 2C and 3B provided in Tables 3 and 4, respectively. As noted, the calculation workflow optimizes on either kg H2/s or kg H2/m3/s, but for data analysis purposes the productivity results are firstly converted to kg H2/m3/s if necessary (i.e., only for Rounds 1 and 2A), and then normalized by the predicted hydrogen productivity (in the same units of kg H2/m3/s) for the conditions and geometry analyzed with DUO and described in the paper by Hettel et al. [11]. It is immediately visible in Fig. 11 that a relatively simple optimization with DUO and DAKOTA yields a design option that has a hydrogen productivity of about 1.6 times the baseline productivity, at around a 120 °C lower peak solid temperature from 818 °C to 700 °C (noting that in the baseline case derived from Hettel et al. [11], a ceramic monolith is used, but for the present study it is assumed a

Fig. 7. Overall geometry used for the Round 1 and Rounds 2 channel optimization study, not to scale. Luncoated = length of the uncoated (no catalyst) region, Dfluid = hydraulic diameter of the fluid channel, U = inlet velocity, T = inlet fluid temperature (T). The solid material is FeCrAlY. 6

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

Table 2 Test matrix for Rounds 3A and 3B. Round

U Range (m/ s)

Inlet T Range (°C)

Zone 1 uncoated length (mm)

Zone 2 uncoated length (mm)

Zone 3 uncoated length (mm)

Zone 4 uncoated length (mm)

Objective Function Method

# Runs

3A 3B

3.3–5.1 3.3–5.1

100–400 100–400

1-4.9 1-49

3–30 3–30

15–30 15–30

15–30 15–30

Same Rd2B-D Fit curve

200 500

FeCrAlY monolith is specified). Alternatively, the maximum productivity can be increased by about 1.8 times the baseline productivity with a reduction in the peak temperature of about 40 °C. Assuming a lower peak solid temperature is more desirable than a slightly higher productivity, then for an industrial application, one could consider the following as a starting point for a more detailed design: Inlet T ≈118 °C, inlet U ≈4.9 m/s, ≈23% of the 0.1 m long channel remains uncoated, the remaining channel is coated with a catalyst with an Fcat/geo ≈27, and channel diameter ≈0.81 mm. It is also immediately noticeable in Fig. 11 that optimizing on hydrogen productivity per reactor volume (kg H2/m3/s) is preferable to optimizing on the mass flow rate of kg H2/s at the reactor outlet, because including the reactor volume in the objective function accentuates the benefits of smaller channel diameter associated with process intensification (e.g., improved mass transfer). Thus, carefully crafting the objective functions can improve the search for optimal design parameters. Interestingly, in Round 3B, which includes multiple zones with varying catalyst loading and an objective function that penalizes performance below that seen in Round 2D, did not yield significant improvements in performance over Round 2 cases. That is, the simplified binary catalyst coating scheme explored in the second round of optimization provided around the same level of performance enhancement as the more complicated scheme explored in the third round. This observation is based on a fixed channel length of 10 cm (an arbitrary length for purposes of this study, but representative of a 10x scale-up from the laboratory experiments). It is possible that as the channel length increases (and correspondingly, as the flow rate increases to maintain the same weight hourly space velocity, WHSV) during scaleup for industrial reactor design, greater discretization of catalyst activity might provide greater control of the peak temperature. Also, comparing the relative loadings for Rounds 2C and 3B (Tables 3 and 4), the 263.6 °C inlet case has a slightly better performance at a lower overall amount of applied catalyst in Round 3B, suggesting that if the optimization objective function were to include catalyst cost, then the optimal configuration may involve multiple zones with different catalyst loadings. However, as seen here, given the present objective functions, sometimes a more complex solution does not yield a better result. Finally, it is noted that given the number of design variables and the range of these variables, the number of runs is probably insufficient to find the true globally optimum design values. Therefore, a large number of calculations were run (in total over 1000 runs) to ensure that a

Fig. 9. The productivity objective function for Round 3B is adjusted so that it is related to the curve fit of relative hydrogen productivity vs. peak solid temperature, according to the SAS JMP analysis (Logistic 3 P model curve fit) shown here.

Fig. 10. Pseudocode for Round 3B, which heavily weights solutions that are better than the best solutions identified in Round 2D (i.e., above the fitted line shown in Fig. 9), and penalizes solutions that are below this fitted line.

narrowly defined local optimum was not missed. However, immediate improvement over the baseline experimental condition is readily seen with just around 100 runs. That is, the approximate optimum results from the study can be approximately identified in any of the 100-run test matrices 2B-2D, as shown in Fig. 11. For example, in Fig. 11 Round

Fig. 8. Overall geometry used for the Round 3A-3B channel optimization studies, not to scale. The key design parameters are clearly outlined: inlet velocity (U), inlet fluid temperature (T), length of the uncoated (no catalyst) region, and the Fcat/geo of the catalyst zones 2, 3 and 4. The solid material is FeCrAlY. 7

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

Fig. 11. Results of the productivity ratio versus peak temperature for all rounds of optimization. “Hand-drawn” curves were added to help guide the eye while interpreting improvement over major rounds of optimization (blue for Round 1, red for Round 2A, and black for Rounds 2B, 2C, 2D, 3A and 3B). A relative value of 1.0 is assigned to the kg H2/m3/s productivity for the conditions based on Hettel et al. [11] with a 0.01 m-long monolith. The right-most (red) circle denotes the initial laboratory result, the center (blue) circle denotes a 1.8 times productivity improvement with a 40 °C decrease in peak solid temperature, and the left-most (green) circle denotes a 1.6 times improvement in hydrogen productivity at a 120 °C lower peak solid temperature. The use of the improved optimization function in Round 3B provided a few marginally improved results, primarily at high temperatures (exceeding 750 °C) and near the inflection point where the productivity ratio rises sharply with temperature between 690 and 710 °C. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Since the 2-D calculation is conservative compared to the 3-D case, the translation from 3-D to 2-D was quite reasonable for reducing the computational cost of the present optimization study since the 2-D results do not overpredict the benefit of proposed design changes at the optimal design condition that yields a peak solid temperature of 700 °C. With further simulation and experimental activities guided by the optimization study, ideally the final product design will operate at much lower temperature and with substantially improved productivity as compared to the laboratory system described in Hettel et al. [11].

Table 3 Representative best results from optimization Round 2C. Inlet T °C

Ch Diam (mm)

Frac Uncoated

Fcat/geo

Inlet U (m/s)

Relative kg H2/m3/s

Peak Solid T °C

330.6 263.6 118.6

0.81 0.81 0.81

0.34 0.23 0.23

26.9 26.9 26.9

4.9 4.9 4.9

2.43 2.30 2.15

775.6 723.4 700.1

2D results (a 100-run test matrix) are within the blue and green circles, denoting the preferred solutions from this study. It is further noted that there is a steep drop-off in productivity below a peak solid T of about 700 °C. On one hand, this is not unexpected because CPOX has a temperature floor below which reaction light-off is unachievable. Also, due to mass transfer resistances, there will be some limit to the benefit of increasing catalyst loading (Fcat/geo ). With the optimization study completed, 3-D simulations can be performed, and experimental data can be collected to verify the results of the simulation studies. For example, for the starting point design conditions described above were evaluated with a 3-D quarter channel grid. As seen in Table 5, the 2-D case performs slightly worse than the 3D case (in part because the 2-D grid, depicted in Fig. 5, is unable to capture the catalyst on both the sides and top/bottom of the channel).

4. Conclusions The availability of computational resources and low-cost or open source software to conduct the CFD and optimization studies will surely unlock new, unforeseen, design options in the fields of reaction engineering and process intensification. From this relatively simple study alone, for example, coupling the optimization toolkit DAKOTA with the heterogeneous OpenFOAM-based CFD solver DUO, led to a greater than 1.6 times improvement in hydrogen productivity at around a 120 °C lower peak solid temperature (which represents a reduction in peak temperature from 818 °C to just under 700 °C) or a 1.8 times productivity improvement with a 40 °C lower wall temperature. This improvement as accomplished by adjusting the channel diameter, allowing a length of uncoated catalyst at the front of the monolith to pre-

Table 4 Representative best results from optimization Round 3B. Inlet T °C

inlet U (m/s)

Zone 1 L (mm)

Zone 2 L (mm)

Fcat/geo Zone 2

Fcat/geo Zone 3

Fcat/geo Zone 4

Relative kg H2/m3/s

Peak Solid T (°C)

330.6 263.6 118.6

5.1 4.9 4.9

16 27 23

34 23 27

23.1 22.5 27.3

28.6 24.1 28.7

28.3 28.3 26.8

2.54 2.37 2.15

781.2 757.0 699.6

8

Chemical Engineering & Processing: Process Intensification 142 (2019) 107582

E. Daymo, et al.

Table 5 Comparison of 3-D and 2-D grid results from Round 2C at the proposed optimized design conditions (high productivity at a peak solid temperature of 700 °C). 2D or 3D

Inlet T (°C)

Ch Diam (mm)

Frac Uncoated

Fcat/geo

Inlet U (m/s)

Relative kg H2/m3/s

Peak Solid T (°C)

2D 3D

118.6 118.6

0.81 0.81

0.23 0.23

26.9 26.9

4.9 4.9

2.15 2.75

700.1 694.8

heat the incoming feed gases, and adjusting the inlet velocity and temperature to more optimal values. At this stage, optimization was not performed on the inlet feed gas composition (e.g., level of nitrogen dilution) and total reactor length. Also, this methodology, when coupled with sufficient computational resources, could be used to evaluate more realistic geometries (e.g., a quarter monolith, including a header and footer), or even the inclusion of gas phase mechanisms. Further, such optimization could be applied to parameter fitting for finding kinetic rate parameters (e.g., varying the pre-exponential factor and the activation energy in order to minimize the sum-of-squares differences between experimentally measured conversion and selectivity and the corresponding model predictions). The opportunity to improve designs is practically boundless when the appropriate, cost-effective reactive CFD and optimization tools are available to engineers.

[8]

[9]

[10]

[11]

[12] [13]

Acknowledgments A cost-free academic license of DETCHEM™ by the Steinbeis GmbH & Co. KG für Technologietransfer (STZ 240 Reaktive Strömung) is gratefully acknowledged. Tonkomo, LLC self-funded the work in this paper as its own research and development, with the exception of a cost-free academic license of DETCHEM.

[14]

[15] [16] [17]

References [1] H.G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys. 12 (1998) 620–631. [2] J. Byrne, P. Cardiff, A. Brabazon, M. O’Neill, Evolving parametric aircraft models for design exploration and optimization, Neurocomputing 142 (2014) 39–47. [3] J. Guerrero, A. Cominetti, J. Pralits, D. Villa, Surrogate-based optimization using an open-source framework: the bulbous bow shape optimization case, Math. Comput. Appl. 23 (2018) 60. [4] P. He, C.A. Mader, J.R.R.A. Martins, K. Maki, An aerodynamic design optimization framework using a discrete adjoint approach with OpenFOAM, Comput. Fluids 168 (2018) 285–303. [5] A. Ohm, H.O..̈ Pétursson, Automated CFD Optimisation of a Small Hydro Turbine for Water Distribution Networks, Master’s thesis in Applied Mechanics Chalmers University of Technology, Gothenburg, Sweden, 2017. [6] N. Umetani, B. Bickel, Learning three-dimensional flow for interactive aerodynamic design, ACM Trans. Graph. 37 (2018) 89. [7] P. Kundu, V. Paul, V. Kumar, I.M. Mishra, Formulation development, modeling and

[18] [19]

[20]

[21] [22] [23] [24]

9

optimization of emulsification process using evolving RSM coupled hybrid ANN-GA framework, Chem. Eng. Res. Des. 104 (2015) 773–790. T. Eppinger, G. Wehinger, M. Kraume, Parameter optimization for the oxidative coupling of methane in a fixed bed reactor by combination of response surface methodology and computational fluid dynamics, Chem. Eng. Res. Des. 92 (2014) 1693–1703. T. Eppinger, G.D. Wehinger, N. Jurtz, R. Aglave, M. Kraume, A numerical optimization study on the catalytic dry reforming of methane in a spatially resolved fixedbed reactor, Chem. Eng. Res. Des. 115 (2016) 374–381. M. Hettel, C. Diehm, H. Bonart, O. Deutschmann, Numerical simulation of a structured catalytic methane reformer by DUO: the new computational interface for OpenFOAM® and DETCHEM™, Catal. Today 258 (2015) 230–240. M. Hettel, E. Daymo, O. Deutschmann, 3D modeling of a CPOX-reformer including detailed chemistry and radiation effects with DUO, Comput. Chem. Eng. 109 (2018) 166–178. Emitec, URL: www.emitec.com. (Accessed 16 November 2018). B.M. Adams, W.J. Bohnhoff, K.R. Dalbey, M.S. Ebeida, J.P. Eddy, M.S. Eldred, J.R. Frye, G. Geraci, R.W. Hooper, P.D. Hough, K.T. Hu, J.D. Jakeman, M. Khalil, K.A. Maupin, J.A. Monschke, E.M. Ridgway, A. Rushdi, J.A. Stephens, L.P. Swiler, D.M. Vigil, T.M. Wildey, J.G. Winokur, A. Dakota, Multilevel Parallel ObjectOriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 6.7 User’s Manual, Sandia Technical Report SAND2014-4633, Updated November 13, 2017 (2017). O. Deutschmann, S. Tischer, C. Correa, D. Chatterjee, S. Kleditzsch, V.M. Janardhanan, N. Mladenov, H.D. Minh, H. Karadeniz, M. Hettel, DETCHEM Software Package, Karlsruhe, Germany: (2014) www.detchem.com. Sandia National Laboratory. URL: https://dakota.sandia.gov/training. (Accessed 13 May 2019). Wolf Dynamics. URL: http://www.wolfdynamics.com/tutorials.html?id=46. (Accessed 13 May 2019). O. Deutschmann, R. Schwiedernoch, L.I. Maier, D. Chatterjee, Natural gas conversion in monolithic catalysts: interaction of chemical reactions and transport phenomena, Stud. Surf. Sci. Catal. 136 (2001) 251–258. O. Deutschmann, Modeling and Simulation of Heterogeneous Catalytic Reactions: From the Molecular Process to the Technical System, Wiley-VCH, Weinheim, 2011. A. Beretta, A. Donazzi, D. Livio, M. Maestri, G. Groppi, E. Tronconi, P. Forzatti, Optimal design of a CH4 CPO-reformer with honeycomb catalyst: combined effect of catalyst load and channel size on the surface temperature profile, Catal. Today 171 (2011) 79–83. A. Bitsch-Larsen, R. Horn, L.D. Schmidt, Catalytic partial oxidation of methane on rhodium and platinum: spatial profiles at elevated pressure, Appl. Catal. A Gen. 348 (2008) 165–172. O. Deutschmann, L.D. Schmidt, Modeling the partial oxidation of methane in a short-contact-time reactor, AICHE J. 44 (1998) 2465–2477. G. Veser, J. Frauhammer, Modelling steady state and ignition during catalytic methane oxidation in a monolith reactor, Chem. Eng. Sci. 55 (2000) 2271–2286. C.C. Coello, G.B. Lamont, D.A. Van Veldhuizen, Evolutionary Algorithms for Solving Multi-Objective Problems, 2nd edition, Springer, New York, 2007. ParaView. URL: www.paraview.org. (Accessed 21 May 2019).