1. Introduction
Fixed bed reactors are tubes filled with catalytic packing material, on which a chemical reaction may occur. Fixed bed reactors are preferred for their simplified technology and operation, and the design of such beds has been based on effective medium models coupled with correlation-based transport parameters. Low tube-to-particle diameter ratio (N) beds are typically selected for highly exothermic or endothermic reactions, such as in partial oxidations or hydrocarbon reforming. The use of larger particles within a tube of relatively small diameter results in a reduced pressure drop across the bed and, thus, a more economically-favorable operation. However, this results in a higher void fraction near the wall of the bed, which induces significant flow and transport effects that are difficult to account for in conventional packed bed models. In this near wall area, boundary layers develop in which heat and mass transfer are dominated by molecular transport, and a high resistance to transport is present, causing strong changes near the wall.
Dispersion is a phenomenon in which a solute in a porous medium spreads as a result of the combined effects of molecular diffusion and convection in the interstices or the space available for fluid flow between particles. Typical applications of dispersion are, for example, in quantifying contaminant transport in soils or in modeling solute transport in fixed beds. The study of dispersion in packed bed reactors is fundamental to their design in that dispersive phenomena characterize reactant and product transport within the bed [
1,
2,
3,
4]. Dispersion may be modeled analogously to Fickian diffusion by replacing the diffusion coefficients with dispersion coefficients as seen in the following partial differential equation for species concentration:
In Equation (1),
DL is the axial or longitudinal dispersion coefficient and
DT is the radial or transverse dispersion coefficient, related to flow in the freestream and cross-stream directions, respectively. When dealing with dispersion, it is often useful to define dimensionless Péclet numbers, representing the ratio of advective to dispersive transport:
and:
where
vi is the interstitial velocity in the bed and
dp is the particle diameter.
The present study is focused on the transverse dispersion coefficient, rather than on the longitudinal coefficient. The latter is important for such applications as chromatography and the understanding of residence time distribution in fixed beds. Our main motivation for studying dispersion, however, is the analogy between convective heat and mass transfer, since radial heat transfer is one of the most important parameters in fixed bed modeling. Studies of mass dispersion can give insight into the modeling of fluid-phase convective heat transfer.
Strong radial temperature profiles are common in fixed beds of low
N. In contrast, both random walk theory and early experimental results suggested that radial concentration profiles would be fairly flat, although radial dispersion near the tube wall could be reduced due to flow channeling there, leading to radial concentration profiles. Experimental studies of radial dispersion in fixed beds have been conducted for over 70 years (see, for example, [
5,
6,
7]) and have for the most part been conducted by a centerline pulse injection of a tracer substance, with a small number of radial concentrations measured at the bed exit to avoid disturbing the bed packing [
5,
6,
7,
8,
9]. These measurements are made on beds of relatively short length, before the solute has spread over the entire bed radius, to obtain non-uniform radial concentration profiles. The consequence is that for most studies, the wall effects are deliberately excluded, so as to obtain infinite-medium
DT values.
Opinions have been divided on whether there is a wall effect on fixed bed radial dispersion. An early study by Fahien and Smith [
6] used a centerline point source, but allowed the dispersing solute to reach the tube walls. They found that
DT was not constant across the tube radius and that the decrease in
DT near the tube wall depended on the tube-to-particle diameter ratio,
N. The authors correlated their results for radial average transverse dispersion coefficient by:
Wall effects on transverse dispersion were also found in the presence of thermal gradients by Schertz and Bischoff [
10] who attributed the effect to changes in fluid viscosity with temperature. Gunn and Pryce [
8], on the other hand, found no wall effect, although their tracer did not reach the tube walls. Han et al. [
9] also concluded that the values of
DT did not depend on bed location, as they were not time dependent in their transient study. Foumeny et al. [
11] from tracer experiments with internal bed measurements at different bed heights suggested a correlation with a term that explicitly accounted for wall effects and remarked upon the unreliability of data taken from the bed exit, as the flow patterns change rapidly there, and upon the high correlation between the longitudinal and transverse dispersion coefficients in standard models.
In order to investigate the wall effects on dispersion and obtain wall transport coefficients that could be used in heat transfer studies by analogy, several research groups performed experiments in which the tube wall was coated with a dissolving substance [
12,
13,
14,
15,
16,
17,
18,
19]. A wall mass transfer coefficient was introduced, either through a modified boundary condition on Equation (1) [
12,
15,
16,
17,
19] or by the theory for limiting current measurements in an electrochemical system [
13,
14,
18]. A detailed understanding of the transport processes near the tube wall could not be obtained by these methods, due to the small number of radial concentration measurements that could be made.
In an experimental and modeling study, Coelho and Guedes de Carvalho [
20] studied transverse dispersion in granular beds. They used a two-region plug-flow model in which a near-wall region of thickness
δ may be thought of as a laminar sub-layer in which molecular diffusion is the dominant mechanism. Outside of this layer, the transport in the bed is due solely to dispersion. They performed experiments with a soluble wall for different bed lengths and obtained values of the overall mass transfer coefficient
kav to which they were able to fit the two parameters
DT and
δ. Using these parameters, they were able to predict overall mass transfer coefficients for other values of
Re and
L with reasonable agreement. An interesting result of their work was that for a long enough bed, the
δ-region saturates with the soluble species and a one-parameter model suffices. The criterion for a one-parameter model to be valid is:
where
is the bed voidage and
Re is based on the interstitial velocity.
In following experimental studies, Guedes de Carvalho and Delgado [
21,
22] studied lateral dispersion in fixed beds of spheres by including in the bed a soluble cylinder of benzoic acid aligned with the flow of water entering the tube. They neglected any near-wall effects on dispersion, stating that the section of packing contacting the soluble slab actually indents onto this soluble section, removing the near-wall, high-voidage region that is typical in a fixed bed. Solute concentration in the outlet section of the bed was measured as a means to determine the rate of dissolution of the cylinder. Diffusion was treated as occurring in one dimension, given that the bed followed the criterion in Equation (5). The authors noted the strong dependence of
DT/
Dm on
Sc in the range 140 <
Sc < 500 and showed that dispersion behavior was generally independent of particle size. In later work by Delgado [
4], both longitudinal and transverse dispersion in porous media were studied over a wide range of the Schmidt and Péclet numbers. The author supplied a set of equations each for the longitudinal and transverse dispersion coefficients, subdivided by different regimes of dispersion based on the molecular Péclet number: (1) diffusion regime (
Pem < 0.1); (2) predominant diffusional regime (0.1 <
Pem < 4); (3) predominantly mechanical dispersion (4 <
Pem and
Re < 10); (4) pure mechanical dispersion (10 <
Re,
Pem < 10
6); and (5) dispersion beyond the validity of Darcy’s law (
Pem > 10
6). While this paper presents important correlative equations for dispersion, the experimental data were largely from beds above the low
N range (2 <
N < 8). The difficulties of obtaining detailed experimental measurements in complicated systems such as fixed beds has suggested the use of computational tools as a complementary method of obtaining insight into transport properties such as dispersion.
Computational fluid dynamics (CFD) is a numerical method-based computing approach for simulating fluid flow, mass and heat transport and reaction kinetics, among other physical phenomena. In general, the use of CFD requires converting the geometry of interest, here a fixed bed of spheres, into a number of small control volumes, collectively called the computational mesh, or grid. After supplying the appropriate boundary conditions related to flow and species transport and an initial, estimated solution for the system of interest, a complete numerical solution can then be obtained by iterative, convergence-guided, numerical techniques. With the recent introduction of high performance computing capabilities and the continued growth of such computing resources, CFD can now be used as an important tool for the design and simulation of fixed beds. The use of CFD allows the elucidation of certain flow and heat characteristics that cannot be obtained from experimental methods, such as the velocity distribution around and between particles or the temperature at any point within the bed.
An early use of CFD to simulate flow and dispersion in computer-generated random sphere packs was made by Schnitzlein [
23]. The motivation for the study was his observation that radial dispersion should be the sum of molecular diffusion and fluid mechanical phenomena (eddy diffusion in local voids, which act as mixing cells, branching in the packing and channeling due to the placement of particles and the wall effect) and should not be driven by a concentration gradient. This point of view was similar to one suggested by Gunn [
1] who described axial and radial dispersion in terms of probability theory, accounting for dispersion in the fast stream (convective-dominated) and the slow stream near the tube wall (diffusion-dominated). Gunn estimated the probability of a particle existing in the diffusion boundary layer or moving into the fast stream area of the bed. He stated that diffusion in the bed had little effect on convective radial dispersion and introduced a fluid-mechanical Péclet number based on convective-dispersion alone, which in the limit of a high Reynolds number tended to 11.
In Schnitzlein’s work, a two-dimensional map of the void fraction was obtained from the computer-generated packing and used in the porous medium 2D Navier–Stokes equations to obtain a two-dimensional velocity field. The velocity field was then used in a 2D mass balance to generate a simulated dispersion experiment. This was then fit by a standard dispersion model to estimate
DT, which gave a limiting fluid-mechanical Péclet number of 28, which was probably due to low variation in local voidage, leading to underestimation of radial displacement of fluid. A second study by Schnitzlein [
24] used the sphere pack to obtain a 3D network model with inertial terms, but no diffusion or wall friction, which gave a velocity distribution. Using this in the 2D tracer model for
Re > 100 gave
PeT ≅ 12.
Magnico [
25] simulated beds of
N = 5.96 and 7.8 consisting of 326 and 620 spheres, respectively, with a focus on near-wall effects on radial dispersion. He used particle tracing without molecular diffusion and found that particles in a boundary layer at the tube wall of thickness
dp/4 moved longitudinally and tangentially, but not radially. He concluded that exchange with this layer was mainly diffusive. At higher flows,
Re = 200, tracer particles were not easily transported through the boundary layer along the wall, signifying a molecular diffusion-dominated mechanism of mass transfer in this region.
Lagrangian particle tracking along with lattice Boltzmann methods (LBM) were used by Freund et al. [
26] for a tube with
N = 9.3 and
L/
dp = 27. Asymptotic behavior was not reached as the tube was too short. Their computed values of
DT/
Dm were in good agreement with literature experimental results [
8,
9] for 0.1 <
Pem < 10
4 with a tortuosity of
τ = 1.45. Soleymani et al. [
27] simulated dispersion in laminar flow for a packing of 3000 spheres for various larger
N values to get
DT in the post-Darcy region. Results were slightly lower than in [
9]. Jafari et al. [
28] used packings of >3000 particles with
N ≅ 60, with solute continuously injected at the inlet on the centerline, but they estimated
DL only.
Augier et al. [
29] simulated liquid flow in a bed of 440 particles in the form of a square volume taken from the center of a larger bed, to eliminate wall effects. They simulated the radial dispersion experiment of Han et al. [
9] in which a radial step profile is gradually flattened out, for
Re < 100. Their results for
PeT were larger than literature correlations [
3,
4,
21,
22]. This was attributed to the change in void fraction that occurred when the authors reduced their particle sizes by 2% to avoid meshing problems near the contact points. A correction factor was applied that improved the agreement.
In a recent study by Jourak et al. [
30], radial dispersion coefficients were derived by supplying 3D concentration profiles to a 2D effective porous medium model and fitting the coefficients to the data. Their simulation experiments used a bed of regular and randomly packed particles, with a fixed concentration boundary condition at the tube wall, to mimic earlier tracer injection experiments. Results were presented for laminar flow with 0.1 <
Re < 100. The researchers recommended using beds of large width and length, noting that radial dispersion coefficients showed some length dependency. That is, the radial dispersion coefficient was found to decrease as the length of the bed increased. They also noted that a 2D effective medium approach typically predicts dispersion coefficients higher in magnitude than the 3D counterpart. This is due to higher intercellular fluid motion in the lateral direction and, therefore, a higher radial dispersion coefficient.
Dispersion in packed beds of low N is an area of research that demands attention. Previous studies of radial dispersion typically neglect the presence of the tube wall, and these studies are generally conducted in beds of large N. The development of the diffusion-dominated boundary layer induces significant flow effects in the bed, and research that includes this additional dispersion mechanism is important. In a review of previous studies, liquid phase flows were typically the fluid of interest in experimental situations. It is of interest to study dispersion in a gas-phase system, over a range of flow rates, covering different flow regimes from steady-state laminar at low Re through transition regimes, up to so-called turbulent flow at Re > 600, approximately. Experimental tracer injections are highly erroneous when deriving dispersion coefficients: concentration data points are selected from few radial positions, and beds of long length are studied, which yield concentration profiles with small gradients. It is therefore important to study dispersion in beds of developing flow also, in which a clear and total picture of the radial concentration profile is included.
The objective of the present work is therefore to evaluate the ability of standard dispersion models to describe developing concentration profiles in fixed beds of low
N where wall effects are important. This is difficult to do experimentally because of the small number of measurements that can be taken, typically at the exit of the bed and not inside it. We shall do this by computational generation of a number of 3D fixed beds of spheres, with more than 1000 particles in each, over a range of 5.04 <
N < 9.3, followed by resolved-particle simulation of flow and diffusion of a solute from a coated wall in each bed, over a range of
Re. These simulated “data” sets will then be used to provide a rigorous test of 2D effective medium models with a single dispersion parameter
DT, two parameters
DT and a wall mass transfer coefficient and, finally, a model derived from mixing length concepts that gives
DT as a function of radial position [
31,
32].
2. Computational Methodology
The finite volume method implemented in ANSYS Fluent v. 16.2 (ANSYS, Canonsburg, PA, USA) was used for the analysis in the present work. In finite volume methods, the computational domain is divided or discretized into a number of small control volumes. The collection of all control volumes, referred to as the mesh or computational grid, represents the geometry being simulated. The governing equations are integrated over the control volumes to produce algebraic equations for the dependent variables of interest (i.e., pressure, velocity, species mass fraction in this case). In conjunction with the necessary boundary conditions, a complete numerical solution can then be obtained by solving the resulting large set of nonlinear algebraic equations. This solution method is repeated iteratively until specified convergence criteria have been met and a converged solution is reached. Detailed descriptions of the governing equations and the numerical methods can be found in the ANSYS Fluent manual [
33], so only a brief description of the main equations and methods is provided here.
2.1. Fluid Mechanics
The steady-state equation for continuity of mass is given by the following partial differential equation:
The source term
Sm is mass added to the continuous phase through user-defined sources or through phase changes. For the simulations carried out in this work, this term was zero.
The steady-state momentum balance for the fluid in three dimensions is given by the Navier–Stokes equation:
In the above equation,
p is the static pressure, and the term
represents the gravitational body force. The last term
is used for external body forces acting on the fluid and in this study was zero.
The stress tensor
is defined as:
where
is the molecular viscosity,
is the identity tensor and the second term on the right gives the effect of volume dilation due to fluid motion.
2.2. Chemical Species Transport
The chemical species mass fractions
Yi are given by the steady-state convection-diffusion equations:
The last two terms on the right-hand side
Ri and
Si denote the rate of generation of species
i by chemical reaction and the rate of addition of species
i through user-defined sources, respectively. Both were zero in this study.
The term
accounts for the species diffusion flux that occurs from a concentration gradient and is as follows:
where
Di,m is the diffusion of species
i in the mixture
m. To obtain the diffusion coefficient, Fickian diffusion was assumed using the dilute approximation method.
2.3. Geometry and Discretization
Eight beds of different
N were generated, each packed with spherical particles of 0.0254 m diameter. The 3D random beds of spheres were produced using a modified soft-sphere collective rearrangement Monte Carlo algorithm [
34]. The original algorithm was used along with a method to arrange the layer of spheres around the tube wall more realistically at the bottom of the tube [
35]. Each bed consisted of three sections. The first was an empty section spanning 0.0762 m (3
dp) from the entrance of the bed, denoted as the “calming section.” The following section was the packed section, whose length is specified in
Table 1. The final section was another empty section, 0.254 m (10
dp) in length, used to mitigate backflow effects near the exit of the bed. The bulk void fractions in
Table 1 were calculated using the center section of the bed, from 5
dp above the inlet to 3
dp below the outlet, to exclude end effects, except for
N = 5.04, where a shorter length was used to avoid packing defects.
The model geometry is shown for three of the beds in
Figure 1. These beds represent the smallest, largest and an intermediate value of
N. Since only fluid transport mechanisms were to be simulated, only the fluid domain needed to be meshed. Tetrahedral cells characterized the fluid domain, and the control volume linear size was 0.00127 m, which corresponded to
dp/20. Boundary or prism layers extended from the particle and tube wall surfaces into the fluid and were 2.54 × 10
−5 m thick. This gave mesh counts of 30.2–40.3 M cells depending on the value of
N. Mesh refinement studies were carried out on a smaller test column with
N = 3.96, using tetrahedral cells of sizes
dp/16.7 (0.00152 m),
dp/25 (0.001016 m) and
dp/33.3 (0.000763 m), which gave mesh counts from 4.2–30.0 M cells. Excellent agreement was observed for void fraction and interstitial velocity radial profiles for all mesh sizes, with average differences on the order of 1%–2%. Further mesh refinement tests were made using three prism layers on the tube surfaces, again with excellent agreement between profiles.
The contact points between spheres and between the spheres and the tube wall give very narrow “fillets” of fluid, which cause problems in meshing. Various methods have been compared to deal with this difficulty [
36], and a local modification using the “caps” contact point approach [
37] was selected as being suitable for flow problems with no heat transfer. In our implementation of this method, a small cylinder of height 5.13 × 10
−4 m was aligned with its axis along the vector connecting the particle centers and centered at the contact point between the particles, then subtracted from each particle and replaced by a small amount of fluid. At the particle-wall contacts, a similar procedure was followed; this cylinder was 6.15 × 10
−4 m in height, slightly larger because of the opposite curvature of the tube wall. The bed void fraction change from this modification to the geometry was less than 0.5%. The modifications are shown schematically in
Figure 2.
2.4. Boundary Conditions
To obtain velocity profiles in the bed, a uniform velocity boundary condition was set normal to the inlet of the tube, for each of the four particle Reynolds numbers used in this study, namely Re = 87, 348, 696 and 870. For these flow rates, direct numerical simulation was used. Air was specified as the fluid in the bed, with constant properties for isothermal conditions (i.e., μ = 1.7894 × 10−5 kg/m/s, ρ = 1.225 kg/m3). A constant pressure condition was specified at the exit of the tube as atmospheric pressure. No-slip boundary conditions were assigned to the tube walls and particle surfaces.
The focus of this work was to model dispersion and more specifically near-wall effects in low-N beds, by simulating methane diffusion from a coated tube wall into flowing air. To avoid unrepresentative concentration profiles caused by developing flow, the first 5dp axial length of the bed was uncoated to allow flow to become established. The coated section was implemented by setting the methane species mass fraction as unity on the tube walls.
2.5. Computational Aspects for Computational Fluid Dynamics (CFD)
The semi-implicit pressure-linked equations (SIMPLE) method for pressure-velocity coupling was chosen [
33]. The least squares cells approach was selected for the gradient method. For spatial discretization, the first-order upwind scheme was chosen for pressure, momentum and species for the first 200–300 solution iterations, to take advantage of the stability properties. After iterations stabilized, the discretization was switched to the more accurate second-order upwind for the remainder of the solution.
Solution convergence was attained by ensuring that column pressure drop remained constant over at least 1000 iterations. For species transport, a monitor was added to a 2-mm isosurface at the end of the packed section to track methane concentration on each iteration step. Once the concentration was no longer monotonically increasing, the solution was deemed converged. A complete solution for a single case required approximately 15,000–20,000 iterations, with solution times averaging 1000 iteration steps per 12 h. All simulations were completed on a Dell R620 PowerEdge Server (Dell, Round Rock, TX, USA) running a Windows server 2008 R2 (Microsoft, Redmond, WA, USA) operating system. The server contained 2 Intel Xeon E5-2680 CPUs (Intel, Santa Clara, CA, USA), each with 8 cores, and 128 GB of Random-access memory (RAM).
To obtain velocity and methane concentration profiles, cylindrical surfaces aligned axially with the bed were generated at various radial positions spanning the column diameter. The surfaces for axial velocity were clipped five particle diameters from the start of the bed packing and three particle diameters before the end of the packing, to prevent any end effects on the velocities. The resulting velocity profile was therefore averaged both axially and circumferentially. The surfaces for the methane concentrations were created similarly, with those surfaces clipped to 2 mm in height, at different axial positions, giving angularly-averaged methane concentration profiles across the diameter of the bed, at specified axial positions. Area-weighted averages were performed to collect both data sets for all runs.
2.6. Dispersion Models
To obtain dispersion coefficients based on the CFD 3D resolved-particle simulation results, 2D effective medium dispersion models based on Equation (1) were used and their mass transfer parameters fitted to the CFD-generated radial profiles of methane concentration. The finite element commercial code COMSOL Multiphysics
® (COMSOL Inc., Stockholm, Sweden) was used. The two-dimensional domain was subdivided as shown in
Figure 3 to reflect the corresponding sections in the 3D CFD model.
The radii of the beds simulated in COMSOL were the same as those in the CFD simulations. In
Figure 3, the first two and the last sections of the bed represent those parts that were assigned a wall boundary condition of zero methane flux or zero methane concentration. These include the empty sections at the beginning and end of the tube. The “flux” section represents the portion of the 3D model that was assigned a non-zero wall mass fraction of methane. In the two-dimensional model, either a methane wall concentration or a non-zero wall methane flux condition is instead assigned, as detailed in the following subsections. The four horizontal internal boundaries were positioned at the z-values at which concentration radial profiles were extracted from the CFD simulations and were used as the dataset for fitting radial dispersion parameters. For each bed, the highest z-location was always three particle diameters below the top of the bed, to avoid exit effects, while sampling at the most developed profile location possible. The other three z-locations were originally approximately equidistant along the rest of the bed, but as results were analyzed, it was sometimes necessary to dispense with the lowest z-location and replace it with a location further down the bed, which corresponded to more developed concentration profiles.
The 2D domain was meshed with 32 layers of quadrilateral prisms on the tube wall surface, to capture the steep gradients there. The first layer thickness was 1/20 of the local element size, and a stretching factor of 1.3 was applied. The rest of the domain was meshed with triangular elements of sizes between 7.81 × 10−4 m and a maximum of 0.012 m. The coarse mesh sizes were used only in the bed center region, where the concentration profiles were relatively flat. Different mesh refinement settings were used within COMSOL® to check for mesh-independence, which led to the given mesh sizes.
Optimization studies were conducted for each of the 2D models, coupled to the transport models within the COMSOL Multiphysics® framework. The Levenberg-Marquardt least-squares optimization method was employed, which is commonly used for parameter fitting. The objective function was the sum of squares of deviations between the radial concentration profiles of methane, extracted from the CFD simulations at the four bed heights, and the corresponding methane radial concentration profiles from the 2D dispersion model under study. Each optimization was conducted for a given N and Re. The axial Péclet number remained constant at Pea = 2 in all runs. The study was complete when a minimization of the least-squares value of methane radial concentration deviations was achieved and an optimality tolerance (1 × 10−6) was met.
2.6.1. Model 1: Constant DT
In this model, we tried to reproduce the developing radial concentration profile with a constant transverse dispersion coefficient and a radially-varying interstitial velocity:
and at the coated tube wall
r =
R:
Here,
and
. The only fitted parameter was
.
2.6.2. Model 2: Constant DT and Wall km
In this model, we try to take account of the steep gradient near the tube wall by introducing a step change at the wall, governed by a wall mass transfer coefficient:
and at the tube wall
r =
R:
Here,
(again equal to 2),
and
. The two fitted parameters were
and
.
2.6.3. Model 3: Radially-Varying DT
In this model, we use a simplified model for the transverse dispersion coefficient that decreases sharply near the tube wall to take account of the strong decrease in
there, while being constant and neglecting any minor variations in the center of the packed bed:
and at the tube wall
r =
R:
Here,
, and the simplified transverse dispersion coefficient is given by:
where:
which is based on a mixing length model [
31] and re-written following Winterberg et al. [
32] with some differences. The three fitted parameters were
,
and
.
The original version of the model [
32] was developed to include thermal conductivity and was therefore based on superficial velocity in both the differential equation and the equations for radial variation in the parameters, which used the bed-center value obtained by numerical solution of the Brinkman–Forchheimer–Darcy (BFD) differential equation. In our implementation of the model, the interstitial velocity is used for consistency with Model 1 and Model 2, and the bed center velocity is replaced by the bed-average interstitial velocity, obtained from the CFD velocity profiles, to avoid having to solve the BFD equation. Parameter
δ gives the thickness of the region next to the tube wall (in
dp) over which the dispersion coefficient rapidly decreases;
nm controls the shape of the decreasing function; while
PeT represents the constant bed-center transverse dispersion.