- Email: [email protected]
Enhanced heat transfer by exothermic reactions in laminar flow capillary reactors
Chemical Engineering Science 141 (2016) 356–362
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces
Enhanced heat transfer by exothermic reactions in laminar flow capillary reactors Sebastian Schwolow a, Jing Ying Ko a, Norbert Kockmann b, Thorsten Röder a,n a b
Mannheim University of Applied Sciences, Institute of Chemical Process Engineering, Paul-Wittsack-Straße 10, 68163 Mannheim, Germany TU Dortmund, Biochemical and Chemical Engineering, Equipment Design, Emil-Figge-Straße 68, 44227 Dortmund, Germany
H I G H L I G H T S
G R A P H I C A L
A B S T R A C T
For a capillary reactor, concentration and temperature profiles are discussed. Exothermic reactions can significantly enhance heat transfer. For reactor design, commonly used plug flow models can lead to miscalculations. Simulating both axial and radial conversion profiles can be a key step for reactor scale-up.
art ic l e i nf o
a b s t r a c t
Article history: Received 15 September 2015 Received in revised form 20 November 2015 Accepted 21 November 2015 Available online 29 November 2015
Conversion and temperature distributions in a tubular reactor with an inner diameter of 1 mm were numerically calculated for a second order exothermic reaction in laminar flow of homogeneous liquid. Based on the resulting radial temperature profiles, the local Nusselt numbers and bulk mean temperatures were determined along the reactor tube. Strong effects of the homogeneous reaction on the heat transfer can be observed in the entrance region of the reactor, where a hot spot emerges. Due to the large radial temperature gradients in vicinity of the reactor wall, heat transfer coefficients are significantly higher compared to a non-reactive system. The consequences of this effect on the design and control of exothermic reactions in reactor/heat exchangers are demonstrated by comparison with a simple one-dimensional plug flow model. In the simplified model, neglecting thermal influence of the exothermic reaction results in a significant underestimation of the required reactor length for defined conversion. Accordingly, numerical simulation of both axial and radial transport in the hot spot region can be essential to precisely predict the bulk temperatures and conversion rates in the reaction mixture, even with the small length scales of milli- and micro-structured reactors. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Capillary reactor Heat transfer Simulation Nusselt number Temperature distribution Entrance flow
1. Introduction Milli- and micro-structured flow reactors have a high potential for process intensification. The finely adjustable reaction times and excellent heat removal in these reactors make it possible to run exothermic reactions in safe and defined conditions (Kockmann and Roberge, 2009; Taghavi-Moghadam et al., 2001; Zhang et al., 2004). To calculate the conversion profiles in microreactors, plug
n
Corresponding author. Tel.: þ 49 621 292 6800. E-mail address: [email protected] (T. Röder).
http://dx.doi.org/10.1016/j.ces.2015.11.022 0009-2509/& 2015 Elsevier Ltd. All rights reserved.
flow and often isothermal assumptions are widely used. Both simplifications are based on the intensified transport processes with the small channel diameters and large reactor length-todiameter ratios: 1) With very short diffusion paths, radial mixing by molecular diffusion can be fast enough to approach plug flow behavior, even in the laminar flow regime. Further enhancement of radial mixing can be obtained by inducing a secondary flow with the specific channel design (Klutz et al., 2015). 2) Heat transfer within the reactor is very efficient as a result of the high surface-to-volume ratio and small length scale for thermal conduction. Thus, the rate of heat removal is assumed to be fast compared to the rate of heat generation and, as a result, the
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
reaction mixture may be considered to be nearly isothermal. However, at a closer look, these common idealizations can significantly deviate from the actual reactor behavior: If the hydrodynamic residence time is much shorter than the characteristic time for molecular diffusion, radial concentration gradients exist, and either the convection model for laminar flow or the dispersion model are more accurate representations of the residence time distribution (Levenspiel, 1999). In the case of very fast and strongly exothermic reactions, even the superior heat transfer characteristics of microreactors might not be sufficient to avoid hot spots in the entrance region of the reactor (Barthe et al., 2008). Local temperature peaks accelerate the reaction and significantly influence the yield in the reaction system as a result. While numerous publications describe experimental characterization and modeling of the residence time distribution in microreactors (Boskovic and Loebecke, 2008; Lohse et al., 2008; Trachsel et al., 2005; Wörner et al., 2007), less consideration is given to the temperature distribution in micro- and millireactors. Nevertheless, the formation of hot spots is a crucial concern for reactor design and control. Because temperature measurements in small spaces of microreactors are often difficult to perform in practice, precise calculation of local temperatures has additional significance. For a rough estimation, the internal heat transfer coefficient, α, can be determined from correlations for the dimensionless Nusselt number, Nu ¼ αλd, in the channel flow (Martin, 2002). For a straight tubular reactor (e.g., a stainless steel capillary in an oil bath), the internal heat transfer can be calculated by assuming a constant wall temperature and fully developed flow. With these assumptions, a constant Nusselt number, Nu ¼ 3:6568, can be derived for laminar flow (Baehr and Stephan, 2006). However, radial temperature profiles can be strongly influenced by the reaction itself. This effect is completely neglected by using common correlations or assuming a constant Nusselt number throughout the entire channel. An intensive and comprehensive insight into the phenomena resulting from an exo- or endothermic chemical reaction with concomitant heat transfer has been published by Churchill et al. (Churchill, 2005; Churchill et al., 2008; Yu and Churchill, 2012). Churchill describes strong effects and complex behavior that significantly enhance heat transfer in the presence of a homogeneous exothermic reaction. Because additional parameters are introduced by one or several chemical reactions, he concluded that a general correlating equation cannot be derived practically. Instead, he recommends that numerical calculation of both the radial and axial transport be performed on a case-by-case basis. Radial temperature profiles in tubular reactors have been investigated for some typical gaseous reaction systems, for example, by Rothenberg and Smith (1966), or by Merrill and Hamrin (1970). The objective of the present study is the detailed investigation and discussion of heat transfer in a liquid phase reaction, based on exemplary process specifications that closely approximate the typical operating conditions in milli- and microstructured flow reactors. A fast, exothermic, second-order reaction in a circular duct with 1 mm diameter was simulated using two different models: 1) A simple one-dimensional plug flow model assuming a constant Nusselt number (heat transfer in thermally developed flow) and 2) a two-dimensional model of the rotationally symmetric geometry that accounts for the axial and radial heat and mass transport. Visualization and discussion of radial temperature and concentration profiles in the investigated capillary reactor provide the basis for an improved understanding of the effects previously described by Churchill and coworkers. By comparing results of both models, the applicability of the simplified model can be evaluated with regard to practical aspects in the design and control of reactor/heat exchangers. Furthermore, the influence of the exothermic reaction is demonstrated by comparing the reactor/heat exchanger with a corresponding heat exchanger without chemical reaction.
357
2. Basic model assumptions The exothermic reaction is assumed to be second order and irreversible and to take place in a homogeneous liquid phase (Table 1). The temperature dependence of the reaction rate was described by the Arrhenius equation, in which the activation energy was varied to ascertain the thermal effects. The material properties used were based on toluene and assumed to be constant for all temperatures and conversions. For the purpose of this study, a reaction medium with a density of 867 kg m 3, heat capacity of 1707 J kg 1 K 1, thermal conductivity of 0.134 W m 1 K 1, and a mass diffusivity of 10 9 m2 s 1 for the reactant was considered. The considered reactor is a tubular flow reactor (Table 2). Due to its rotational symmetry, the three-dimensional problem can be reduced to a two-dimensional axisymmetric problem. A constant wall temperature, Tw, was defined and for all simulations with an exothermic reaction, the inlet temperature of the reactor was defined as T0 ¼Tw. The partial differential equations resulting from the mass and energy balance in the tube flow are ∂X A D 1 ∂X A ∂2 X A ∂2 X A rA þ 2 þ 2 ¼0 ∂z u r ∂r ucA;0 ∂r ∂z
ð1Þ
∂T λ 1 ∂T ∂2 T ∂2 T r A ΔHR þ 2þ 2 þ ¼0 ∂z uρcp r ∂r ∂r uρcp ∂z
ð2Þ
The following boundary conditions that were used in this study: T ðr; z ¼ 0Þ ¼ T 0
X A ðr; z ¼ 0Þ ¼ 0
T ðr ¼ R; zÞ ¼ T w
dX A ðr ¼ R; zÞ ¼ 0 dr
dT ðr ¼ 0; zÞ ¼ 0 dr
dX A ðr ¼ 0; zÞ ¼ 0 dr
For all z, the velocity profile in laminar Hagen–Poiseuille flow is defined as r 2 uðr Þ ¼ 2um 1 R
ð3Þ
With these definitions, it was assumed that the reaction was initiated in a hydrodynamic fully developed flow with a homogeneous radial reactant distribution. Table 1 Reaction parameters used for the simulations. Rate constant Reference temperature Initial reactant concentration Reaction enthalpy
0.0004 30 4 40
m3 mol 1 s 1 °C mol L 1 kJ mol 1
Table 2 Reactor parameters used for the simulations. Reactor diameter Reactor length Mean flow velocity Wall temperature Reactor pressurea
1 1 0.1 25 5
mm m m s1 °C bar
a The reactor pressure was chosen in order to avoid exceeding the boiling temperature in the reactor tube for the reference material toluene.
358
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
Fig. 1. Temperature and conversion profiles in the entrance region of the reactor, with an exothermic second order reaction (ΔTad ¼ 108 K and EA ¼ 50 kJ/mol).
To determine the longitudinal temperature profile in the tube from the results, fluid bulk mean temperatures were calculated: Z R 1 Tm ¼ 2 2r πuðr ÞT ðr Þdr ð4Þ R πum 0 This temperature, also referred to as “mixing cup temperature”, is the bulk temperature that theoretically would be measured if the channel was sliced at any cross-section and the outflowing fluid was captured in an adiabatic, well-mixed cup. Bulk mean conversions were determined in an analogous manner. The local Nusselt number in the reactor can be obtained, when the heat transfer coefficient is calculated from the wall heat flux density divided by Tm Tw:
growth rate of 1.031. Thus, a very fine mesh resolution could be obtained in the wall region with a minimum element width of 0.7 mm.
3. Results and discussion 3.1. Heat transfer enhancement
ð10Þ
In the absence of thermal energy sources, the temperature distribution in a fully developed Hagen-Poiseuille flow can be calculated based on the energy balance. A derivation of the temperature distribution for constant wall temperature can be found in the literature (Baehr and Stephan, 2006) and results in a rapidly convergent power series. The fully developed Nusselt number corresponding to this temperature distribution was shown to be Nu¼ 3.66. For the case of a thermal entrance in hydrodynamically developed flow, the local Nusselt numbers in the thermal entrance region are higher compared to those in the fully developed flow region, and can be calculated with correlations given in the literature (Gnielinski, 1989; Shah and London, 1978). With an exothermic chemical reaction, a thermal energy source is added to the heat transfer problem. Due to the dependence of the heat release on the local reactant concentration and temperature, radial and axial transport has to be calculated numerically. The results are shown in Fig. 1, where the development of the radial temperature and conversion profiles in the entrance region is pictured. Herein, fluid enters the tube with the wall temperature (T0 ¼Tw), and the internal heat source in the fluid is determined from the second order reaction equations. On a longitudinal length scale of z 50 Ud, a hot spot is formed. The reaction rate is high near the reactor inlet because of the high reactant concentrations. As the fluid moves forward in the reactor channel, the heat generated by the reaction decreases with the increase in reactant conversion. Simultaneously, the heat removal rate increases with the temperature difference Tm Tw. The temperature and conversion profiles shown in Fig. 1 exhibit thermal effects with mutual interactions. The formation of the hot spot enhances the reaction rate, therefore, a fast increase in local conversion can be observed at the hot spot.
The differential equations of the plug flow model were solved using the ode15s-solver in Matlab (Mathworks Inc., Natick, MA, USA). The two-dimensional axisymmetric problem was solved with the finite element software COMSOL Multiphysics (COMSOL, Inc., Palo Alto, CA, USA). Herein, a structured quadrilateral mesh with a 720,000 elements was used. The highest gradients can be expected at the reactor inlet and close to the reactor wall. For this reason, mesh refinement was applied with element size distributions defined by arithmetic sequences in the z- and r-directions, with an average
3.1.1. Influence of the exothermic reaction on radial temperature profiles To discuss the thermal effects, Fig. 2 shows a comparison between cooling a non-reactive hot inlet stream (2 a) and cooling a fluid which is increasing in temperature due to an exothermic reaction (2b). In Fig. 2a, the selected radial temperature profiles show the effect of the thermal entrance without a chemical reaction and T0 Tw ¼ 108 K. Corresponding temperature profiles from the simulation that accounts for the exothermic reaction are
Nu ¼
q_ w 2R αd ¼ λ ðT m T w Þλ
ð5Þ
Herein, qw was determined from the enthalpy balance in a cylindrical tube element 2 dT m q_ w ¼ um ρcp r A ΔH R R dz
ð6Þ
For the simplified, one-dimensional model, the usual assumptions of plug flow were applied, which reduce Eqs. (1) and (2) to dX A rA ¼0 dz cA;0 u
ð7Þ
dT rA ΔH R þ q_ w R2 þ ¼0 dz uρcp
ð8Þ
In Eq. (8), q_ w R2 is the heat transfer rate to the wall over the reactor length with :
qw ¼ α ðT T w Þ
ð9Þ
The heat transfer coefficient α was defined by assuming fully developed, laminar flow heat transfer with Nu ¼ 3.66. For the case of heat transfer without chemical reactions, the local Nusselt numbers in the thermal entrance region can be calculated by the following correlation (Shah and London, 1978): 2
!3 313 13 d Nu ¼ 43:66 þ 0:7 þ 1:077 RePr 0:7 5 z 3
3
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
359
Fig. 2. Radial temperature profiles near the reactor inlet for cooling at constant wall temperature with (a) T0 Tw ¼ 108 K, no chemical reaction and (b) T0 Tw ¼0 K, exothermic reaction with ΔTad ¼ 108 K and EA ¼ 50 kJ/mol.
shown in Fig. 2b, with T0 ¼Tw and a reaction that would cause an adiabatic temperature increase of 108 K. In the case with a reaction taking place, the temperature profiles close to the reactor inlet (e.g., z/L¼0.015) show two radial temperature maxima near the tube wall. Due to the low velocity of the fluid elements close to the reactor wall, higher local conversion gradients, dXA/dz, are reached, which causes an early release of heat from the reaction. As a result, high radial temperature gradients, dT/dr, and subsequently high heat flux densities, q_ w , occur near the inlet. When the resulting temperature distributions (Fig. 2b) are compared to typical laminar temperature profiles, which have a temperature maximum on the tube center line (Fig. 2a), higher heat transfer coefficients can be predicted for the double-peak profiles: Considering distributions with identical temperature differences Tm Tw for both profile shapes, a higher heat transfer _ coefficient, α ¼ ðT mqwT w Þ, is caused by the enhanced heat flux density ∂T q_ w ¼ λ ∂r jr ¼ R that is generated by the reaction. The closer to the inlet, the higher the enhancement of α because the generated heat flux density is high, whereas Tm Tw is still small. For larger values of z/L, the local heat release close to the reactor wall decreases due to the decrease in the reactant concentration. Thus, the two radial temperature maxima move towards the tube center, where the flow velocity is higher, and less reactant is converted. Finally, with high conversions, the maxima merge, and the temperature profiles approach fully developed laminar flow. 3.1.2. Comparison with the plug flow model Fig. 3a and b shows the profiles of the bulk mean temperatures calculated by Eq. (4) from the temperature distributions shown in Fig. 2. For comparison, the dotted lines represent the simplified calculation with the plug flow model. The thermal entrance effects and influences of the heat source were neglected by assuming thermally developed laminar flow with a Nusselt number of 3.66. In both cases, cooling without chemical reaction (Fig. 3a) and cooling with an exothermic reaction (Fig. 3b), the heat transfer is underestimated if thermal effects at the reactor entrance are neglected as in the simplified one dimensional model. However, the thermal entrance region in absence of a heat source (Fig. 3a) can be accounted for by calculating the local Nusselt number with a correlation for the thermal entrance flow (Eq. 10). Fig. 3c shows good agreement of the correlation with local Nusselt numbers calculated from the two-dimensional simulation. Local Nusselt numbers that result from the simulation with an exothermic reaction, on the other hand, indicate an enhancement of the heat transfer that is significantly higher than the enhancement described by the
correlation for the thermal entrance. In conclusion, twodimensional modeling of the reactor is necessary to correctly describe the temperature in the entrance region where the hot spot forms. In their research, Yu and Churchill (2012) found that the enhancement of heat transfer due to chemical reactions is effectively independent of the fractional adiabatic temperature increase ΔT ad =T 0 . This can be explained with regard to the definition of the heat transfer coefficient by the fact that it is proportional to both the heat flux density from the wall and change in the bulk mean temperature. 3.2. Influence of the activation energy The complex behavior of an exothermic reaction with concomitant heat transfer is caused by the temperature dependence of the conversion rate. Therefore, the activation energy of the considered reaction is a crucial parameter in the reactor/heat exchanger model. 3.2.1. Influence on radial concentration profiles The influence of thermal effects on the conversion profiles is demonstrated in Fig. 4. Conversion profiles for a hypothetic activation energy EA ¼0 kJ/mol (which corresponds to the conversion behavior of an isothermal reactor) are compared to the results for an activation energy of EA ¼ 50 kJ/mol (Fig. 4b). Mutual interactions between the temperature and concentration profiles can be observed in Fig. 4b: Near the reactor inlet (z/ L¼ 0.015), the concentration profile is similar to the corresponding profile in Fig. 4a but is clearly influenced by the formation of temperature maxima close to the reactor wall. The radial locations of the temperature maxima are marked by arrows in the concentration profiles. The curve shape for z/L ¼0.03 with three concentration maxima is a result of the temperature peaks moving from the tube wall towards the center. If the temperature peaks merge to one central hot spot (z/L¼0.05), the highest local conversion is reached in the center of the tube. The development of the curve shape in this region is then determined by three parallel effects: 1) The decrease in the flow velocity from the center toward the wall causes higher local residence times of reactant molecules close to the wall. 2) The decrease in the temperature from the center towards the wall causes higher local reaction rates in the center of the tube. Both effects act simultaneously, which explains the formation of concentration maxima between the wall and center line. 3) Furthermore, high concentrations cause higher local reaction rates. Along
360
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
Fig. 3. Comparison of axial temperature profiles (a, b), calculated with the plug flow model (Nu¼ 3.66) and with the two-dimensional model (bulk mean temperatures). Local Nusselt numbers are compared for the cases of cooling with and without an exothermic reaction (c).
Fig. 4. Radial concentration profiles close to the reactor inlet for an exothermic reaction with EA ¼ 0 kJ/mol (a) and EA ¼ 50 kJ/mol (b).
with molecular diffusion, this leads to an equalization of the concentration profiles along the z-coordinate. 3.2.2. Comparison with the plug flow model In Fig. 5, axial profiles of the bulk mean temperature and conversion are compared for different activation energies. With increasing temperature dependence, the conversion profiles vary from quasi-isothermal behavior (EA ¼0 kJ/mol) to runaway behavior (EA ¼80 kJ/mol). The bulk mean values of the two-dimensional model are compared to the simplified one-dimensional model assuming 1) a uniform radial concentration distribution (plug flow) and 2) heat transfer in a thermally fully developed laminar flow with constant wall temperature (Nu¼3.66). Regarding the conversion profile for EA ¼0 kJ/mol, only the first simplification has to be considered. Because thermal effects do not influence the reaction rate, deviations of the conversion profiles are solely caused by the plug flow assumption in the one-dimensional model. Non-uniform concentration profiles, as shown in Fig. 4a, result from the radial velocity distribution and from molecular diffusion. Generally, deviations between both models decrease with enhanced radial mixing, e.g., due to higher diffusion coefficients or smaller diffusion paths in the microchannels. For higher activation energies, increasing the influence of the described thermal effects also changes the conversion and temperature profiles. The lower temperatures in two-dimensional simulation (Fig. 5) result from faster heat removal close to the reactor inlet as well as from lower heat release due to the decreased reaction rate (for further insight, a comparison of heat release rates
and heat transfer rates between both models can be found in Supplementary Material, Fig. S1). As a result, both the extent and location of the hot spot significantly deviate from the results that are obtained from the plug flow model. In the case of runaway behavior (EA ¼80 kJ/mol), the rate of heat removal becomes negligible compared to the rapidly increasing rate of heat generation. Thus, the driving temperature difference at the hot spot Tm,max–Tw approaches a maximum value, which is the adiabatic temperature increase ΔTad. Consequently, the deviations between both models decreases by approaching adiabatic behavior at high activation energies. 3.3. Consequences for reactor design and control of exothermic reactions with cooling Precise evaluation of the formation of hot spots can be critical for reactor safety and product quality. Deviations between the bulk mean temperatures and temperatures calculated with the onedimensional model result in different rates of reaction, possible side reactions or decomposition reactions. Thus, the prediction of the maximum temperature has to be based on an accurate model. For this study, the calculated maximum temperatures given in Table 3 show a maximum deviation of 22 K (EA ¼50 kJ/mol). However, because the heat transfer coefficients are generally underestimated by using the simplified one-dimensional model, it can be suitable for preliminary safety considerations. Additionally, the two-dimensional model allows for a determination of the local maximum temperature in a reactor. This might be important to ensure that the temperature does not locally exceed the boiling or decomposition temperatures.
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
361
Fig. 5. Axial temperature and conversion profiles for different activation energies, calculated with the plug flow model (Nu ¼3.66) and the two-dimensional model (bulk mean values).
Table 3 Maximum temperatures for both employed reactor models and reactor length ratios for a conversion of 90%. Activation energy EA Plug flow model -Maximum temperature Two-dimensional model -Max. bulk mean temperature -Max. local temperature Reactor length ratio for XA ¼0.9 L2d/LPFR
[kJ/mol]
0
30
50
80
[°C]
58
80
105
126
[°C] [°C]
48 62
62 87
83 120
111 155
[-]
1.30
2.57
2.36
1.62
In the example shown in this study, the calculated maximum local temperatures in the radial profiles are up to 44 K above the maximum bulk mean temperature. Differences in the calculated axial conversion profiles, as shown in Fig. 5, can significantly influence the design of a reactor. The reactor length that is required for a conversion of 0.9 was determined with both models. When the radial heat and mass transport is taken into account, generally larger reactor sizes result due to a broader residence time distribution and enhanced heat transfer. In Table 3, these deviations are quantified by the length ratios resulting from both models. For EA ¼ 0 kJ/mol, the influence of thermal effects on the reaction rate is neglected and a ratio of 1.30 basically shows the deviation between isothermal laminar flow and isothermal plug flow. In literature, this value is given with 1.32 for laminar flow without molecular diffusion (Churchill et al., 2008). Thus, molecular diffusion only has small influence on the residence time distribution in the considered scenario. For smaller microchannels, however, the influence will be greater due to the quadratic dependence of the diffusion time on the channel diameter.
If thermal effects are taken into account (EA 40 kJ/mol), the ratio of required reactor length increases significantly. Considering the case accounting for the residence time distribution while neglecting the interactions between heat transfer and reaction, the required reactor volume for XA ¼ 0.9 would still be up to twice as large as the volume calculated with the simplified model. The increase of the length ratio L2d/LPFR mainly depends on the local temperature increase at the hot spot and on the temperature dependency of the exothermic reaction. A preliminary onedimensional simulation can be helpful to estimate the influence of the described thermal effects in the capillary reactor. Analogously to considerations on reactor stability (Semenov, 1961; Baerns et al., 2006; Gelhausen et al., 2015), a maximum dimensionless temperature difference can be calculated with the hot spot temperature in plug flow simulation: EA T max T w ð11Þ Tw RT w Simulation results for varying reaction parameters (see Supplementary Material, Fig. S2) indicate that interactions between heat transfer and reaction can be neglected in good approximation for small values of ΔT 0max (approximately ΔT 0max o 0.5). In this case, one-dimensional models can be applied for predicting the required reactor length and solely the residence time distribution has to be taken into account (e.g. by the convection model for laminar flow (Levenspiel, 1999)). For higher values of ΔT 0max , only the two-dimensional model is able to provide accurate results. ΔT 0max ¼
4. Conclusions For modeling milli- and microstructured reactors, the commonly used plug flow assumption, which neglects radial temperature and conversion profiles, can lead to significant miscalculations. In the
362
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
demonstrated example, high local deviations from the twodimensional model occur, above 30% in conversion and over 20 K in the maximum temperature at the hot spot. Part of these deviations can be explained from the fact that the radial conversion distribution associated with the residence time distribution is neglected. In addition to this non-temperature-related effect, complex thermal behavior has to be considered. With an exothermic reaction, high radial temperature gradients in the inlet region of the reactor can lead to a significant increase in the heat transfer coefficients. Due to the mutual interactions between the local reaction rates and radial temperature and concentration distributions, the thermal behavior strongly depends on the reaction and its activation energy EA. Ignoring these effects in the design of a reactor/heat exchanger would result in the reactor sizes being too small and incorrect predictions of the selectivity, which might strongly depend on local and mean maximum temperatures. Therefore, simulating both the axial and radial conversion profiles can be a key step for reactor design or reactor scale-up. Particularly for the transfer from small microchannels to larger millistructured reactor channels, accurate calculation of the temperature distributions in the hot spot region gives essential information about the reactor safety and scalability. Solving the differential models with conventional simulation tools (e.g., finite element software) is a common practice and can be done with reasonable effort. In practice, where channel designs are often more complex than a straight tube, CFD simulation of the hot spot region close to the reactor inlet is essential for determining the heat and mass transfer rates influenced by an exothermic reaction.
Nomenclature cA cA;0 cp d D EA ΔH R L Nu Pr q_ w r rA R R Re T T max ΔT 0max ΔT ad Tm Tw T0 u um XA z
concentration of component A (mol m 3) initial concentration of component A at the reactor inlet (mol m 3) heat capacity (J kg 1 K 1) tube diameter (m) molecular diffusion coefficient (m2 s 1) activation energy (J mol 1) enthalpy of reaction (J mol 1) reactor length (m) Nusselt number (dimensionless) Prandtl number (dimensionless) wall heat flux density (W m 2) radial distance in wall direction (m) rate of disappearance of A due to reaction (mol m 3 s 1) tube radius (m) universal gas constant (J mol 1 K 1) Reynolds number (dimensionless) temperature (K) maximum temperature (K) maximum dimensionless temperature difference adiabatic temperature increase (K) bulk mean temperature (K) wall temperature (K) temperature at the reactor inlet (K) axial velocity (m s 1) mean axial velocity (m s 1) conversion of component A (dimensionless) axial distance from the reactor inlet (m)
Greek Symbols
α
heat transfer coefficient (W m 2 K 1)
λ ρ
heat conductivity of the fluid (W m 1 K 1) density (kg m 3)
Acknowledgments This work was funded by the German Federal Ministry of Education and Research (BMBF, funding code 03FH012I2). The authors would like to thank Annette Meiners (COMSOL Multiphysics GmbH, Germany) for valuable comments and suggestions.
Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ces.2015.11.022.
References Baehr, H.D., Stephan, K., 2006. Heat and Mass Transfer, 2nd ed. Springer, Berlin, New York. Baerns, M., et al., 2006. Technische Chemie. Wiley-VCH, Weinheim. Barthe, P., Guermeur, C., Lobet, O., Moreno, M., Woehl, P., Roberge, D.M., Bieler, N., Zimmermann, B., 2008. Continuous multi-injection reactor for multipurpose production-Part I. Chem. Eng. Technol. 31 (8), 1146–1154. Boskovic, D., Loebecke, S., 2008. Modelling of the residence time distribution in micromixers. Chem. Eng. J. 135, S138–S146. Churchill, S.W., 2005. Interaction of chemical reactions and transport. 1. An overview. Ind. Eng. Chem. Res. 44 (14), 5199–5212. Churchill, S.W., Yu, B., Kawaguchi, Y., 2008. Interaction of reactions and transport: homogeneous reactions in a round tube. IV. The first-order effects of flow and heat of reaction. Ind. Eng. Chem. Res. 47 (15), 5132–5145. Gelhausen, M.G., Kurt, S.K., Kockmann, N., 2015. Parametrische Empfindlichkeit einer stark exothermen Reaktion im Kapillarwendelreaktor. Chem. Ing. Tech. 87 (6), 781–790. Gnielinski, V., 1989. Zur Wärmeübertragung bei laminarer Rohrströmung und konstanter Wandtemperatur. Chem. Ing. Tech. 61 (2), 160–161. Klutz, S., Kurt, S.K., Lobedann, M., Kockmann, N., 2015. Narrow residence time distribution in tubular reactor concept for Reynolds number range of 10–100. Chem. Eng. Res. Des. 95, 22–33. Kockmann, N., Roberge, D.M., 2009. Harsh reaction conditions in continuous-flow microreactors for pharmaceutical production. Chem. Eng. Technol. 32 (11), 1682–1694. Levenspiel, O., 1999. Chemical Reaction Engineering, 3rd ed. Wiley, New York. Lohse, S., Kohnen, B.T., Janasek, D., Dittrich, P.S., Franzke, J., Agar, D.W., 2008. A novel method for determining residence time distribution in intricately structured microreactors. Lab. Chip 8 (3), 431. Martin, H., 2002. The generalized Lévêque equation and its practical use for the prediction of heat and mass transfer rates from pressure drop. Chem. Eng. Sci. 57 (16), 3217–3223. Merrill, L.S., Hamrin, C.E., 1970. Conversion and temperature profiles for complex reactions in laminar and plug flow. AIChE J. 16 (2), 194–198. Rothenberg, R.I., Smith, J.M., 1966. Heat transfer and reaction in laminar tube flow. AIChE J. 12 (2), 213–220. Semenov, N.V., 1961. Einige Probleme der chemischen Kinetik und Reaktionsfähigkeit. Akademie Verlag, Berlin. Shah, R.K., London, A.L., 1978. Laminar Flow Forced Convection in Ducts: A Source Book for Compact Heat Exchanger Analytical Data. Academic Press, New York, San Francisco, London. Taghavi-Moghadam, S., Kleemann, A., Golbig, G., 2001. Microreaction Technology as a Novel Approach to Drug Design, Process Development and Reliability. Org. Process. Res. Dev. 5 (6), 652–658. Trachsel, F., Günther, A., Khan, S., Jensen, K.F., 2005. Measurement of residence time distribution in microfluidic systems. Chem. Eng. Sci. 60 (21), 5729–5737. Wörner, M., Ghidersa, B., Onea, A., 2007. A model for the residence time distribution of bubble-train flow in a square mini-channel based on direct numerical simulation results. Int. J. Heat. Fluid Flow 28 (1), 83–94. Yu, B., Churchill, S.W., 2012. Prediction of the influence of energetic chemical reactions on forced convective heat transfer, 1st ed. In: Cho, Y.I., Greene, G.A. (Eds.), Advances in Heat Transfer, vol. 44. Academic, Boston, pp. 1–117. Zhang, X., Stefanick, S., Villani, F.J., 2004. Application of microreactor technology in process development. Org. Process. Res. Dev. 8 (3), 455–460.
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces
Enhanced heat transfer by exothermic reactions in laminar flow capillary reactors Sebastian Schwolow a, Jing Ying Ko a, Norbert Kockmann b, Thorsten Röder a,n a b
Mannheim University of Applied Sciences, Institute of Chemical Process Engineering, Paul-Wittsack-Straße 10, 68163 Mannheim, Germany TU Dortmund, Biochemical and Chemical Engineering, Equipment Design, Emil-Figge-Straße 68, 44227 Dortmund, Germany
H I G H L I G H T S
G R A P H I C A L
A B S T R A C T
For a capillary reactor, concentration and temperature profiles are discussed. Exothermic reactions can significantly enhance heat transfer. For reactor design, commonly used plug flow models can lead to miscalculations. Simulating both axial and radial conversion profiles can be a key step for reactor scale-up.
art ic l e i nf o
a b s t r a c t
Article history: Received 15 September 2015 Received in revised form 20 November 2015 Accepted 21 November 2015 Available online 29 November 2015
Conversion and temperature distributions in a tubular reactor with an inner diameter of 1 mm were numerically calculated for a second order exothermic reaction in laminar flow of homogeneous liquid. Based on the resulting radial temperature profiles, the local Nusselt numbers and bulk mean temperatures were determined along the reactor tube. Strong effects of the homogeneous reaction on the heat transfer can be observed in the entrance region of the reactor, where a hot spot emerges. Due to the large radial temperature gradients in vicinity of the reactor wall, heat transfer coefficients are significantly higher compared to a non-reactive system. The consequences of this effect on the design and control of exothermic reactions in reactor/heat exchangers are demonstrated by comparison with a simple one-dimensional plug flow model. In the simplified model, neglecting thermal influence of the exothermic reaction results in a significant underestimation of the required reactor length for defined conversion. Accordingly, numerical simulation of both axial and radial transport in the hot spot region can be essential to precisely predict the bulk temperatures and conversion rates in the reaction mixture, even with the small length scales of milli- and micro-structured reactors. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Capillary reactor Heat transfer Simulation Nusselt number Temperature distribution Entrance flow
1. Introduction Milli- and micro-structured flow reactors have a high potential for process intensification. The finely adjustable reaction times and excellent heat removal in these reactors make it possible to run exothermic reactions in safe and defined conditions (Kockmann and Roberge, 2009; Taghavi-Moghadam et al., 2001; Zhang et al., 2004). To calculate the conversion profiles in microreactors, plug
n
Corresponding author. Tel.: þ 49 621 292 6800. E-mail address: [email protected] (T. Röder).
http://dx.doi.org/10.1016/j.ces.2015.11.022 0009-2509/& 2015 Elsevier Ltd. All rights reserved.
flow and often isothermal assumptions are widely used. Both simplifications are based on the intensified transport processes with the small channel diameters and large reactor length-todiameter ratios: 1) With very short diffusion paths, radial mixing by molecular diffusion can be fast enough to approach plug flow behavior, even in the laminar flow regime. Further enhancement of radial mixing can be obtained by inducing a secondary flow with the specific channel design (Klutz et al., 2015). 2) Heat transfer within the reactor is very efficient as a result of the high surface-to-volume ratio and small length scale for thermal conduction. Thus, the rate of heat removal is assumed to be fast compared to the rate of heat generation and, as a result, the
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
reaction mixture may be considered to be nearly isothermal. However, at a closer look, these common idealizations can significantly deviate from the actual reactor behavior: If the hydrodynamic residence time is much shorter than the characteristic time for molecular diffusion, radial concentration gradients exist, and either the convection model for laminar flow or the dispersion model are more accurate representations of the residence time distribution (Levenspiel, 1999). In the case of very fast and strongly exothermic reactions, even the superior heat transfer characteristics of microreactors might not be sufficient to avoid hot spots in the entrance region of the reactor (Barthe et al., 2008). Local temperature peaks accelerate the reaction and significantly influence the yield in the reaction system as a result. While numerous publications describe experimental characterization and modeling of the residence time distribution in microreactors (Boskovic and Loebecke, 2008; Lohse et al., 2008; Trachsel et al., 2005; Wörner et al., 2007), less consideration is given to the temperature distribution in micro- and millireactors. Nevertheless, the formation of hot spots is a crucial concern for reactor design and control. Because temperature measurements in small spaces of microreactors are often difficult to perform in practice, precise calculation of local temperatures has additional significance. For a rough estimation, the internal heat transfer coefficient, α, can be determined from correlations for the dimensionless Nusselt number, Nu ¼ αλd, in the channel flow (Martin, 2002). For a straight tubular reactor (e.g., a stainless steel capillary in an oil bath), the internal heat transfer can be calculated by assuming a constant wall temperature and fully developed flow. With these assumptions, a constant Nusselt number, Nu ¼ 3:6568, can be derived for laminar flow (Baehr and Stephan, 2006). However, radial temperature profiles can be strongly influenced by the reaction itself. This effect is completely neglected by using common correlations or assuming a constant Nusselt number throughout the entire channel. An intensive and comprehensive insight into the phenomena resulting from an exo- or endothermic chemical reaction with concomitant heat transfer has been published by Churchill et al. (Churchill, 2005; Churchill et al., 2008; Yu and Churchill, 2012). Churchill describes strong effects and complex behavior that significantly enhance heat transfer in the presence of a homogeneous exothermic reaction. Because additional parameters are introduced by one or several chemical reactions, he concluded that a general correlating equation cannot be derived practically. Instead, he recommends that numerical calculation of both the radial and axial transport be performed on a case-by-case basis. Radial temperature profiles in tubular reactors have been investigated for some typical gaseous reaction systems, for example, by Rothenberg and Smith (1966), or by Merrill and Hamrin (1970). The objective of the present study is the detailed investigation and discussion of heat transfer in a liquid phase reaction, based on exemplary process specifications that closely approximate the typical operating conditions in milli- and microstructured flow reactors. A fast, exothermic, second-order reaction in a circular duct with 1 mm diameter was simulated using two different models: 1) A simple one-dimensional plug flow model assuming a constant Nusselt number (heat transfer in thermally developed flow) and 2) a two-dimensional model of the rotationally symmetric geometry that accounts for the axial and radial heat and mass transport. Visualization and discussion of radial temperature and concentration profiles in the investigated capillary reactor provide the basis for an improved understanding of the effects previously described by Churchill and coworkers. By comparing results of both models, the applicability of the simplified model can be evaluated with regard to practical aspects in the design and control of reactor/heat exchangers. Furthermore, the influence of the exothermic reaction is demonstrated by comparing the reactor/heat exchanger with a corresponding heat exchanger without chemical reaction.
357
2. Basic model assumptions The exothermic reaction is assumed to be second order and irreversible and to take place in a homogeneous liquid phase (Table 1). The temperature dependence of the reaction rate was described by the Arrhenius equation, in which the activation energy was varied to ascertain the thermal effects. The material properties used were based on toluene and assumed to be constant for all temperatures and conversions. For the purpose of this study, a reaction medium with a density of 867 kg m 3, heat capacity of 1707 J kg 1 K 1, thermal conductivity of 0.134 W m 1 K 1, and a mass diffusivity of 10 9 m2 s 1 for the reactant was considered. The considered reactor is a tubular flow reactor (Table 2). Due to its rotational symmetry, the three-dimensional problem can be reduced to a two-dimensional axisymmetric problem. A constant wall temperature, Tw, was defined and for all simulations with an exothermic reaction, the inlet temperature of the reactor was defined as T0 ¼Tw. The partial differential equations resulting from the mass and energy balance in the tube flow are ∂X A D 1 ∂X A ∂2 X A ∂2 X A rA þ 2 þ 2 ¼0 ∂z u r ∂r ucA;0 ∂r ∂z
ð1Þ
∂T λ 1 ∂T ∂2 T ∂2 T r A ΔHR þ 2þ 2 þ ¼0 ∂z uρcp r ∂r ∂r uρcp ∂z
ð2Þ
The following boundary conditions that were used in this study: T ðr; z ¼ 0Þ ¼ T 0
X A ðr; z ¼ 0Þ ¼ 0
T ðr ¼ R; zÞ ¼ T w
dX A ðr ¼ R; zÞ ¼ 0 dr
dT ðr ¼ 0; zÞ ¼ 0 dr
dX A ðr ¼ 0; zÞ ¼ 0 dr
For all z, the velocity profile in laminar Hagen–Poiseuille flow is defined as r 2 uðr Þ ¼ 2um 1 R
ð3Þ
With these definitions, it was assumed that the reaction was initiated in a hydrodynamic fully developed flow with a homogeneous radial reactant distribution. Table 1 Reaction parameters used for the simulations. Rate constant Reference temperature Initial reactant concentration Reaction enthalpy
0.0004 30 4 40
m3 mol 1 s 1 °C mol L 1 kJ mol 1
Table 2 Reactor parameters used for the simulations. Reactor diameter Reactor length Mean flow velocity Wall temperature Reactor pressurea
1 1 0.1 25 5
mm m m s1 °C bar
a The reactor pressure was chosen in order to avoid exceeding the boiling temperature in the reactor tube for the reference material toluene.
358
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
Fig. 1. Temperature and conversion profiles in the entrance region of the reactor, with an exothermic second order reaction (ΔTad ¼ 108 K and EA ¼ 50 kJ/mol).
To determine the longitudinal temperature profile in the tube from the results, fluid bulk mean temperatures were calculated: Z R 1 Tm ¼ 2 2r πuðr ÞT ðr Þdr ð4Þ R πum 0 This temperature, also referred to as “mixing cup temperature”, is the bulk temperature that theoretically would be measured if the channel was sliced at any cross-section and the outflowing fluid was captured in an adiabatic, well-mixed cup. Bulk mean conversions were determined in an analogous manner. The local Nusselt number in the reactor can be obtained, when the heat transfer coefficient is calculated from the wall heat flux density divided by Tm Tw:
growth rate of 1.031. Thus, a very fine mesh resolution could be obtained in the wall region with a minimum element width of 0.7 mm.
3. Results and discussion 3.1. Heat transfer enhancement
ð10Þ
In the absence of thermal energy sources, the temperature distribution in a fully developed Hagen-Poiseuille flow can be calculated based on the energy balance. A derivation of the temperature distribution for constant wall temperature can be found in the literature (Baehr and Stephan, 2006) and results in a rapidly convergent power series. The fully developed Nusselt number corresponding to this temperature distribution was shown to be Nu¼ 3.66. For the case of a thermal entrance in hydrodynamically developed flow, the local Nusselt numbers in the thermal entrance region are higher compared to those in the fully developed flow region, and can be calculated with correlations given in the literature (Gnielinski, 1989; Shah and London, 1978). With an exothermic chemical reaction, a thermal energy source is added to the heat transfer problem. Due to the dependence of the heat release on the local reactant concentration and temperature, radial and axial transport has to be calculated numerically. The results are shown in Fig. 1, where the development of the radial temperature and conversion profiles in the entrance region is pictured. Herein, fluid enters the tube with the wall temperature (T0 ¼Tw), and the internal heat source in the fluid is determined from the second order reaction equations. On a longitudinal length scale of z 50 Ud, a hot spot is formed. The reaction rate is high near the reactor inlet because of the high reactant concentrations. As the fluid moves forward in the reactor channel, the heat generated by the reaction decreases with the increase in reactant conversion. Simultaneously, the heat removal rate increases with the temperature difference Tm Tw. The temperature and conversion profiles shown in Fig. 1 exhibit thermal effects with mutual interactions. The formation of the hot spot enhances the reaction rate, therefore, a fast increase in local conversion can be observed at the hot spot.
The differential equations of the plug flow model were solved using the ode15s-solver in Matlab (Mathworks Inc., Natick, MA, USA). The two-dimensional axisymmetric problem was solved with the finite element software COMSOL Multiphysics (COMSOL, Inc., Palo Alto, CA, USA). Herein, a structured quadrilateral mesh with a 720,000 elements was used. The highest gradients can be expected at the reactor inlet and close to the reactor wall. For this reason, mesh refinement was applied with element size distributions defined by arithmetic sequences in the z- and r-directions, with an average
3.1.1. Influence of the exothermic reaction on radial temperature profiles To discuss the thermal effects, Fig. 2 shows a comparison between cooling a non-reactive hot inlet stream (2 a) and cooling a fluid which is increasing in temperature due to an exothermic reaction (2b). In Fig. 2a, the selected radial temperature profiles show the effect of the thermal entrance without a chemical reaction and T0 Tw ¼ 108 K. Corresponding temperature profiles from the simulation that accounts for the exothermic reaction are
Nu ¼
q_ w 2R αd ¼ λ ðT m T w Þλ
ð5Þ
Herein, qw was determined from the enthalpy balance in a cylindrical tube element 2 dT m q_ w ¼ um ρcp r A ΔH R R dz
ð6Þ
For the simplified, one-dimensional model, the usual assumptions of plug flow were applied, which reduce Eqs. (1) and (2) to dX A rA ¼0 dz cA;0 u
ð7Þ
dT rA ΔH R þ q_ w R2 þ ¼0 dz uρcp
ð8Þ
In Eq. (8), q_ w R2 is the heat transfer rate to the wall over the reactor length with :
qw ¼ α ðT T w Þ
ð9Þ
The heat transfer coefficient α was defined by assuming fully developed, laminar flow heat transfer with Nu ¼ 3.66. For the case of heat transfer without chemical reactions, the local Nusselt numbers in the thermal entrance region can be calculated by the following correlation (Shah and London, 1978): 2
!3 313 13 d Nu ¼ 43:66 þ 0:7 þ 1:077 RePr 0:7 5 z 3
3
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
359
Fig. 2. Radial temperature profiles near the reactor inlet for cooling at constant wall temperature with (a) T0 Tw ¼ 108 K, no chemical reaction and (b) T0 Tw ¼0 K, exothermic reaction with ΔTad ¼ 108 K and EA ¼ 50 kJ/mol.
shown in Fig. 2b, with T0 ¼Tw and a reaction that would cause an adiabatic temperature increase of 108 K. In the case with a reaction taking place, the temperature profiles close to the reactor inlet (e.g., z/L¼0.015) show two radial temperature maxima near the tube wall. Due to the low velocity of the fluid elements close to the reactor wall, higher local conversion gradients, dXA/dz, are reached, which causes an early release of heat from the reaction. As a result, high radial temperature gradients, dT/dr, and subsequently high heat flux densities, q_ w , occur near the inlet. When the resulting temperature distributions (Fig. 2b) are compared to typical laminar temperature profiles, which have a temperature maximum on the tube center line (Fig. 2a), higher heat transfer coefficients can be predicted for the double-peak profiles: Considering distributions with identical temperature differences Tm Tw for both profile shapes, a higher heat transfer _ coefficient, α ¼ ðT mqwT w Þ, is caused by the enhanced heat flux density ∂T q_ w ¼ λ ∂r jr ¼ R that is generated by the reaction. The closer to the inlet, the higher the enhancement of α because the generated heat flux density is high, whereas Tm Tw is still small. For larger values of z/L, the local heat release close to the reactor wall decreases due to the decrease in the reactant concentration. Thus, the two radial temperature maxima move towards the tube center, where the flow velocity is higher, and less reactant is converted. Finally, with high conversions, the maxima merge, and the temperature profiles approach fully developed laminar flow. 3.1.2. Comparison with the plug flow model Fig. 3a and b shows the profiles of the bulk mean temperatures calculated by Eq. (4) from the temperature distributions shown in Fig. 2. For comparison, the dotted lines represent the simplified calculation with the plug flow model. The thermal entrance effects and influences of the heat source were neglected by assuming thermally developed laminar flow with a Nusselt number of 3.66. In both cases, cooling without chemical reaction (Fig. 3a) and cooling with an exothermic reaction (Fig. 3b), the heat transfer is underestimated if thermal effects at the reactor entrance are neglected as in the simplified one dimensional model. However, the thermal entrance region in absence of a heat source (Fig. 3a) can be accounted for by calculating the local Nusselt number with a correlation for the thermal entrance flow (Eq. 10). Fig. 3c shows good agreement of the correlation with local Nusselt numbers calculated from the two-dimensional simulation. Local Nusselt numbers that result from the simulation with an exothermic reaction, on the other hand, indicate an enhancement of the heat transfer that is significantly higher than the enhancement described by the
correlation for the thermal entrance. In conclusion, twodimensional modeling of the reactor is necessary to correctly describe the temperature in the entrance region where the hot spot forms. In their research, Yu and Churchill (2012) found that the enhancement of heat transfer due to chemical reactions is effectively independent of the fractional adiabatic temperature increase ΔT ad =T 0 . This can be explained with regard to the definition of the heat transfer coefficient by the fact that it is proportional to both the heat flux density from the wall and change in the bulk mean temperature. 3.2. Influence of the activation energy The complex behavior of an exothermic reaction with concomitant heat transfer is caused by the temperature dependence of the conversion rate. Therefore, the activation energy of the considered reaction is a crucial parameter in the reactor/heat exchanger model. 3.2.1. Influence on radial concentration profiles The influence of thermal effects on the conversion profiles is demonstrated in Fig. 4. Conversion profiles for a hypothetic activation energy EA ¼0 kJ/mol (which corresponds to the conversion behavior of an isothermal reactor) are compared to the results for an activation energy of EA ¼ 50 kJ/mol (Fig. 4b). Mutual interactions between the temperature and concentration profiles can be observed in Fig. 4b: Near the reactor inlet (z/ L¼ 0.015), the concentration profile is similar to the corresponding profile in Fig. 4a but is clearly influenced by the formation of temperature maxima close to the reactor wall. The radial locations of the temperature maxima are marked by arrows in the concentration profiles. The curve shape for z/L ¼0.03 with three concentration maxima is a result of the temperature peaks moving from the tube wall towards the center. If the temperature peaks merge to one central hot spot (z/L¼0.05), the highest local conversion is reached in the center of the tube. The development of the curve shape in this region is then determined by three parallel effects: 1) The decrease in the flow velocity from the center toward the wall causes higher local residence times of reactant molecules close to the wall. 2) The decrease in the temperature from the center towards the wall causes higher local reaction rates in the center of the tube. Both effects act simultaneously, which explains the formation of concentration maxima between the wall and center line. 3) Furthermore, high concentrations cause higher local reaction rates. Along
360
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
Fig. 3. Comparison of axial temperature profiles (a, b), calculated with the plug flow model (Nu¼ 3.66) and with the two-dimensional model (bulk mean temperatures). Local Nusselt numbers are compared for the cases of cooling with and without an exothermic reaction (c).
Fig. 4. Radial concentration profiles close to the reactor inlet for an exothermic reaction with EA ¼ 0 kJ/mol (a) and EA ¼ 50 kJ/mol (b).
with molecular diffusion, this leads to an equalization of the concentration profiles along the z-coordinate. 3.2.2. Comparison with the plug flow model In Fig. 5, axial profiles of the bulk mean temperature and conversion are compared for different activation energies. With increasing temperature dependence, the conversion profiles vary from quasi-isothermal behavior (EA ¼0 kJ/mol) to runaway behavior (EA ¼80 kJ/mol). The bulk mean values of the two-dimensional model are compared to the simplified one-dimensional model assuming 1) a uniform radial concentration distribution (plug flow) and 2) heat transfer in a thermally fully developed laminar flow with constant wall temperature (Nu¼3.66). Regarding the conversion profile for EA ¼0 kJ/mol, only the first simplification has to be considered. Because thermal effects do not influence the reaction rate, deviations of the conversion profiles are solely caused by the plug flow assumption in the one-dimensional model. Non-uniform concentration profiles, as shown in Fig. 4a, result from the radial velocity distribution and from molecular diffusion. Generally, deviations between both models decrease with enhanced radial mixing, e.g., due to higher diffusion coefficients or smaller diffusion paths in the microchannels. For higher activation energies, increasing the influence of the described thermal effects also changes the conversion and temperature profiles. The lower temperatures in two-dimensional simulation (Fig. 5) result from faster heat removal close to the reactor inlet as well as from lower heat release due to the decreased reaction rate (for further insight, a comparison of heat release rates
and heat transfer rates between both models can be found in Supplementary Material, Fig. S1). As a result, both the extent and location of the hot spot significantly deviate from the results that are obtained from the plug flow model. In the case of runaway behavior (EA ¼80 kJ/mol), the rate of heat removal becomes negligible compared to the rapidly increasing rate of heat generation. Thus, the driving temperature difference at the hot spot Tm,max–Tw approaches a maximum value, which is the adiabatic temperature increase ΔTad. Consequently, the deviations between both models decreases by approaching adiabatic behavior at high activation energies. 3.3. Consequences for reactor design and control of exothermic reactions with cooling Precise evaluation of the formation of hot spots can be critical for reactor safety and product quality. Deviations between the bulk mean temperatures and temperatures calculated with the onedimensional model result in different rates of reaction, possible side reactions or decomposition reactions. Thus, the prediction of the maximum temperature has to be based on an accurate model. For this study, the calculated maximum temperatures given in Table 3 show a maximum deviation of 22 K (EA ¼50 kJ/mol). However, because the heat transfer coefficients are generally underestimated by using the simplified one-dimensional model, it can be suitable for preliminary safety considerations. Additionally, the two-dimensional model allows for a determination of the local maximum temperature in a reactor. This might be important to ensure that the temperature does not locally exceed the boiling or decomposition temperatures.
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
361
Fig. 5. Axial temperature and conversion profiles for different activation energies, calculated with the plug flow model (Nu ¼3.66) and the two-dimensional model (bulk mean values).
Table 3 Maximum temperatures for both employed reactor models and reactor length ratios for a conversion of 90%. Activation energy EA Plug flow model -Maximum temperature Two-dimensional model -Max. bulk mean temperature -Max. local temperature Reactor length ratio for XA ¼0.9 L2d/LPFR
[kJ/mol]
0
30
50
80
[°C]
58
80
105
126
[°C] [°C]
48 62
62 87
83 120
111 155
[-]
1.30
2.57
2.36
1.62
In the example shown in this study, the calculated maximum local temperatures in the radial profiles are up to 44 K above the maximum bulk mean temperature. Differences in the calculated axial conversion profiles, as shown in Fig. 5, can significantly influence the design of a reactor. The reactor length that is required for a conversion of 0.9 was determined with both models. When the radial heat and mass transport is taken into account, generally larger reactor sizes result due to a broader residence time distribution and enhanced heat transfer. In Table 3, these deviations are quantified by the length ratios resulting from both models. For EA ¼ 0 kJ/mol, the influence of thermal effects on the reaction rate is neglected and a ratio of 1.30 basically shows the deviation between isothermal laminar flow and isothermal plug flow. In literature, this value is given with 1.32 for laminar flow without molecular diffusion (Churchill et al., 2008). Thus, molecular diffusion only has small influence on the residence time distribution in the considered scenario. For smaller microchannels, however, the influence will be greater due to the quadratic dependence of the diffusion time on the channel diameter.
If thermal effects are taken into account (EA 40 kJ/mol), the ratio of required reactor length increases significantly. Considering the case accounting for the residence time distribution while neglecting the interactions between heat transfer and reaction, the required reactor volume for XA ¼ 0.9 would still be up to twice as large as the volume calculated with the simplified model. The increase of the length ratio L2d/LPFR mainly depends on the local temperature increase at the hot spot and on the temperature dependency of the exothermic reaction. A preliminary onedimensional simulation can be helpful to estimate the influence of the described thermal effects in the capillary reactor. Analogously to considerations on reactor stability (Semenov, 1961; Baerns et al., 2006; Gelhausen et al., 2015), a maximum dimensionless temperature difference can be calculated with the hot spot temperature in plug flow simulation: EA T max T w ð11Þ Tw RT w Simulation results for varying reaction parameters (see Supplementary Material, Fig. S2) indicate that interactions between heat transfer and reaction can be neglected in good approximation for small values of ΔT 0max (approximately ΔT 0max o 0.5). In this case, one-dimensional models can be applied for predicting the required reactor length and solely the residence time distribution has to be taken into account (e.g. by the convection model for laminar flow (Levenspiel, 1999)). For higher values of ΔT 0max , only the two-dimensional model is able to provide accurate results. ΔT 0max ¼
4. Conclusions For modeling milli- and microstructured reactors, the commonly used plug flow assumption, which neglects radial temperature and conversion profiles, can lead to significant miscalculations. In the
362
S. Schwolow et al. / Chemical Engineering Science 141 (2016) 356–362
demonstrated example, high local deviations from the twodimensional model occur, above 30% in conversion and over 20 K in the maximum temperature at the hot spot. Part of these deviations can be explained from the fact that the radial conversion distribution associated with the residence time distribution is neglected. In addition to this non-temperature-related effect, complex thermal behavior has to be considered. With an exothermic reaction, high radial temperature gradients in the inlet region of the reactor can lead to a significant increase in the heat transfer coefficients. Due to the mutual interactions between the local reaction rates and radial temperature and concentration distributions, the thermal behavior strongly depends on the reaction and its activation energy EA. Ignoring these effects in the design of a reactor/heat exchanger would result in the reactor sizes being too small and incorrect predictions of the selectivity, which might strongly depend on local and mean maximum temperatures. Therefore, simulating both the axial and radial conversion profiles can be a key step for reactor design or reactor scale-up. Particularly for the transfer from small microchannels to larger millistructured reactor channels, accurate calculation of the temperature distributions in the hot spot region gives essential information about the reactor safety and scalability. Solving the differential models with conventional simulation tools (e.g., finite element software) is a common practice and can be done with reasonable effort. In practice, where channel designs are often more complex than a straight tube, CFD simulation of the hot spot region close to the reactor inlet is essential for determining the heat and mass transfer rates influenced by an exothermic reaction.
Nomenclature cA cA;0 cp d D EA ΔH R L Nu Pr q_ w r rA R R Re T T max ΔT 0max ΔT ad Tm Tw T0 u um XA z
concentration of component A (mol m 3) initial concentration of component A at the reactor inlet (mol m 3) heat capacity (J kg 1 K 1) tube diameter (m) molecular diffusion coefficient (m2 s 1) activation energy (J mol 1) enthalpy of reaction (J mol 1) reactor length (m) Nusselt number (dimensionless) Prandtl number (dimensionless) wall heat flux density (W m 2) radial distance in wall direction (m) rate of disappearance of A due to reaction (mol m 3 s 1) tube radius (m) universal gas constant (J mol 1 K 1) Reynolds number (dimensionless) temperature (K) maximum temperature (K) maximum dimensionless temperature difference adiabatic temperature increase (K) bulk mean temperature (K) wall temperature (K) temperature at the reactor inlet (K) axial velocity (m s 1) mean axial velocity (m s 1) conversion of component A (dimensionless) axial distance from the reactor inlet (m)
Greek Symbols
α
heat transfer coefficient (W m 2 K 1)
λ ρ
heat conductivity of the fluid (W m 1 K 1) density (kg m 3)
Acknowledgments This work was funded by the German Federal Ministry of Education and Research (BMBF, funding code 03FH012I2). The authors would like to thank Annette Meiners (COMSOL Multiphysics GmbH, Germany) for valuable comments and suggestions.
Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ces.2015.11.022.
References Baehr, H.D., Stephan, K., 2006. Heat and Mass Transfer, 2nd ed. Springer, Berlin, New York. Baerns, M., et al., 2006. Technische Chemie. Wiley-VCH, Weinheim. Barthe, P., Guermeur, C., Lobet, O., Moreno, M., Woehl, P., Roberge, D.M., Bieler, N., Zimmermann, B., 2008. Continuous multi-injection reactor for multipurpose production-Part I. Chem. Eng. Technol. 31 (8), 1146–1154. Boskovic, D., Loebecke, S., 2008. Modelling of the residence time distribution in micromixers. Chem. Eng. J. 135, S138–S146. Churchill, S.W., 2005. Interaction of chemical reactions and transport. 1. An overview. Ind. Eng. Chem. Res. 44 (14), 5199–5212. Churchill, S.W., Yu, B., Kawaguchi, Y., 2008. Interaction of reactions and transport: homogeneous reactions in a round tube. IV. The first-order effects of flow and heat of reaction. Ind. Eng. Chem. Res. 47 (15), 5132–5145. Gelhausen, M.G., Kurt, S.K., Kockmann, N., 2015. Parametrische Empfindlichkeit einer stark exothermen Reaktion im Kapillarwendelreaktor. Chem. Ing. Tech. 87 (6), 781–790. Gnielinski, V., 1989. Zur Wärmeübertragung bei laminarer Rohrströmung und konstanter Wandtemperatur. Chem. Ing. Tech. 61 (2), 160–161. Klutz, S., Kurt, S.K., Lobedann, M., Kockmann, N., 2015. Narrow residence time distribution in tubular reactor concept for Reynolds number range of 10–100. Chem. Eng. Res. Des. 95, 22–33. Kockmann, N., Roberge, D.M., 2009. Harsh reaction conditions in continuous-flow microreactors for pharmaceutical production. Chem. Eng. Technol. 32 (11), 1682–1694. Levenspiel, O., 1999. Chemical Reaction Engineering, 3rd ed. Wiley, New York. Lohse, S., Kohnen, B.T., Janasek, D., Dittrich, P.S., Franzke, J., Agar, D.W., 2008. A novel method for determining residence time distribution in intricately structured microreactors. Lab. Chip 8 (3), 431. Martin, H., 2002. The generalized Lévêque equation and its practical use for the prediction of heat and mass transfer rates from pressure drop. Chem. Eng. Sci. 57 (16), 3217–3223. Merrill, L.S., Hamrin, C.E., 1970. Conversion and temperature profiles for complex reactions in laminar and plug flow. AIChE J. 16 (2), 194–198. Rothenberg, R.I., Smith, J.M., 1966. Heat transfer and reaction in laminar tube flow. AIChE J. 12 (2), 213–220. Semenov, N.V., 1961. Einige Probleme der chemischen Kinetik und Reaktionsfähigkeit. Akademie Verlag, Berlin. Shah, R.K., London, A.L., 1978. Laminar Flow Forced Convection in Ducts: A Source Book for Compact Heat Exchanger Analytical Data. Academic Press, New York, San Francisco, London. Taghavi-Moghadam, S., Kleemann, A., Golbig, G., 2001. Microreaction Technology as a Novel Approach to Drug Design, Process Development and Reliability. Org. Process. Res. Dev. 5 (6), 652–658. Trachsel, F., Günther, A., Khan, S., Jensen, K.F., 2005. Measurement of residence time distribution in microfluidic systems. Chem. Eng. Sci. 60 (21), 5729–5737. Wörner, M., Ghidersa, B., Onea, A., 2007. A model for the residence time distribution of bubble-train flow in a square mini-channel based on direct numerical simulation results. Int. J. Heat. Fluid Flow 28 (1), 83–94. Yu, B., Churchill, S.W., 2012. Prediction of the influence of energetic chemical reactions on forced convective heat transfer, 1st ed. In: Cho, Y.I., Greene, G.A. (Eds.), Advances in Heat Transfer, vol. 44. Academic, Boston, pp. 1–117. Zhang, X., Stefanick, S., Villani, F.J., 2004. Application of microreactor technology in process development. Org. Process. Res. Dev. 8 (3), 455–460.