Contact Modifications for CFD Simulations of Fixed-Bed Reactors

Article pubs.acs.org/IECR

Contact Modifications for CFD Simulations of Fixed-Bed Reactors: Cylindrical Particles Gregor D. Wehinger,* Carsten Fütterer, and Matthias Kraume Chemical and Process Engineering, Technische Universität Berlin, Fraunhoferstr. 33-36, 10587 Berlin, Germany ABSTRACT: One of the major issues for particle-resolved CFD simulations of fixed beds are particle−particle and particle−wall contacts. A fixed bed of cylinders is studied with different local contact modifications (i.e., caps and bridges method for line and area contacts, and caps and united method for overlaps resulting from composite DEM particles). Effects are analyzed regarding pressure drop, local velocity, and local heat transfer for 191 < Rep < 763. Modeling particle−wall contacts is becoming more important with increasing Rep. Contact modifications inside the bed are becoming less important because convective heat transfer is dominant. The stable caps method shows promising results in comparison with heat-transfer predictions. However, it overestimates convective heat transfer in the contact regions for high flow rates. The bridges method shows good results for pressure drop and heat-transfer predictions. A difficulty remains the choice of the thermal conductivity of these bridges.

1. INTRODUCTION Fixed-bed reactors are commonly used in the chemical and process industry.1 Especially for highly endothermic (or exothermic) reactions, the requirements for the design of fixedbed reactors are often somehow special. The choice of tube diameter (D) and particle dimensions (i.e., particle diameter dP and particle length lP) have several constraints. Heat has to be transferred efficiently into the reactor, which leads to small tube diameters. A reasonable pressure drop and a high gas throughput constrain the minimum particle diameter. Thus, a small tube-toparticle-diameter ratio (1 < N = D/dP < 15) is recommended.2 For these special fixed-bed reactor configurations, conventional descriptions based on plug-flow and pseudohomogeneous models are questionable. Local backflow and stagnation zones occur; channeling in the near-wall region affects the radial heat transfer; radiation plays a significant role at the present temperatures; local kinetics influence local transport phenomena, and vice versa; reaction rates are limited due to film diffusion or pore processes. As a consequence, simplified flow and kinetic models fail to predict severe situations in these multiscale, multiphenomena reactors. Over the last two decades, computational fluid dynamics (CFD) has been used to model the actual flow field in such low N packed beds. One of the pioneering works was the contribution of Derkx and Dixon.3 With their short paper on modeling with CFD flow between two spheres, the idea of resolving the interstitial fluid flow inside a fixed-bed reactor was realized for the first time. The concept is to take into account the actual shape of each individual particle inside the reactor. As a consequence, the fluid flow field is determined by the particles it has to pass by. In the following, we call this approach particle-resolved CFD simulations of fixed-bed reactors to emphasis the resolution of the actual bed shape. The advantage is the absence of transport correlations, which are necessary in pseudohomogeneous or © XXXX American Chemical Society

heterogeneous models to describe the transport of momentum, heat, and mass. The particle-resolved approach is in line with multiscale modeling, or first-principles approach,4 which pursues to describe entirely the system by theory of the actual phenomena. In reality, randomly packed particles inside a tube touch their neighboring particles or the tube wall. Between spherical particles, contact points occur. Contrarily, nonspherical particles stay in contact with other particles or the wall with contact points, lines, or areas, see Figure 1. These contacts have been a great challenge in modeling particle-resolved fixed-bed reactors with the finite-volume method. Close to these regions, the computational cells, with which the calculation region is discretized, show bad quality (i.e., low aspect ratio, large skewness, etc.). These ”bad” cells can lead to a divergence of the numerical solution. Only few authors stated to apply three-dimensional CFD in particle-resolved fixed beds without editing the actual particle shape. Magnico used a fine structured mesh for direct numerical simulations of low flow through a bed of spheres, which he claims ”avoids contact point problems”.5 However, the interstitial flow between 326 and 620 particles was discretized only by 9.6−17.7 million cells. According to the author, mesh independence was not achieved.5 Other authors admitted that particle contacts lead to difficulties of numerical convergence. As a consequence, several strategies have been established over the last years to modify these contacts. The majority of particle-resolved CFD simulations were carried out in beds of spherical particles. Hence, contact-point modifications were investigated. Received: Revised: Accepted: Published: A

September 17, 2016 December 6, 2016 December 13, 2016 December 13, 2016 DOI: 10.1021/acs.iecr.6b03596 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

between experimental results and predictions of local porosity and velocity fields, as well as pressure drop.27 This method was applied in subsequent works studying heat transfer,28 as well as catalytic fixed-bed reactors.29−31 Only few researchers have compared different contact point modifications and their effect toward pressure drop, velocity, and temperature fields. Dixon et al. investigated with CFD two touching spheres regarding drag coefficient and heat flow for the four contact-point modifications for a wide range of flow rates (500 ≤ Rep ≤ 10,000).6 They conclude that global methods result in erroneous drag coefficients. Local methods give much better results. The bridges method is to be preferred for heat transfer problems with cylinder radii of 0.05 ≤ r/dp ≤ 0.1. Moreover, Bu et al. compared the four modifications in different forms of unitcell configurations, i.e., simple cubic, body-centered cubic and face-centered cubic. Flow and heat transfer characteristics have been investigated for flow rates up to Rep = 3000. The authors suggest to use the bridges method with cylinder radii of 0.08 ≤ r/ dp ≤ 0.1.32 The results are in line with ref 6. Spherical particles are only one shape among others for packed beds in industrial applications. More commonly, nonspherical particles are used (i.e., cylinders, one-hole or multihole cylinders, etc.). In contrast to spheres, nonspherical particles show contact points, contact lines, and contact areas with each other and with confining walls (see Figure 1, top right for cylinders). The different contact types can occur in varying extent in a packed bed. Moreover, their effect toward the local transport phenomena can be variable depending, among others, on the individual position, flow regime, and thermal properties. Figure 1 illustrates the four contact-point modifications, which are enlarged for clarity, applied to cylindrical particles. In this study, contact modifications, based on the local caps method and the bridges method, are investigated critically in a bed of cylinders. The influence of the geometrical modifications is evaluated with respect to pressure drop, local velocity fields, heat transfer, and local temperature fields. This study is oriented toward the experiments of Bey and Eigenberger,33 who investigated local velocity in a packed beds of cylinders. Finally, recommendations are given which contact-modification approach should be chosen to model adequately packed beds of cylindrical particles.

Figure 1. Modes of contacts in packed beds (top) and modification strategies (bottom).

Until now, four distinct contact modifications have been applied to avoid low quality cells in the proximity of contact points between spherical particles and between spherical particles and walls: gaps, overlaps, bridges, and caps.6 These four methods can be further classified into global methods (gaps and overlaps) and local methods (bridges and caps). In the first work describing a global contact-point modification, the size of each of the 44 spheres was shrunk by 1%.7 This contact-point modification was later called the gaps method. Until now, this method was used by many authors,7−19 where the shrinking value can range from 0.5%20 to 2%.21 The opposite of the global gaps method is enlarging each particle by a certain amount. Guardo et al. expanded each of the 44 spheres by up to 1%. The spheres were arranged in a structured manner with N = 3.9. Flows with Rep up to 912 were studied.22,23 This global expansion method is called the overlaps method. Also other authors applied this method.23,24 However, the global methods manipulate the entire bed structure in a drastic way. Concerning overall porosity, global contact-point modification leads to an error of approximately 4% for the 99% gaps and 101% overlaps size.6 However, pressure drop is also affected. A rule of thumb was elaborated that 1% change in porosity leads to a 3% deviation of pressure drop.6 Under theses circumstances, local modifications have been considered by other researchers. The so-called bridges method was elaborated by Ookawara et al.25 The authors connected neighboring spherical particles by placing small cylinders in the contact area. The authors claimed that this local modification did not influence significantly the macroporous flow properties over the studied range of flow rate (1 ≤ Rep ≤ 1000).25,26 For heat transfer problems, the thermal conductivity of the bridges can be set independently from the solid and the gas phase. For spherical particles, Dixon et al. suggested that the thermal conductivity of the bridges depends on the distance from the particle taking the Smoluchowski effect into account.6 The counterpart of the bridges method is the so-called caps method. Eppinger et al. were the first authors presenting this kind of approach. They used an elaborated workflow, where the geometry is manipulated dependent on surface proximity. Surfaces close to the touching point are flattened, resulting in a gap filled with fluid cells with a specified width. The authors achieved reasonable agreements

2. MODELING PARTICLE-RESOLVED FIXED-BED REACTORS 2.1. Generation of Randomly Packed Beds. In this study, the randomly packed bed of cylinders is generated with a workflow based on discrete-element method (DEM) simulations. It is described in detail in our previous work.30 The DEM simulation mimics the filling process of a packed bed by including interparticle contact forces into the equations of motion. It is carried out with the CAD software tool STAR-CCM+ 9.06 by CD-adapco.34 In the first step of the workflow, the tube is filled with DEM particles (i.e., spheres or cylinders). Because DEM is formulated for spheres, cylindrical particles are approximated by so-called composite particles. The desired shape of the cylindrical particle is approximated by a user-defined amount of spheres which retain their position with each other. The arrangement of the amount of particles is computed automatically by STARCCM+. However, the more DEM particles, the higher the computational time for simulations. In this study, the cylindrical shape is approximated by 100 DEM particles. This number represents a trade-off between close approximation of the shape and computational time of the DEM simulation. An inadequate B

DOI: 10.1021/acs.iecr.6b03596 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


modifications can be carried out, and the resulting geometry is meshed. Finally, CFD simulations are performed with desired physical models and boundary conditions.30 2.2. Contact Modifications. In this study, local modifications of particle−particle and particle−wall contacts are applied. In Table 1, different contacts of cylindrical particles with distances between centroids and angles between face normals are listed. This table shows that unlike spherical particles, the contacts are located in a certain range of distances and angles. Hence, tolerances have to be introduced defining the type of contact. The choice of these tolerances represents a trade-off in terms of the amount of geometrical modifications. These tolerances are listed in Table 1. The position data and face normals are exported from the DEM simulations. The distance L between centroids of particles and angles between face normals α and angles between a centroid and a face normal β are subsequently calculated by an Octave code, which decides which contact type is present. The workflow of the entire contact modification is shown in Figure 3. In the literature, the different contact scenarios might have different names or tolerances (e.g., see refs 37,38). However, all possible contact cases have been considered.

particle-shape representation can lead to erroneous packing structures, as a recent study showed.35 Figure 2 shows the actual

Figure 2. (A) Actual CAD particle shape, (B) approximated with 100 DEM spheres.

CAD particle shape (A), as well as a composite particle consisting of 100 DEM spheres (B). In principle, the composite DEM approach can be applied to any pellet shape. On the contrary, a novel approach based on the graphics software Blender36 has been used to generate nonspherical packings.24 When the DEM particles are settled, the local position and orientation of each particle is extracted. With this information, a CAD model of the packed bed is built. Afterward, the surfaces are assigned, contact

Table 1. Contact Types of Cylindrical Particles with Distances between Centroids and Angles between Face Normals

C


Article


united at overlaps, which results in several particle clusters, or agglomerations. Second, the caps method can be applied at overlaps. The united method has the disadvantage that heat can be transferred from one particle to another without a resistance. Furthermore, under certain circumstances, there are thin regions to be meshed, which provokes bad-quality cells. The caps method avoids these difficult regions by introducing defined gaps between particles. Then, heat from one particle to another is transferred only through the gas phase, because there is no direct contact between the particles. As a consequence, the bridges method is applied for contact types described in Table 1 in combination with the united method for overlaps. 2.3. Governing Equations. The set of governing equations consists of conservation of total mass, conservation of momentum, and conservation of energy in terms of specific enthalpy. The detailed description can be found elsewhere.39 With the help of the particle Reynolds number Rep, the flow in fixed-bed reactors can be characterized: Rep =

vind p (1)

ν

with vin as the superficial velocity and ν as the kinematic viscosity. The Reynolds number is based on the diameter of a sphere of Vp

equal specific surface: d p = 6· A . For Rep > 300, the flow behaves p

highly unsteady, chaotic, and qualitatively resembling turbulent.40 Turbulent flow was modeled with the Realizable k − ε model, developed by the research in ref 41, combined with a Two-Layer all-y+ Wall Treatment driven by shear. For a complete formulation, see the manual of STAR-CCM+.34

3. EXPERIMENTAL AND NUMERICAL SETUP 3.1. Experimental Setup from Bey and Eigenberger (1997). The numerical study in this work is geared to the experiments from Bey and Eigenberger.33 The authors measured velocity profiles below fixed beds of different particle shapes (i.e., spheres, cylinders, and one-hole cylinders). The packed beds were characterized by their low tube-to particle-diameter ratio, which varied from 3.3 ≤ N ≤ 11 by using different particle sizes. The working fluid was air at ambient pressure with superficial velocities covering 0.5 ≤ vin ≤ 1.5 m/s. The tube had an inner diameter of 25 mm. The bed height is not mentioned in the manuscript. The particles were placed on a monolith of 3.5 mm in length with a channel width of 1 mm and a total porosity of 85%. An anemometer was placed 5 mm below this monolith to measure the local fluid velocity. The position of the sensor was varied in terms of radial and annular position, see Figure 6. Velocity measurements were averaged after repacking the bed to obtain a statistical representative flow profile. For more information on the experimental setup, see ref 33. The authors assume that the monolith has no influence on the velocity profile. However, strong recirculation zones can be noticed below a packing, which was observed for example by CFD simulations.27,29 Still, the setup itself is interesting because low N beds were investigated experimentally with particle dimensions relevant for industrial applications. 3.2. Numerical Setup. A packed bed of cylinders with dp = lp = 12 mm was studied in line with the experiments from 33. The bed of cylinders was created with the workflow described in section 2.1. The resulting tube-to-particle-diameter-ratio is N = 0.05/0.012 = 4.17. The automatically generated fixed bed with all dimensions is illustrated in Figure 7A. There is no bed height or

Figure 3. Workflow of the contact modifications.

Besides the caps method, which flattens touching particles locally, the bridges method is implemented. In Figure 4, this method is shown schematically for the different contact types stated in Table 1. The dimensions of the bridges are enlarged for visualization reasons. The geometrical modifications are all based on inserting extra geometries and applied Boolean operations (i.e., subtract and intersect). First, a larger bridge (blue in Figure 4) is inserted between the particles or between particles and the confining wall. In the next step, this larger bridge is subtracted from the two particles or in the case of wall-side and wall-end contact from the particle. An enlarged bridge geometry is needed only due to the Boolean operation. In step three, extensions are created (light magenta in Figure 4). For the area contact: end− end and the line contact: side−side, an extension of each particle is created, either on the side or at the end. For the other contact types, only one particle is extended. In the final step the extensions are intersected. This leads to the actual bridge shape (magenta in Figure 4). The shape of a bridge is a cylinder at area contacts, and a cuboid at line contacts. This procedure is carried out in STAR-CCM+ by a Java script. Because the composite DEM particles consist of multiple spheres, sharp edges of cylinders cannot be represented. Consequently, overlaps occur where DEM cylinders touch each other with edges. In Figure 5, a schematic of an overlap between two cylindrical particles is shown. For visualization purposes, only the total amount of DEM particles is very low, and their dimensions are enlarged. The resulting overlap regions can be treated generally by two methods. First, the particles can be D


Article


Figure 4. Schematic of the bridges-method workflow for the different contact types.

turned out that a value of 0.015 results in a porosity that is close to the value predicted by the equation above. It has to be mentioned that this value is unrealistically low. However, this value lumps effects like additional compression due to vibrations during filling and operation, as well as additional forces due to fluid flow through the bed. Because the arrangement of the bed does not change in the CFD simulations, these aspects are accounted for in the DEM simulation. Sliney and DellaCorte studied friction and wear of ceramic/ceramic and ceramic/metal combinations over a wide range of temperature. The authors found typical friction coefficients for monolithic alumina in the range of 0.5−0.7.44 However, composite ceramics are characterized by lower values (e.g., mullite showed values between 0.35−0.5). Three different superficial velocities were chosen in line with the experiments (i.e., 0.25, 0.5, and 1 m/s). For the evaluation of pressure drop and heat transfer, the properties of air were set constant to ρf = 1.18 kg/m3, μf = 1.855 × 10−5 Pa·s, cp,f = 1003. 62 J/kg·K, kf = 2.603 × 10−2 W/m·K. This results in particle

Figure 5. Schematic of overlap modification for cylindrical particles.

porosity given in the description of the experiments.33 Because the axial porosity profile is uniform after approximately 3 particle diameters,42 a bed height of ≈10·dp was chosen. The amount of particles in this bed was approximated with the overall porosity ε̅ of a bed of cylinders described by the following equation:43 ε ̅ = 0.36 + 0.1/N + 0.7/N 2

(2)

To match the predicted porosity, the packing generation was tuned by adjusting the static friction coefficient of the DEM particles. This method was suggested earlier in the literature.25 It E


Article


to be smaller than the solid thermal conductivity. Two different values were tested (i.e., kbridge = 0.07 and 0.5 W/m·K). The lower value is close to the gas-phase value. The latter represents 10% of the solid value. The composite-DEM method results in certain overlaps, which can also occur between particles and the confining wall. Hence, in all cases of contact-area modification, a tube with 0.99·D was subtracted from all of the particles touching the tube wall. This ensures a well-defined gap between particles and confining wall. In Table 2 a detail is given of touching cylinders and a wall contact visualizing the different modifications of contacts. As can be seen, the locally modified regions are small in comparison with an overall modification of the particles. Bridges are highlighted in magenta. These bridges are treated as an additional region type in the CFD simulation, which enables the specification of thermal properties of the bridges. The location of contacts is shown in magenta in Figure 7B. Only four cylinders touch the wall with their plane face. In these cases, area contacts were inserted applying the bridges method. More often, line contacts between the confining wall and the shell surface of cylinders can be found. In addition, inside the bed, there exists a number of point contacts and line contacts, respectively. Polyhedral meshes were used with 2−3 prism layers. The thickness of the prism layers was approximated with the thickness of the boundary layer δBL in the stagnation point of a sphere with dp depending on Rep:46

Figure 6. Schematic of experimental setup of velocity field measurement. Adapted with permission from Bey and Eigenberger.33 Copyright 1997 Elsevier.

δ BL = 1.13·Rep−0.5 dp

Reynolds numbers of 191, 382, and 763, respectively. Conjugate heat transfer inside the packed bed was investigated numerically, although there is no corresponding experiment. The wall temperature was set constant at Tw = 373.15 K, and the inlet temperature was set to Tf,in = 300 K. The particles were modeled with the following constant properties: ρp = 1300 kg/m3, cp,p = 1000 Pa·s, kp = 5 W/m·K. The contacts were modified with different approaches described already in section 2.2, that is, the caps method and the bridges method, as well as the united method only for contact points. Thermal conductivity of bridges was kept constant because varying values, as described by ref 45, are too complex to calculate for cylindrical particles. However, the value was chosen

(3)

In Figure 8A, a detail is given of the mesh close to two connecting cylinders treated with the bridges method. On the left-hand side, a line contact can be seen. On the right-hand side the cylinder and the wall are connected with an area-contact bridge. In Figure 8B, the same location is shown for the caps method. At contact areas, the particles are slightly flattened. The different colors represent the different calculation regions (i.e., fluid region in blue, particle region in gray, and bridge region in magenta). The total amount of computational cells varies with contact-area modification (see Table 3). Much finer meshes are needed for

Figure 7. (A) Automatically generated packed bed of cylinders. (B) Contacts highlighted in magenta. (C) View at z = H/2. F


Article

Industrial & Engineering Chemistry Research Table 2. Details of Touching Cylinders and Wall Contact for Different Contact-Area Modifications

Table 3. Porosity and Pressure Drop for Packed Beds of Different Particle Shapes and Different Contact Modifications Δp/ΔH for Rep = particle shape

no. of cells [·10 ]

ε ̅ [-]

1.65

0.415

−1.9

103

331

1110

1.64

0.415

−1.9

103

332

1112

1.66 3.66

0.416 0.413

−1.7 −2.4

102 108

328 383

1075 1309

3.67

0.413

−2.4

117

405

1317

6

united w/o wall contact united w/wall contact caps bridges w/o wall contact bridges w/wall contact

Figure 8. (A) Detail of mesh with bridges method. (B) Same location with caps method.

applying the bridges method because more complex geometries and thinner regions occur.

dev. to eq 2 [%]

191

382

[Pa/m] [Pa/m]

763 [Pa/m]

near-wall region, where the bridges method with wall contact shows the lowest porosity. It becomes clear that these local modification approaches have a low influence on the averaged porosity profile and hence on the averaged local velocity (see Figure 9B). Deviations occur only at the first maximum in the near-wall region. A first impression of the flow field inside the fixed bed gives Figure 10A. Velocity vectors are animated on a plane cut through the bed to get a better idea of changing directions of the flow. The gas enters the bed from the top. Stagnation zones can be identified, as well as channeling close to the wall. Moreover, regions with accelerated flow and recirculation zones occur inside

4. RESULTS 4.1. Porosity, Pressure Drop, and Velocity Profiles. The global porosity of the bed is a first indication of the influence of geometrical modifications. As can be seen in Table 3, the low static friction coefficient leads to a porosity, which is close to the predicted value by eq 2 (≈ − 2%). Furthermore, the difference between the individual contact modifications toward global porosity is low (i.e., < 1%). In Figure 9A, the radial porosity profiles are compared between different contact modifications. The values are obtained by averaging the void fraction over height and over circumference. Differences occur only in the G


Article


Figure 9. (A) Radial porosity profile and (B) velocity of packed-bed of cylinders averaged over height and circumference. Rep = 382.

experimental pressure drop of low N beds and predictions from the Ergun equation in another study.15 Particle orientation and void fraction played a significant role. Integral values like pressure drop, among others, are a first indication for validating the different modification strategies. However, in certain regions, the local flow field can show large deviations between the different approaches. Hence, a detailed analysis of local velocity and temperature is presented for two different positions (i.e., z = H/2 and z = H) and two different Reynolds numbers (Rep = 191 and 763), respectively. A view of the bed at z = H/2 is given in Figure 7C. The position z = H is just below the packed bed. Consequently, there is no solid phase at this position. In Figure 12a and Figure 13a, the velocity in the axial direction is shown. The figure on top shows values after half of the bed height. The figure at the bottom shows values below the bed. All values are averaged in the circumferential direction. In addition, radial porosity is given only for two modification approaches (i.e., caps and bridges method with wall contacts). As can be seen in this figure and also in Figure 7C, the modified region is very small. Only in the center of the bed does radial porosity differ significantly. However, the local velocities vary by up to a factor of 4 for the different modifications at the position z = H/2. In the figure, the velocity in z direction is shown. The deviation is due to the change in geometry at individual positions, which results in a change of velocity direction. A different picture shows the section below the bed (i.e., at z = H). The flow field is mainly dominated by the position of the last layer of particles, which is identical for all of the investigated cases. For the lowflow situation, there exists a stagnation zone at the center of the bed. The axial velocity is close to zero for (R − r)/dp > 1.6. On the contrary, for the high Re situation a back-flow region occurs in the center. This is indicated by negative axial velocities. Except for the near wall region, there are no large deviations between the different contact modifications below the bed. This is also due to the low number of particles at the bottom of the bed. Hence, the influence of contact modifications is low. 4.2. Heat Transfer. An impression of the temperature field inside the fixed bed gives Figure 10B for Rep = 382. The interactions between local flow and local heat transfer can be recognized. Channeling occurs close to the wall, which has a strong effect on the transferred heat from the wall into the bed. The wall Nusselt number Nuw = hwdp/kf is suitable for comparing near-wall heat-transfer characteristics. It can be governed from the wall heat-transfer coefficient hw extracted from the CFD simulations, which is defined as

the bed. The largest velocity on this section through the bed is 8− 9 times larger than the superficial velocity. In Figure 10C, streamlines indicate the flow direction inside the bed, which is from left to right. Below the bed recirculation zones are visible. A comparison of the velocity below the bed between CFD and experiments from ref 33 is shown in Figure 11A. The CFD results are shown only for the caps method because the difference between the methods is very small at this distance below the bed. They are in reasonable agreement over the radius in the range of 0 ≤ (R − r)/dp ≤ 1.3. Close to the center of the tube, the velocities of the CFD simulation differ from the experimental values. The CFD simulation predicts a recirculation zone indicated by negative velocity components in the z direction. In the experiments the particles are placed on a monolith, which first influences the flow direction below the bed, and second prevents recirculation zones. Furthermore, the monolith tends to depress the peaks in velocity below the bed. In addition, the anemometer cannot distinguish between the direction of the flow. Moreover, the velocity below the bed is dominated by the position of the last row of particles. Bey and Eigenberger repacked the bed and repeated the measurement a certain number of times. The CFD simulation was carried out only for one bed geometry. Still, the agreement between CFD and experiments is satisfactory. The averaged velocity profiles for the different modifications do not show large deviations. However, the pressure drop differs by approximately 15% between the lowest value (i.e., the caps method) and the highest value (i.e., the bridges method with cylinders inserted at particle−wall contacts). As can be seen in Table 3, wall cylinders have only little influence on the entire pressure drop. It is more influenced by the type of contact modification. Two groups can be identified in Figure 11B. For three different Reynolds numbers, a parity plot is shown for pressure drop of CFD simulations against predictions from the Ergun equation:47 ⎞ ρ ·vin2 (1 − ε ̅ ) ⎛ Δp (1 − ε ̅ ) ⎜ ⎟⎟ + = f 150 1.75 ⎜ H dp Rep ε̅3 ⎝ ⎠

(4)

with ̅ε as the mean bed porosity and the Reynolds number Rep based on the particle diameter. The bridges method with and without wall contacts shows high pressure drops. The caps method and the united method with and without wall contacts show similar lower pressure drops. Over the investigated velocity range the two groups differ by approximately 10%. For the lowest and moderate flow rate, CFD pressure drop is in a range of ±15% in comparison with the Ergun equation. For the high Reynolds number, the CFD values are lower than the values from the equation. Moreover, large deviations were observed between

hw = H

Q̇ w A w ·ΔTlog

(5) DOI: 10.1021/acs.iecr.6b03596 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Figure 10. (A) Velocity contour and (B) temperature on a plane cut through the bed of cylinders. (C) Streamlines inside the bed. Contacts were modified with the caps method. Rep = 382.

where Q̇ w is the heat transfer rate from the wall in [W], Aw is the surface area of transfer, and ΔTlog is the logarithmic temperature difference defined as ΔTlog =

In Figure 14, the wall Nusselt numbers of the CFD simulations are plotted against Nuw from the correlation of Martin and Nilles:48 Nu w = Nu w,0 + 0.19Pr1/3Rep3/4

(Tw − Tf )in − (Tw − Tf )out ln

(Tw − Tf )in (Tw − Tf )out

(7)

where the Prandtl number is defined as Pr = cpμ/k. Dixon49 recommends to calculate Nuw,0 with

(6)

with Tw being the wall temperature and Tf being the fluid temperature, which is averaged over the cross section at the inlet and outlet, respectively. The heat transfer rate is given in the CFD simulations.

0 ⎛ 5 ⎞⎛ k ⎞ Nu w,0 = ⎜1.3 + ⎟⎜ r ⎟ ⎝ N ⎠⎝ k f ⎠

I

(8) DOI: 10.1021/acs.iecr.6b03596 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Figure 11. (A) Comparison between experiments and CFD simulations of local velocity below fixed bed of cylinders (Rep = 382). (B) Parity plot of pressure drop of CFD simulations and prediction from Ergun equation.

Figure 12. (a) Circumferentially averaged velocity and (b) temperature over radius for different contact region modifications. (A) at z = H/2 and (B) at z = H for Re = 191.

with N = dp/D, and k0r is the thermal conductivity of the bed, which can be calculated by the Zehner-Schlünder eq (see, e.g., ref 50). Three different Reynolds numbers are evaluated. In the figure, intervals of ±10% of the predicted values are represented by dotted lines. Five out of seven contact modifications are located in this interval. Contrarily, the united method with wall contacts and the bridges method with wall contacts and kbridge = 0.5 W/m· K overestimate the wall Nusselt number by approximately 50% and 30%, respectively. The heat transferred in the bed with bridges depends significantly on the effective thermal con-

ductivity of the bridges. Using a very low thermal conductivity of the bridges (i.e., kbridge = 0.07 W/m·K), Nuw is even lower than with the caps method. This is due to convective heat transfer inside the gaps by applying the caps method. This phenomenon is getting more dominant with higher Re. Inside the bridges, only conduction is assumed. However, the largest effect toward transferred heat in the near-wall region is connected with wall contacts. By applying wall contacts with the united method, Nuw is increased by approximately 50%, in comparison with gaps between the wall and particles. The heat transfer is that high because a large portion of the particles are connected directly J


Article


Figure 13. (a) Circumferentially averaged velocity and (b) temperature over radius for different contact region modifications. (A) at z = H/2 and (B) at z = H for Re = 763.

modeling heat transfer in fixed-bed reactors. The bridges method can be extended by cylinders connecting the particles and the wall. The Nuw model, as it is described by 48, shows some weakness for low flow rates and cylindrical particles (cf. the discussion in ref 49). The predicted heat transfer is underestimated by a certain extent, especially for low and moderate flow rates. Furthermore, the model is disadvantageous for low N beds in general. As a consequence, the absolute values of the deviation between Nuw from CFD and from the Martin and Nilles model for cylinders can be misleading. Still, comparing the different contact modifications with each other gives insights into the modeling of low N fixed-bed reactors with particle-resolved CFD simulations. Besides near-wall regions, radial heat transfer in the core of the bed is also affected by the different contact modifications. In Figure 12b and Figure 13b normalized temperatures circumferentially averaged are shown for Rep= 191 and 763. For both flow rates, the united method with wall contacts shows the highest radial temperatures. The bridges method with the lower kbridge and the caps method exhibit similar radial temperatures. As expected, the higher the flow rate, the lower the influence of heat transfer by conduction in the gas phase. Consequently, for high Reynolds numbers, the influence of contact modification is getting less important for radial heat transfer. Convective heat transfer is dominant.

Figure 14. Parity plot of wall Nusselt number from CFD simulations of different contact-area modifications over Nuw from the correlation of ref 48.

with the wall. That means there is no thermal resistance between wall and neighboring particles. Such a situation in reality is only imaginable for sintered particles, for example. Thus, the united method in combination with direct wall contacts is not suited for K


Article


5. CONCLUSIONS Two different local contact modifications were applied and analyzed for a packed bed of cylinders, that is, the caps method and the bridges method for contact lines and contact areas, and the united method, as well as the caps method for overlaps resulting from the composite DEM particles. Furthermore, particle−wall contacts have been considered separately. The investigated bed is characterized by a low tube-to-particlediameter ratio (N = 4.17). Averaged porosity and velocity profiles show low dependencies of the different contact modifications. However, pressure drop gives some idea of the influence of a changing flow pattern. Whereas the caps method shows the lowest pressure drop, the bridges method in combination with cylinders at the wall−particle contacts shows the highest values. This trend can be explained in more detail with local velocity fields. Especially, the bridges method provokes flow separation, which results in higher pressure drop. The results for the heat transfer study show that the higher Rep, the more important the modeling of particle−wall contacts. Contact modifications inside the bed are becoming less important because convective heat transfer is dominant. Over the entire investigated range of Rep the caps method shows encouraging results in comparison with predictions from heat-transfer correlations. This model is in most cases the easiest method to implement and the most stable in terms of calculation. All CFD simulations are compared with the Nuw model from Martin and Nilles, which itself shows disadvantages for cylindrical particles. Still, the differences between the different contact modifications are apparent. A recommendation can be given for modeling transport of momentum and heat in particle-resolved fixed-bed reactors of cylindrical particles. On the one hand, the bridges method with cylinders between wall−particle contacts shows good results for pressure drop and heat transfer prediction. However, the choice of the thermal conductivity of these bridges needs further clarification. On the other hand, the easy-toimplement caps method should be considered for preliminary calculations due to its low time consumption, numerical stability, and straightforward parameter selection. Because most of the nonspherical particles in industrial applications are based on a cylindrical shape, the results and recommendations of this study can be used as a basis for forthcoming investigations. Still, there is a lack of detailed experiments of heat and mass transfer in low N packed beds with which CFD can be validated in an adequate way. Especially, more complex particle shapes (i.e., multihole cylinders) and shaped surfaces are of interest. Furthermore, the influence of contact modifications toward local kinetics should be studied.

■

Additional thanks go to Thomas Eppinger and Nico Jurtz from CD-adapco for fruitful discussions.

■

■

NOMENCLATURE Ap = particle surface area [m2] Aw = wall surface area [m2] cp = specific heat capacity [J/kg·K] dp = particle diameter [m] D = tube diameter [m] hw = wall heat-transfer coefficient [W/m2·K] H = bed height [m] k = thermal conductivity [W/m·K] lp = particle height [m] L = distance between centroids of particles [m] N = tube-to-particle-diameter ratio [-] Nuw = wall Nusselt number Nuw = hwdp/kf [-] Pr = Prandtl number Pr = μcp/kf [-] Q̇ w = heat transfer rate [W] q̇w = specific heat flux [W/m2] r = radial coordinate [m] Rep = particle Reynolds number Rep = vindp/ν [-] T = temperature [K] Vp = volume of particle [m3] vin = superficial velocity [m/s] α = angle between face normals [°] β = angle between centroids [°] δBL = boundary layer thickness [m] ΔP = pressure drop [Pa] ΔTlog = logarithmic temperature difference [K] ε = porosity [-] ε = mean bed porosity [-] ν̅ = kinematic viscosity [m2/s] ρ = density [kg/m3] REFERENCES

(1) Eigenberger, G. Handbook of Heterogeneous Catalysis; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, 2008; Chapter 10.1 Catalytic Fixed-Bed Reactors; pp 2075−2106. (2) Dixon, A. G. Heat Transfer in Fixed Beds at Very Low (

Contact Modifications for CFD Simulations of Fixed-Bed Reactors

Recommend Documents