OpenFOAM Computational Fluid Dynamic Simulations of Single

Article pubs.acs.org/IECR

OpenFOAM Computational Fluid Dynamic Simulations of SinglePhase Flows in an Advanced-Flow Reactor María José Nieves-Remacha,† Amol A. Kulkarni,†,‡ and Klavs F. Jensen*,† †

Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States CEPD, CSIR-National Chemical Laboratory, Pune 411008, India

‡

Downloaded by STOCKHOLM UNIV on August 27, 2015 | http://pubs.acs.org Publication Date (Web): July 27, 2015 | doi: 10.1021/acs.iecr.5b00232

S Supporting Information *

ABSTRACT: Computational fluid dynamic (CFD) simulations are carried out for single-phase flow in an advanced-flow reactor (AFR) using the open source software OpenFOAM. Excellent agreement of simulations with experimental pressure drop and residence time distribution (RTD) is obtained. Streamlines, stagnant zones, velocity profiles, and pressure fields are obtained at different flow rates ranging from 5 to 100 mL/min. A change in the flow behavior with the presence of recirculation zones is observed with a 40 mL/min flow rate. The extent of the recirculation zones increases with increasing flow rate from 40 to 60 mL/min and is limited further by the presence of a second cylindrical post inside the heart cell, remaining almost constant in the flow rate range of 60−100 mL/min. The RTD is also determined for all flow rates, and a comparison between different reactor designs (two-post, single-post, and low-flow-reactor-like single-post) is presented. The AFR shows a plug-flow behavior with a small degree of dispersion, which broadens the RTD. Symmetric RTD curves are obtained for the single-post designs, whereas the Gen 1 AFR design experiences asymmetry in the RTD at flow rates in the range between 20 and 60 mL/min.

■

INTRODUCTION Microreactors have been demonstrated to provide advantages over conventional process technologies for the synthesis of chemical compounds and kinetic studies at laboratory scale.1 High heat and mass transfer rates, rapid mixing, and higher selectivities and conversions can be achieved in these microdevices thanks to the small characteristic dimensions, enabling the synthesis of compounds that cannot be synthesized in conventional reactors. In the past years, efforts have been directed toward the application of microreactor technology for production purposes, especially in the pharmaceutical and fine chemicals industry.2−4 The challenge is how to get the benefit of the transport rates inherent to microreactors while increasing the throughput for production applications. Two approaches to increase production rate are possible: (a) scale-out by parallelization of units and (b) scaleup by increase in channel size and flow rates. Scale-out requires thousands of units to achieve kilogram per minute of production rates and development of very expensive and complex control systems to ensure identical operating conditions in each unit for a perfect and predictable overall reactor performance. In contrast, scale-up by increase in channel size risks losing mass and heat transfer performance. The advanced-flow reactor (AFR) manufactured by Corning, Inc., combines both approaches and is able to yield production rates of 10−300 g/min per module.5 The AFR is formed by three layers, and the reaction mixture flows through the middle one. The reaction layer has a special design in the form of several rows of heart-shaped cells in series (Figure 1a) that form a sequence of convergent−divergent segments that enhance mass transfer rates. The height for this layer is 1.1 mm, and the minimum width is 1 mm (Figure 1b). The two outer layers are devoted to the flow of the heating or cooling fluid (Figure 1c). Although each module has a volume of 8.7 © 2015 American Chemical Society

mL, parallelization of several units helps to achieve production rates on the order of a kilogram per minute (Figure 1d). If the AFR is demonstrated to perform efficiently and to be easily scalable, then it could become an alternative for process intensification and transition from batch to continuous production in the pharmaceutical and fine chemicals industry. Additional advantages include shorter process development times thanks to the scalability of the reactor modules, higher selectivities and yields, greener production processes, and the possibility of introducing new chemistries.6−10 Knowledge of the hydrodynamic characteristics of a reactor is essential for reactor design, development, and optimization. The traditional approach for process development is experimentally based. However, experiments provide limited information about the behavior inside the reactor. Computational fluid dynamics (CFD) is a useful tool that allows the prediction of the hydrodynamics without the need to perform experiments once the model has been validated. In particular, the complex geometry of the AFR and its design makes it difficult to obtain information locally from experiments. The availability of local information helps one to understand how the system works and to make decisions about how to modify the reactor design to improve its performance by overcoming the limitations of the existing devices. Moreover, CFD simulations can be used to study the effect of the fluid properties on the flow characteristics and its influence on reactor performance. In this way, the CFD tool can be further used to optimize reaction conditions and maximize yields, minimize energy requirements, and identify safety concerns. Received: Revised: Accepted: Published: 7543

January 17, 2015 June 21, 2015 July 7, 2015 July 27, 2015 DOI: 10.1021/acs.iecr.5b00232 Ind. Eng. Chem. Res. 2015, 54, 7543−7553

Article

Industrial & Engineering Chemistry Research


where the overbar indicates time-averaged values and v′ is the fluctuating velocity. Turbulence models allow for the calculation of the mean flow without the need to calculate the full time-dependent flow field. We refer the reader to the literature on turbulent models for more details.13−18 The CFD simulations shown here are carried out at steady-state for single-phase flow at flow rates ranging from 5 to 100 mL/min. Within this interval, only laminar flow is required to accurately describe the system. The AFR is able to operate at flow rates up to 300 mL/min, so turbulence models may be needed at these flow rates. Determination of the RTD is implemented using a convection−diffusion equation for the tracer (CT = tracer concentration), which is treated as a passive scalar on the basis of the velocity profiles obtained at steady-state (v)⃗ . The transport equation is shown in eq 5, where DT is the diffusion coefficient of the tracer and is assumed to be constant. ∂C T 2 + ∇·(vC ⃗ T) − ∇ (DTC T) = 0 ∂t

Numerical Approach. OpenFOAM is based on a finite volume discretization method applied to arbitrary-shaped cells and uses a segregated iterative solution.6 The solver used for our simulations to study the flow in the AFR is simpleFoam, implemented for incompressible steady-state single-phase flow. Once the flow field is known, the RTD is determined using the scalarTransportFoam solver that solves a transport equation for a passive scalar, which in this case is the tracer injected at the reactor inlet. To solve the equations for velocity and pressure, these solvers use the pressure-implicit with splitting of operators (PISO) algorithm for transient problems or semiimplicit method for pressure-linked equations (SIMPLE) algorithm for steady-state problems. Both algorithms in OpenFOAM can introduce nonorthogonality correctors (0 for an orthogonal mesh and 20 for the most nonorthogonal mesh).6,19 Velocity and concentrations are solved with preconditioned biconjugate gradient (PBiCG) and the preconditioner diagonal-incomplete LU (DILU). For pressure, a preconditioned conjugate gradient (PCG) solver with the preconditioner diagonal-incomplete Cholesky (DIC) is used. Relaxation factors are used to stabilize the computation, especially for steady-state problems. They work by limiting the amount by which the variable changes between iteration steps. Small values enhance stability but make the computation slower. Maximum and minimum values for relaxation factors are 0.2 and 0.9. On the basis of CFD experience, relaxation factors have been set to 0.3 for the pressure field and 0.7 for the velocity. The numerical schemes used in this work are Gauss linear scheme for gradient, Gauss upwind/linear/linearUpwind divergence calculations (upwind and linear are first- and second-order accuracy, respectively), and Gauss linear corrected for Laplacian. For steady-state simulations, the time derivative is not solved and the keyword steadystate within the ddtSchemes is used. Upwind (first-order) schemes usually provide physical results, but they have the problem of creating instabilities and numerical diffusion, especially when the flow is not parallel to the grid. For the scalarTransportFoam simulations, Gauss linear limited was used in order to reduce numerical diffusion.6 Mesh. The computational domain for the entire AFR is discretized in a 3D-structured mesh composed of 378 460 hexahedral cells as shown in Figure 2. The mesh is built using the software Pointwise V.17. The structured mesh provides

Figure 1. Details of the structure of the advanced-flow reactor (AFR) Gen 1: (a) front view of the AFR mixing module, (b) detail of a heartshaped cell, (c) lateral cross section of a single AFR module showing reaction and heat-transfer zones, and (d) typical AFR system with multiple connected plates.5

CFD simulations are carried out here to simulate single-phase incompressible flow through the AFR, using the open source software OpenFOAM (field operation and manipulation).11 The response variables analyzed are pressure drop, velocity profiles, streamlines, and swirling strength, all as a function of flow rate. For solver validation, CFD results are compared with experimental measurements of pressure drop and residence time distributions (RTDs). Although not included here, further validation of CFD velocity profiles could be carried out by microparticle image velocity (micro-PIV) measurements.12 This would probably require a special AFR module with the two outer heating/cooling layers removed. Basic Physical Equations. The governing equations for single-phase flow of an incompressible Newtonian fluid under isothermal conditions with only gravity as body force are given by the continuity equation (eq 1) and the conservation of linear momentum equation (eq 2). ∇·v ⃗ = 0

(1)

⎛ ∂v ⃗ ⎞ ⎯→ ⎯ ρ⎜ + v ⃗ ·∇v ⃗⎟ = −∇P + ρg ⃗ + μ∇2 v ⃗ ⎝ ∂t ⎠

(2)

If turbulence needs to be considered, then additional equations for the Reynolds stress or turbulent stress (v′⃗ v′⃗ ) that provide closure to the Reynolds-averaged Navier−Stokes equations are required. ∇·v ⃗ = 0

ρ

(5)

(3)

⎛ ∂v ⃗ ⎞ ⎯⎯→ Dv ⃗ = ρ⎜ + v ⃗ ·∇v ⃗ ⎟ = −∇P + ρg ⃗ + μ∇2 v ⃗ + v ⃗′v ′⃗ ⎝ ∂t ⎠ Dt (4) 7544

DOI: 10.1021/acs.iecr.5b00232 Ind. Eng. Chem. Res. 2015, 54, 7543−7553

Article

Industrial & Engineering Chemistry Research Table 1a. Numerical Methods Used in OpenFOAM Simulations solver

operator

simpleFoam

solver

method

gradient divergence

scalarTransportFoam

Table 2. Boundary Conditions for OpenFOAM Simulations

Laplacian time derivative gradient divergence Laplacian time derivative

Gauss linear Gauss upwind/linear/ linearUpwind Gauss linear corrected steadystate Gauss linear Gauss limited linear Gauss linear corrected Euler

velocity simpleFoam pressure

scalarTransportFoam


pressure velocity concentration

solver preconditioned conjugate gradient preconditioner diagonal-incomplete Cholesky preconditioned biconjugate gradient preconditioner diagonal-incomplete LU

concentration

boundary location

boundary condition

inlet outlet walls inlet outlet walls inlet outlet walls

fixed flat profile zero gradient zero velocity zero gradient zero pressure zero gradient fixed scalar zero gradient zero gradient

profiles, streamlines, stagnant zones, and pressure drop are shown in the Results and Discussion section. Pressure drop for the entire AFR is computed for different flow rates and compared with experimental data. The swirling strength is computed to identify stagnant zones. This is quantified by the imaginary part of the complex eigenvalue of the local velocity gradient tensor. Vortices are identified by plotting isoregions of imaginary part of the complex conjugate. Because the computation of the swirling strength is carried out in the center plane of the reactor (parallel to the z axis), we can assume symmetry of the heart cell in the z direction. This allows for simplification of the calculation of the velocity gradient tensor D as represented by eq 6.20

Table 1b. Solvers Used in OpenFOAM Simulations variables

variable

abbreviation PCG DIC PBiCG DILU

D=

∂vx ∂x ∂x1 ∂x 2 ∂vy ∂x1

Figure 2. Structured 3D mesh composed by 378 460 cells for the entire advanced-flow reactor.

∂y ∂x 2

(6)

Velocity profiles are used to obtain the RTD by response to a step size in concentration at the inlet. The monitored concentration at the outlet is differentiated with respect to time to obtain E(t).

■

more efficiency, faster convergence, and better resolution as compared to unstructured mesh. Note that the small characteristic dimensions of the AFR make it necessary to use a 3D mesh versus a 2D mesh in order to accurately represent the system (comparison in the Supporting Information). Structured meshes show a better numerical performance than unstructured meshes, especially when the flow is aligned with the cells. Numerical diffusion occurs more significantly when there is an angle between the flow direction and the mesh grid. When the flow fields are parallel to the grid, the numerical diffusion is minimized. Our mesh was carefully built with hexahedral cells keeping the aspect ratio and orthogonality and reducing skewness whenever possible. Boundary Conditions. All simulations share the following boundary conditions as shown in Table 2: the velocity is fixed at the reactor inlet with a flat profile, nonslip condition at the walls, and zero gradient velocity at the outlet. The pressure is fixed at the outlet (zero, reactor outlet opened to the atmosphere), whereas the gradient pressure at the inlet and at the walls is zero. For RTD simulations, the concentration at the inlet is well mixed and constant, the flux at the walls is zero, and the gradient concentration at the outlet of the reactor is zero. Analysis of Simulations. Steady-state simulations for single-phase flow using water as incompressible fluid for flow rates ranging from 5 to 100 mL/min are carried out. Velocity

RESULTS AND DISCUSSION Postprocessing of the CFD results obtained with simpleFoam is carried out with Tecplot and paraView. The CFD results are shown for flow rates ranging from 5 to 100 mL/min in terms of velocity profiles, streamlines, stagnant zones, swirling strength, and pressure drop. Simulation results obtained with scalarTransportFoam are shown in terms of the dimensionless RTD E(θ). Comparison among different reactor designs is also presented. Velocity fields, Streamlines, and Swirling Strength. The objectives of the CFD simulations carried out for the entire AFR are to (a) determine velocity field along the reactor and the position at which the flow becomes fully developed, (b) determine total pressure drop in the AFR and compare with experimental data, and (c) compute RTD and compare with experimental data. The velocity profiles within each heart-shaped cell are shown in Figure 3 for different flow rates. The maximum velocity is achieved at the inlet of each heart-shaped cell. Reynolds numbers calculated at the inlet of the heart for 5 and 100 mL/ min are 79 and 1580. The streamlines show that the fluid encounters the first obstacle and splits into two streams. Each stream travels at approximately half speed next to the walls, leading to a low7545


Article



velocity region between the two obstacles. Recirculation zones are formed between these two obstacles, increasing the amplitude with increasing flow rates. The velocity profiles in the middle plane of the z coordinate (at z = 0.55 mm) along the cross section at one side of the center of the heart cells is shown in Figure 4. The same plot for

Figure 4. Velocity profile along a cross section on side of heart cell in first row of AFR (black line crossing as shown at top). Water flow rate (mL/min): solid blue trace, 5; solid orange trace, 10; solid red trace, 20; solid black trace, 40; solid green trace, 60; solid violet trace, 80; and dashed orange trace, 100.

the center line is shown in Figure SI-7 of the Supporting Information. In both cases, the velocity profiles are periodic along the series of hearts. The maximum velocities for both locations of the cross sections (side and center) fall under the same parallel lines, only with slight differences encountered in the first heart. This shows that nearly a fully developed flow is already reached within the first heart of the reactor and all heart cells behave in the same way. The swirling strength has been computed for different flow rates for the center plane, and the results are shown in Figure 5.

Figure 5. Swirling strength for different water flow rates (mL/min). Computed from simpleFoam velocity profiles in 3D mesh. Postprocessed in Tecplot. Swirling strength identifies stagnant spots and recirculation zones not detected by other methods.

The amount of swirling strength increases with flow rate and creates stagnation points in those regions where the swirling strength is higher. These stagnation points are located right before the flow encounters the first obstacle within the heart and also at the edges of the long post. Pressure Fields. Figure 6 includes the pressure profile for an AFR segment that illustrates the pressure change within the

Figure 3. Velocity profiles and streamlines for 3D meshes in advancedflow reactor. Simulations using simpleFoam, postprocessing in Tecplot. Velocity profiles are nondimensionalized by the maximum velocity (red), where Vmin = 0 (blue) for all cases and Vmax (red) varies with the flow rate, accordingly (m/s): 0.07 (5 mL/min); 0.15 (10 mL/ min); 0.33 (20 mL/min); 0.60 (40 mL/min); 0.88 (60 mL/min); 1.18 (80 mL/min); 1.44 (100 mL/min). 7546


Article



Figure 6. Pressure drop profile for AFR segment.

heart-shaped cells. A confined annular configuration yields a coflow, followed by a converging/diverging segments that generates a low-pressure zone, followed by a pressure recovery. This helps break the flow in multiphase systems and increase the specific interfacial area. Table 3 includes results for pressure drop in each heart-cell and the entire AFR for different flow rates.The CFD results obtained with OpenFOAM are validated by comparison with experimental results between 20 and 80 mL/min. The CFD results are in very good agreement with the experiments, with relative errors with respect to the experiments below 4% in absolute value. It is known that in classical theory a liquid flowing in laminar flow regime through a pipe has a pressure drop that is linearly proportional to the flow rate or Reynolds number as shown by the Hagen−Poiseuille equation (eq 7).21 ΔP =

8μLQ πr 4

Figure 7. Pressure drop in AFR for incompressible water at different Reynolds numbers using simpleFoam. The nonlinear behavior maybe be due to presence of obstacles in the reactor.

occurs can range from 300 to 2000.23 In addition, several authors have reported in the literature deviations in the predictions of pressure drop in microchannels from the classical theory. This is probably due to the effects of surface roughness at the walls that become very significant at such small scales. Several authors obtained experimental values of the Poiseuille number f Re larger than the conventional theory, such as Harley et al.,24 Rahman and Gui,25 Wilding et al.,26 Peng et al.,27,28 Jiang et al.,29 Mala and Li (0−40%),30 Papautsky et al. (10− 20%),31 Qu et al.,32 Ding et al.,33 Pfund et al.,34 Ren et al.,35 Kandlikar et al.,36 Li et al.,37 and Urbanek et al. (5−30%).38 In other cases, this value was lower than the expected f Re using the classical theory. This was found by Pfalher et al. (0−30%),39 Yu et al. (19%),40 Jiang et al. (50−100%),29 and Judy et al.41 In other cases, no significant differences were found with respect to the classical theory, for example, Harms et al.,42 Pfund et al.,34 and Webb and Zhang.43 Different critical Reynolds numbers have been reported in the literature by several authors to determine the transition regime for different materials, geometries, and channel sizes (Table SI-1 in Supporting Information). Wu and Little reported values for Rec between 350 and 900.44 Other values include those reported by Peng et al. (Rec = 200−700),27,28 Mala and Li (Rec =300−900).45 Higher ranges were observed by Choi et al. (Rec = 2300),46 Yu et al. (Rec = 2000),40 Pfund et al. (Rec = 1700−2200),34 and Wu and Cheng (Rec = 1500− 2000).47 In 2004, Morini23 experimentally investigated the pressure drop through microchannels and observed deviations

(7)

For a fully turbulent fluid, the pressure drop can be estimated using the Darcy−Weisbach22 equation (eq 8), in which the friction factor is calculated using eq 9. In this flow regime, the pressure drop is proportional to the square of the flow rate or the Reynolds number, instead of being linearly proportional as in the laminar flow regime.

ΔP = f

L V2 ρ D 2

(8)

⎛ ε 1 2.51 ⎞ ⎟⎟ = −2 log⎜⎜ r + f Re f ⎠ ⎝ 3.7D h

(9)

The maximum Reynolds numbers at the inlet of the heartcells in the AFR range from 79 to 1580; thus, a laminar flow behavior is expected. However, the nonlinear dependence of pressure drop with the Reynolds number, shown in Figure 7, with a power law ΔP ≈ Re1.44 suggests higher friction at small scales linked to surface roughness or/and the presence of obstacles along the flow path. It is reported in the literature that the critical Reynolds numbers (Rec) at which the transition

Table 3. Comparison of Experimental and CFD Results for Pressure Drop in AFRa

a

QL (mL/min)

Reynolds (−)

Vinlet heart (m/s)

Vmax (m/s)

ΔPAFR (kPa)

ΔPtotal (kPa)

ΔPexperimental (kPa)

relative error (%)

5 10 20 40 60 80 100

79 158 317 634 951 1268 1584

0.076 0.152 0.303 0.606 0.909 1.212 1.515

0.108 0.197 0.360 0.677 0.998 1.362 1.710

0.7 1.6 4.1 11.4 21.7 35.0 51.1

0.81 1.9 5.0 14.0 26.8 43.2 63.1

N/A N/A 4.82 13.8 26.2 44.8 N/A

N/A N/A 3.7 1.4 2.3 −3.6 N/A

Total pressure drop includes additional tubing at the inlet and outlet. 7547


Article


arrives earlier to the outlet of the reactor because of the bypass effect but later than for 60 and 80 mL/min because the extent of the recirculation zones is smaller (and thus the effective volume larger). The RTD curve is not symmetric, and this is due to the enclosed recirculation zones. Not only is the tracer transported from the main stream into the recirculation zones, but also it takes more time to get out because the resistance to the transport is higher as a result of the enclosed double recirculation zones. In practice, it is desirable to have a reactor with well-defined or easily predictable RTD curves. This is not the case for the AFR, in which there are two different behaviors with symmetric RTDs at the lowest (10 and 20 mL/min) and highest flow rates (60 and 80 mL/min) and a transition region with nonsymmetric RTDs for intermediate flow rates (40 mL/min). A detailed study of the influence of reactor design on the velocity profiles and RTD is described later in this work (Effect of Reactor Design). The next step is to compare the RTD obtained from CFD simulations with experimental results. The CFD overall RTD is obtained by convolution of the RTD of the AFR and the RTD of the tubing at the inlet (170 mm length, 1/16 in. O.D.) and the outlet (130 mm length, 1/8 in. O.D., and 160 mm length, 1/16 in. O.D.) of the reactor with the measuring system similar to the experimental setup. The dimensionless RTD curves for different flow rates are shown in Figure 9. As observed, the convoluted RTD obtained with CFD is in agreement with the experimental RTD.


from conventional behavior for hydraulic diameters less than 1 mm. The transition from laminar to turbulent occurs at Reynolds numbers lower than those of conventional channels and is explained in the case of the Obot−Jones model and the increase in the surface roughness. Morini also concluded that the fluid flow in microchannels still remains a topic of future research for a full understanding of this phenomenon. Residence Time Distribution. One of the essential steps in the scale-up of reactors is the determination of the RTD. The dimensionless RTD curves for the AFR are shown in Figure 8.

Figure 8. Dimensionless RTD E(θ) for different flow rates. Water flow rate (mL/min): solid orange trace, 10; solid red trace, 20; solid lightblue trace, 40; solid dark-blue trace, 60; and solid green trace, 80. A plug flow behavior with certain degree of dispersion is observed, with a change of regime above 40 mL/min where recirculation zones are present.

We observe a different regime at the lowest flow rates (10 and 20 mL/min) than at the largest flow rates (60 and 80 mL/min), with a transition regime at 40 mL/min. At the lowest flow rates, the RTD is sharper than at the highest flow rates, whereas the intermediate flow rate causes a nonsymmetric RTD curve. This reflects the complexity of the geometry and thus of the flow in the AFR. In a single tube, one expects that larger flow rates increase the Peclet number and reduce the axial dispersion, with the result of a narrower E(θ) curve.48 This is not what is observed in the AFR according to the CFD simulations. An explanation is found on the basis of the streamlines obtained and how they change for the different flow rates. Looking back at Figure 3, flow rates from 5 to 20 mL/min do not show any recirculation zone, which yields a narrow RTD. Increasing the flow rate to 60 and 80 mL/min broadens the RTD. This is due to two main contributions: (a) The recirculation zones that occupy the region between the two obstacles cause the tracer to bypass these regions and arrive at the outlet of the reactor earlier than if no recirculation zones are present. (b) The mass transfer between the recirculation zones and the main flow stream causes the tracer to be transported first from the main stream to the recirculation zones and then to leave the recirculation zones and follow the main stream toward the outlet of the reactor. This causes a delayed arrival of the tracer and, together with the bypass effect, broadens the RTD. A different behavior is observed for 40 mL/min. Here, recirculation zones are also present, but there are two smaller recirculation zones enclosed in a larger one. First, the tracer

Figure 9. Comparison of CFD and experimental RTD, E(θ). Simulation results obtained for 3D mesh using scalarTransportFoam and differentiation of signal obtained at the outlet using a step function for the inlet concentration. Flow rates (mL/min): (a) 60 and (b) 80. Legend: dashed blue trace, experimental 1; dashed green trace, experimental 2; and dashed red trace, CFD. 7548


Article



Figure 10. Geometry for three AFR designs studied.

Figure 11. Comparison of velocity fields and streamlines among three different AFR designs (dot, no dot, low-flow reactor-like) for incompressible flow using simpleFoam and a 3D mesh. Simulations using simpleFoam, postprocessing in Tecplot. Nondimensionalized velocity fields with maximum velocity (Vmax), where Vmin = 0 (blue) for all cases and Vmax (red) varies with the flow rate, accordingly (m/s): 0.15 (10 mL/min); 0.33 (20 mL/min); 0.60 (40 mL/min); 0.88 (60 mL/min); 1.18 (80 mL/min).

Effect of Reactor Design. There are different versions of the AFR that Corning, Inc., has created with slightly different designs as shown in Figure 10. Modifications in the reactor design can have a significant impact on the reactor performance; thus, determination of the hydrodynamic characteristics and comparison between designs can provide insightful information on this matter. Here we analyze velocity fields and pressure drops for the three different designs and

determine through CFD simulations the RTD for singlephase flow. Simulations are carried out in 3D meshes for the entire AFR reactors. The velocity profiles and streamlines for the three different reactor designs at 10, 20, 40, 60, and 80 mL/min are included in Figure 11. A common characteristic at low flow rates is the absence of recirculation zones for the three designs. The largest difference appears at 40 mL/min, when recirculation zones start 7549


Article


extension of recirculation zones increase. The overall effect is to broaden the RTD symmetrically as Re increases for design 2. In design 3, increasing the flow rate from 60 to 80 mL/min does not increase the extent of recirculation zones; it has the contrary effect instead. This causes a narrower RTD. Because of the aerodynamic design of the middle section, it does not enhance wake formation and hence allows the fluid to move directly to the next segment. The pressure drop across the entire reactor for the different AFR designs is shown versus the Reynolds number at the inlet of the heart cell in Figure 13. It is observed that the three


to be present in designs 1 and 2, whereas they are almost absent in design 3. Design 1 has two enclosures of recirculation zones, whereas design 2 only has two single recirculation zones at both sides of the heart cell. Thus, one would expect the RTD to be more symmetric for design 2 as compared to design 1. Regarding design 3, the RTD is expected to be narrower because no significant recirculation zones are present. The RTDs for the three designs are shown in Figure 12.

Figure 13. Effect of reactor design on pressure drop. Results obtained with simpleFoam using a 3D mesh. Legend: open black circle, design 1, 8 mL; open blue square, design 2, 8 mL volume; open red triangle, design 3, 8 mL volume; and filled green diamond, LFR, 0.45 mL volume. Larger pressure drops are observed for the same Reynolds number in the LFR than those in the AFR.

designs have similar pressure drops, especially at the lowest Re, but differences start to increase at larger Re. Design 3 has an average channel width smaller than that of designs 1 and 2, and this may be the reason why pressure drops for this reactor design are higher. If we compare the AFR with the LFR of 0.45 mL volume, then the total pressure drops are even larger. This is due to the smaller reactor height and width. The height dimension for the AFR of 8.7 mL volume is 1.1 mm, whereas the LFR has a height of 0.35 mm. The maximum recommended operating flow rate for the LFR is 10 mL/min, for which the pressure drop almost equals the pressure drop in the AFR at 100 mL/min. An improvement in the design has been made in the LFR, removing the space between the two posts in the original AFR design and using a single post. As shown in Figure 14, this design removes the formation of large recirculation zones for the recommended flow rates, between 0.5 and 10 mL/min. It is expected that the RTD for this reactor design will be sharper and closer to that of plug flow behavior with less axial dispersion, especially at the largest flow rates. This is shown in Figure 15, where it is seen also that the RTDs are symmetric, as expected from the streamlines in the LFR. A small tail in the RTD is observed at the largest flow rate. This is due to the appearance of small recirculation zones, which makes the tracer bypass these zones while mass transfer occurs from the main stream. Later on, mass transfer transports the tracer from the recirculation zones to the main stream, taking more time to

Figure 12. Dimensionless RTDs E(θ) for different AFR designs: (a) design 2 and (b) design 3. Simulation results obtained for 3D mesh using scalarTransportFoam and differentiation of signal obtained at the outlet using a step function for the inlet concentration. Water flow rate (mL/min): solid red trace, 10; dashed orange trace, 20; solid blue trace, 40; dashed green trace, 60; and solid violet trace, 80.

The extent of recirculation zones increases significantly from 40 to 60 mL/min for all designs. A further increase in the flow rate increases the recirculation strength in design 2. However, the second obstacle in design 1 prevents the wake from extending further downstream and makes the streamlines to be very similar to all flow rates above 60 mL/min. This causes the RTD to be the same for design 1 at flow rates above this flow rate. In design 2, however, the extent of recirculation zone still increases with flow rate. This causes the tracer first to bypass these zones of increasing volume (and thus reducing the effective reactor volume), causing the tracer to arrive earlier at the outlet of the reactor as the flow rate increases. At the same time, the tracer takes longer time to arrive at the outlet as the 7550


Article



Figure 14. Velocity profiles and streamlines in LFR. Results obtained with simpleFoam for 3D mesh. Water flow rates (mL/min): a) 0.5; b) 1; c) 5; d) 8; e) 10; Velocity profiles are nondimensionalized with respect to maximum velocity (Vmax), where Vmin = 0 m/s (blue) for all cases and Vmax (red) varies with the flow rate, accordingly (m/s): 0.096 (0.5 mL/min); 0.18 (1 mL/min); 0.80 (5 mL/min); 1.20 (8 mL/min); 1.47 (10 mL/min).

certain degree of dispersion, bypass, and mass transfer effects, the degree of each phenomenon is different depending on the flow rates and the design. To quantify the importance of the differences between designs in the actual performance of the reactor when reaction is occurring, the average reagent concentration is calculated from the RTD and compared with two ideal cases: (a) plug flow (PFR) and (b) single stirred tank (CSTR). The results are shown in Figure 16 for different Damköhler numbers and designs. The model reaction to study is an irreversible reaction A + B → C, with first-order (n = 1) kinetics r = kCA, and second-order (n = 2) kinetics, r = kCA2. The Damkhöler number is represented as eq 10. Da = τkC A 0 n − 1 Figure 15. Residence time distribution in LFR. Simulation results obtained for 3D mesh using scalarTransportFoam and differentiation of signal obtained at the outlet using a step function for the inlet concentration.

(10)

From the results, it can be seen that the conversion for all reactor designs is not significantly different, and all of them are very close to the plug flow conversions. A single CSTR gives the minimum conversion that can be achieved in the system and is below the actual conversion achieved in the AFR (8 mL volume) and the LFR (0.45 mL volume). A slightly better performance is achieved in the LFR, which has sharper RTD curves and this reflects a closer conversion to PFR. However, the difference is not significant, especially at the smallest Da values (high flow rates, small residence time, and narrower RTD). This shows once more that the different generations of AFR have been designed to provide a sharp RTD regardless of the flow rates and slightly different designs. The design of individual segments enhances mixing locally, which is similar to a CSTR, and a sequence of such segments resembles a series of CSTRs finally leading to a performance similar to a PFR.

arrive at the outlet of the reactor. Because of the very small extent of recirculation zones, this effect is not as important as in the AFR design. One important thing to consider here is the ease of scalability of the different Corning reactors. The best case scenario would be to have the same RTD for the same Reynolds numbers (or other variable of reference for scale up) on the different Corning scales. In this way, the overall reaction kinetics would be affected in the same way at the different scales, and the scalability between different sizes would be straightforward. However, this is not the case for the designs shown in this work. Although the reactors have plug-flow behavior with a

Figure 16. Conversion of reagent A for an irreversible reaction A + B → C. (a) First-order kinetics: r = kCA. (b) Second-order kinetics: r = kCA2. Legend: dashed red line, plug flow; dashed black line, CSTR; filled violet diamonds, LFR 0.45 mL volume; filled green triangles, design 1; open red circles, design 2; and open black squares, design 3. 7551


Article



■

g⃗ L k n P ΔP Psim,i Pexp,i Q QL t v⃗ v′ v̅ V Vmax Vmin Vinlet heart r xAF z

CONCLUSIONS Computational fluid dynamic (CFD) simulations are carried out in the advanced-flow reactor (AFR) for single-phase incompressible flow using the open source software OpenFOAM. Results showed a very good agreement with experimental results of pressure drop and RTD analyses. The velocity profiles reveal the presence of stagnant zones between the two obstacles within the heart-shaped Gen 1 reactor and recirculation zones starting at 40 mL/min. The RTDs for the AFR are plug flow with certain degree of back mixing, with two different flow regimes. The extent of recirculation zones, which increases with flow rates, broadens the RTD above 40 mL/min. Introduction of slight differences in the design in the zone between the two obstacles slightly changes the flow patterns and recirculation zones, therefore having an effect on the RTD. However, predictions for a first-order irreversible reaction in all designs yielded conversions that are not significantly different. The low-flow reactor (LFR) of 0.45 mL volume is also simulated here. Results showed a higher pressure drop for the same Reynolds numbers than those of the Gen 1 AFR, which limits the operating flow rate to a maximum of 10 mL/min. In contrast, sharper RTD curves are obtained in the LFR with respect to the AFR, yielding conversions slightly larger for large Damkhöler numbers. Although the AFR is able to achieve 95− 99% times the PFR conversion, conversions above 99% of the PFR conversion is achieved in the LFR. However, the LFR has the limitation of low throughput, whereas the AFR can handle flow rates on the order of 100 mL/min.

■

DimensionlessNumbers

Da Damkhöler number (−) Re Reynolds number (−) Rec critical Reynolds number (−) Greek Letters

ρ εr θ τ μ

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b00232.

■

gravity vector (m/s2) characteristic length (m) kinetic constant (Ln−1 mol1−n s−1) reaction order (−) pressure (Pa) pressure drop (Pa) result from simulation result from experiment flow rate (mL/min) liquid flow rate (mL/min) time (s) velocity vector (m/s) velocity fluctuation (m/s) averaged velocity (m/s) velocity magnitude (m/s) maximum velocity (m) minimum velocity (m) velocity at heart inlet (m) channel radius (m) reagent conversion (%) coordinate in z axis

■

AUTHOR INFORMATION

fluid density (kg/m3) material rugosity (m) dimensionless time (−) residence time (−) fluid viscosity (Pa s)

REFERENCES

(1) Jensen, K. F. Microreaction engineering − is small better? Chem. Eng. Sci. 2001, 56, 293−303. (2) Roberge, D. M.; Ducry, L.; Bieler, N.; Cretton, P.; Zimmermann, B. Microreactor technology: a revolution for the fine chemical and pharmaceutical industries? Chem. Eng. Technol. 2005, 28 (3), 318− 323. (3) Roberge, D. M.; Gottsponer, M.; Eyholzer, M.; Kockmann, N. Industrial design, scale-up, and use of microreactors. Chim. Oggi 2009, 27, 8−11. (4) Kockmann, N.; Gottsponer, M.; Roberge, D. M. Scale-up concept of single-channel microreactors from process development to industrial production. Chem. Eng. J. 2011, 167 (2−3), 718−726. (5) Corning Inc. www.corning.com/WorkArea (accessed on April 3, 2014). (6) Weller, H. OpenFOAM; CFD Direct: Reading, United Kingdom, 2015; http://www.openfoam.org (accessed on April 3, 2014). (7) McMullen, J. P.; Stone, T. M.; Buchwald, S. L.; Jensen, K. F. An Integrated Microreactor System for Self-Optimization of a Heck Reaction: From Micro- to Mesoscale Flow Systems. Angew. Chem., Int. Ed. 2010, 49 (39), 7076−7080. (8) Monbaliu, J-C. M.; Winter, M.; Chevalier, B.; Schmidt, F.; Jiang, Y.; Hoogendoorn, R.; Kousemaker, M. A.; Stevens, C. V. Bioresour. Technol. 2011, 102 (19), 9304−9307. (9) Rossi, E.; et al. Scalable in Situ Diazomethane Generation in Continuous-Flow Reactors. Org. Process Res. Dev. 2012, 16 (5), 1146− 1149. (10) Buisson, B.; Donegan, S.; Wray, D.; Parracho, Ana.; Gamble, J.; Caze, P. Slurry hydrogenation in a continuous flow reactor for pharmaceutical application. Chim. Oggi 2009, 27, 12−16. (11) Woehl, P. Advanced-Flow Reactors: High throughput continuous flow reactor technology; Corning, Inc.: Corning, NY, 2010; http://www.

Corresponding Author

*E-mail: [email protected]. Present Address

M.J.N.-R.: The Dow Chemical Company, 2301 North Brazosport Blvd., Building B-1603, Freeport, TX 77541, USA. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge Novartis-MIT Center for Continuous Manufacturing for financial support and Corning, Inc., for the Advanced-Flow Reactor. M.J.N.-R. thanks the Fundación Cajamadrid (Spain) for the fellowship to pursue postgraduate studies, and A.A.K. thanks the Indo-US Science and Technology Fellowship (IUSSTF) for funding during stay at the Massachusetts Institute of Technology.

■

NOTATIONS CA concentration of reagent A (mol/L) CT tracer concentration (mol/L) Dh hydraulic diameter (m) DT tracer diffusion coefficient (m2/s) D determinant to determine vortices from velocity field E(t) residence time distribution from pulse experiment (s−1) E(θ) dimensionless residence time distribution of E(t) (−) f Darcy’s friction factor (−) 7552


Article



International Mechanical Engineering Congress and Exposition: Dynamic Systems and Control; ASME: New York, 1998; Vol. 66, pp 193−198. (35) Ren, L.; Qu, W.; Li, D. Interfacial electrokinetic effects on liquid flow in microchannels. Int. J. Heat Mass Transfer 2001, 44, 3125−3134. (36) Kandlikar, S. G.; Joshi, S.; Tian, S. Effect of channel roughness on heat transfer and fluid flow characteristics at low Reynolds numbers in small diameter tubes. In Proceedings of the 35th National Heat Transfer Conference, Anaheim, California, 2001; ASME: New York, 2001; paper 12134. (37) Li, Z. X.; Du, D. X.; Guo, Z. Y. Experimental study on flow characteristics of liquid in circular microtubes. In Proceedings of the International Conference on Heat Transfer and Transport Phenomena in Microscale; Celata, G.P., Ed.; Begell House, New York, 2000; pp 162− 168. (38) Urbanek, W.; Zemel, J. N.; Bau, H. An investigation of the temperature dependence of Poiseuille numbers in microchannel flow. J. Micromech. Microeng. 1993, 3, 206−208. (39) Pfahler, J.; Harley, J.; Bau, H.; Zemel, J. N. Gas and liquid flow in small channels. In Micromechanical Sensors Actuators and Systems, Proceedings of the Winter Annual Meeting of the American Society of Mechanical Engineers, Atlanta, Georgia, December 1-6, 1991; Cho, D., Ed.; ASME: New York, 1991; Vol. 32, pp 49−58. (40) Yu, D.; Warrington, R.; Barron, R.; Ameel, T. Experimental and theoretical investigation of fluid flow and heat transfer in microtubes. In Proceedings of the 1995 ASME/JSME Thermal Engineering Joint Conference, Maui, Hawaii, 1995; ASME: New York, 1994; Vol. 1, pp 523−530. (41) Judy, J.; Maynes, D.; Webb, B. W. Characterization of frictional pressure drop for liquid flows through microchannels. Int. J. Heat Mass Transfer 2002, 45 (17), 3477−3489. (42) Harms, T. M.; Kazmierczak, M.; Gerner, F. M.; Holke, A.; Henderson, H. T.; Pilchowski, J.; Baker, K. Experimental investigation of heat transfer and pressure drop through deep microchannels in a (110) silicon substrate. In Proceedings of the ASME Heat Transfer Division-1997; American Society of Mechanical Engineers: New York, 1997; Vol. 1, pp 347−357. (43) Webb, R. L.; Zhang, M. Heat transfer and friction in small diameter channels. Microscale Thermophys. Eng. 1998, 2, 189−202. (44) Peiyi, W.; Little, W. A. Measurement of friction factors for the flow of gases in very fine channels used for microminiaturize Joule− Thomson refrigerators. Cryogenics 1983, 23 (8), 273−277. (45) Mohiuddin Mala, Gh.; Li, D. Flow characteristics of water in microtubes. Int. J. Heat Fluid Flow 1999, 20 (2), 142−148. (46) Choi, S. B.; Barron, R. F.; Warrington, R. O. Fluid flow and heat transfer in microtubes. In Micromechanical Sensors Actuators and Systems, Proceedings of the Winter Annual Meeting of the American Society of Mechanical Engineers, Atlanta, GA, December 1−6, 1991; Vol. 32, pp 123−134. (47) Wu, H. Y.; Cheng, P. Friction factors in smooth trapezoidal silicon microchannels with different aspect ratios. Int. J. Heat Mass Transfer 2003, 46, 2519−2525. (48) Levenspiel, O. Chemical Reaction Engineering, 3rd ed.; John Wiley & Sons: New York, 1999.

ifpac.com/cortona/Cortona10-Presentations/Woehl.pdf (accessed on June 1, 2015). (12) Meinhart, C. D.; Wereley, S. T.; Santiago, J. G. PIV measurements of a microchannel flow. Exp. Fluids 1999, 27, 414. (13) Pope, S. B. Turbulent flows; Cambridge University Press: Cambridge, U.K., 2000. (14) Bardina, J. E.; Huang, P. G.; Coakley, T. J. Turbulence Modeling Validation, Testing, and Development; NASA Technical Memorandum 110446; National Aeronautics and Space Administration: Moffett Field, CA,1997; pp 94035−1000. (15) Lam, C. K. G.; Bremhorst, K. A modified form of the k-epsilon model for predicting wall turbulence ASME, Transactions. J. Fluids Eng. 1981, 103 (3), 456−460. (16) Nagano, Y.; Hishida, M. Improved form of the k-ε model for wall turbulent shear flows. J. Fluids Eng. 1987, 109 (2), 156−160. (17) Myong, H. K.; Kasagi, N. A new approach to the improvement of k-epsilon turbulence model for wall-bounded shear flows. JSME Int. J., Ser. II 1990, 33, 63−72. (18) Bakker, A. The Colorful Fluid Mixing Gallery; 2012; http://www. bakker.org/cfm (accessed on March 15, 2014). (19) Jasak, H. Error analysis and estimation for the Finite Volume method with applications to fluid flows. Ph.D. Thesis. Imperial College, University of London,1996. (20) Adrian, R. J.; Christensen, K. T.; Liu, Z. C. Analysis and interpretation of instantaneous turbulent velocity fields. Exp. Fluids 2000, 29, 275−290. (21) Deen, W. M. Analysis of transport phenomena; Oxford University Press: New York, 1998. (22) Brown, G. The History of the Darcy-Weisbach Equation for Pipe Flow Resistance. Environ. Water Resour. Hist. 2003, 34−43. (23) Morini, G. L. Laminar-to-turbulent flow transition in microchannels. Microscale Thermophys. Eng. 2004, 8 (1), 15−30. (24) Harley, J.; Huang, Y.; Bau, H. H.; Zemel, J. N. Gas flow in microchannels. J. Fluid Mech. 1995, 284, 257−274. (25) Rahman, M. M.; Gui, F. J. Experimental measurements of fluid flow and heat transfer in microchannel cooling passages in a chip substrate. In Advances in Electronic Packaging 1993, Proceedings of the 1993 ASME International Electronics Packaging Conference, Binghamton, New York, September 29−October 2, 1993; ASME: New York, 1993; Vol. 2, pp 685−692. (26) Wilding, P.; Pfalher, J.; Zemel, J. N.; Bau, H. H.; Kricka, L. J. Manipulation and flow of biological fluids in straight channels micromachined in silicon. Clin. Chem. 1994, 40, 43−47. (27) Peng, X. F.; Peterson, G. P. Convective heat transfer and flow friction for water flow in microchannel structures. Int. J. Heat Mass Transfer 1996, 39, 2599−2608. (28) Peng, X. F.; Peterson, G. P.; Wang, B. X. Frictional flow characteristics of water flowing through rectangular microchannels. Exp. Heat Transfer 1994, 7, 249−265. (29) Jiang, X. N.; Zhou, Z. Y.; Huang, X. Y.; Liu, C. Y. Laminar flow through microchannels used for microscale cooling systems. In Proceedings of the 1997 1st Electronic Packaging Technology Conference, Singapore, October 10, 1997; Tay, A. A. O., Lim, T. B., Eds.; Institute of Electrical and Electronics Engineers, Inc.: Singapore, 1997; pp 119− 122 (30) Mohiuddin Mala, G. M.; Li, D.; Dale, J. D. Heat transfer and fluid flow in microchannels. Int. J. Heat Mass Transfer 1997, 40, 3079− 3088. (31) Papautsky, I.; Gale, B. K.; Mohanty, S.; Ameel, T. A.; Frazier, A. B. Effects of rectangular microchannel aspect ratio on laminar friction constant. Proc. SPIE 1999, 3877, 147−158. (32) Qu, W.; Mala, M.; Li, D. Pressure-driven water flows in trapezoidal silicon microchannels. Int. J. Heat Mass Transfer 2000, 43, 3925. (33) Ding, L. S.; Sun, H.; Sheng, X. L.; Lee, B. D. Measurement of friction factors for R134a and R12 through microchannels. Proc. Symp. Energy Eng. 21st Century 2000, 2, 650−657. (34) Pfund, D.; Shekarriz, A.; Popescu, A.; Welty, J. R. Pressure drop measurements in a microchannel. In Proceedings of the 1998 ASME 7553


OpenFOAM Computational Fluid Dynamic Simulations of Single

Recommend Documents