Analysis of Particle-Resolved CFD Results for Dispersion in Packed Beds

: Dispersion is the spreading of a solute while it is moved by a ﬂowing medium. The study of dispersion in catalytic chemical reactors is fundamental to their design, since dispersion inﬂuences the reactant and product transport within the bed. In this paper, longitudinal and transverse dispersion of an inert tracer in slender packed beds of spheres and spherocylinders is studied using Computational Fluid Dynamics simulations. The focus is on the analysis of dispersion from full ﬁeld data. The purpose is to develop a methodology that can later also be used to characterize dispersion from full ﬁeld experimental data such as MRI measurements. Results obtained by means of particle-resolved CFD simulations are discussed. Spatial distributions and residence times are analyzed and the results are interpreted by comparison to results obtained from 1D and 2D convection-diffusion equations.


Introduction
Hydrodynamic dispersion is the spread of a solute due to flow.In chemical reactors it is very important because it determines the residence time distributions of reactions and products.This in turn influences conversion and selectivity.In this paper, we study dispersion in packed beds.The regime of the dispersion is hydrodynamic in nature, i.e., molecular diffusion is negligible.Hydrodynamic dispersion is caused by advection of the solute by the flow through the pore space formed by the interstitial volume between the spheres.The variation of local velocities within the pores and different path lengths in a tortuous pore network give rise to the spread.
The method that is used in this study is particle-resolved CFD.For systems with many particles, say hundreds or more, the grid count becomes very high.Therefore, such calculations are computationally demanding and have only become feasible in recent decades.Table 1 lists some of the more recent CFD studies that deal with dispersion in packed beds.These mostly deal with packed beds of spherical particles, but we will also consider a packing of non-spherical particles, namely spherocylinders.Icardi et al. [4] Periodic random packing of irregular and polydisperse objects.
CFD with high-order numerical schemes and advanced meshing techniques.
theoretical results for RTDs measured at many axial positions inside the column and thus an evolution of the RTD while the solute is advected through the pore network.Dispersion 'coefficients' are always linked to an interpretation model.Therefore, we use three interpretation models, namely a one-dimensional convection-diffusion model and two versions of a two-dimensional convection-diffusion model.In all cases, dispersion is modelled mathematically as a diffusion term.The difference in the models is, of course, the radial direction.Variations of porosity and axial velocities are significant near the wall region, especially in slender packed beds, such as those considered here, leading to channeling.In the most extended 2D model, radial porosity and velocity profiles, obtained from the CFD results, are taken into account.In the simplified 2D model uniform porosity and velocities are considered.The 2D models will be used to determine transverse dispersion coefficients.The interplay between transverse dispersion and velocity profiles is also discussed.

Packing Generation
Two random packed beds of mono-disperse polypropylene spheres (d = 3 and 5 mm) were generated with a ballistic deposition approach by means of LIGGGHTS, a DEM code developed by Kloss et al. [13].The parameters shown in Table 2 were used to generate the packings.
A random packed bed of mono-disperse nylon spherocylinders (d = 3 mm, h = 9 mm) was also generated by means of Woo, an open-source DEM code developed by Smilauer [14], originally written in C++ and accessible through Python scripts, which list the geometrical and physical features of the random packed bed to be generated.The parameters used in the Python script are shown in Table 2.

CFD for Flow and Passive Transport
The output from these DEM codes is the 3D position of the particles (and also orientation for spherocylinders).These generated particle configurations provide the porous structure for the CFD code.The Navier-Stokes equations for flow were solved together with the scalar transport equation for the tracer species.A total variation diminishing scheme was used for convective terms and central differencing for diffusive and viscous terms.A projection method was applied to enforce the incompressibility condition of the flow field.
A second order spatial discretization was used on a staggered uniform Cartesian grid.At the particle surfaces, a no-slip boundary condition was imposed by means of a ghost-cell immersed boundary method.This means there is no need for a body conforming grid.The immersed boundary method is second order in space and therefore no-slip boundary conditions are imposed accurately.For that reason, we did not use a refinement of the grid near particle boundaries.A value of 20 grid cells per particle diameter was used to obtain accurate results in a reasonable computational time.This resolution was justified by the convergence study presented in Appendix A of Das et al. [11].
The methods, geometry, grid resolution, and so forth, are very similar to those used in the paper of Das et al. [11].We refer to that paper for the full details on the methods and numerical implementation.In that paper one of the cases is flow (and heat transport) in a slender bed with a diameter of 21 mm filled with 3 mm particles at the same resolution as in the current paper.The only difference is that our bed is a bit longer, containing 900 instead of 680 particles.An impression of the geometry of the beds used here can be obtained from the slices shown in Figure 1.The parameters used in the simulations for spheres with diameter d = 3 mm, d = 5 mm and spherocylinders with d = 3 mm and shaft-length h = 9 mm are summarized in Table 3.
The no-flux boundary conditions for the passive scalar on the particle surfaces and cylinder wall were also imposed by an immersed boundary method.The approach used was developed by Das et al. [11], in accordance with the principles detailed in Deen et al. [15] and Deen et al. [16].
From time t = t 0 an inert tracer was continuously injected at a constant concentration and flow rate.At the entrance of the column, water was introduced with a constant velocity over the whole cross-section of the column except for a small hole (with a diameter of 1 mm) at the center of the bottom plate where water was injected at a higher velocity.After 0.2 s, when the velocity field throughout the column was developed, a tracer was added in this central jet with a fixed and uniform concentration.
Depending on the geometry and velocity we define a particle Reynolds number as: where for the spherocylinders an equivalent particle diameter computed from the particle volume is used.For all three beds, one Reynolds number was considered, namely Re p = 50 based on v equal to the background velocity (excluding the volumetric flow of the central jet).For the choice v = v sup (including tracer flow), the average superficial velocity, the values are 57.8,68.8, and 64.0, respectively.Using an average interstitial velocity, v = v sup /ε, gives Reynolds number values of 135.0, 148.7 and 130.2.Experimentally, dispersion in a packed bed is commonly quantified by means of a residence time distribution (RTD).The RTD, E(t), is the probability density that particles leave a reactor at a time t after they entered.It can be measured by applying a pulse injection of a tracer at t = 0 and measuring the concentration at the exit.Experimentally, and also in simulations, the so-called cumulative RTD, F(t), is easier to measure because it is the response to a step-change of concentration at the entrance rather than a pulse.It is the fraction of molecules that remains in a reactor for a duration less than t.The two quantities are related as: As any probability density, E(t) is normalized such that lim t→∞ F(t) = 1.
Here, we introduce a somewhat extended definition, E(z, t), that denotes the residence time distribution at axial position z.It is the normalized cup-averaged flow at axial position z at time t for a concentration pulse applied at z = 0 for time t = 0. Note that inside a reactor there can be back-flow, so there is a subtle difference in definition with the ordinary residence distribution that characterizes only exiting molecules.The extended definition E(z, t) is useful in cases where there is full field information such as simulations or full-field non-invasive monitoring of dispersion, e.g., with MRI.Additionally, for this generalized definition, the cumulative F(z, t) is related to E(z, t) via Equation (2).
The residence time distribution, E(z, t), has a peaked shape.An important characterization of E(z, t), interpreted as a probability density function, is by its moments: where the square brackets are used to indicate expectation values with respect to the probability density, E(z, t).The equalities among integrals can be obtained by using Equation (2).Note that here (1 − F) is used instead of F to facilitate the partial integration and that the last equality is not valid for n = 0. Computation of moments of E(t) from the accumulative RTD was proposed by Yu et al. [12].In our simulations, F(z, t i ) is generated as output (at discrete axial positions z).To compute the moments, we used a finite sum to approximate the integral where t n dF ≈ (0.5 t i + 0.5 t i+1 ) n (F i+1 − F i ).
From the moments one can construct so-called cumulants that are more convenient to characterize a peak, Here, τ is the mean residence time and σ the standard deviation.The shape of the peak of E(z, t) can be further characterized by the skewness, k 3 /σ 3 , and excess kurtosis, k 4 /σ 4 .
Because of our ability to determine the RTDs as a function of the axial position we will here provide some observation on the z-dependence.Some insight can be obtained by considering the Laplace transform of E(z, t), One use of this Laplace transform is as a way to compute cumulants as: One might expect that for large enough segments, ∆z, the dispersion in a segment is independent of that in the dispersion already present.If this is the case, dispersion over two segments can be obtained by feeding the output of the first section as input to the second section.This can be expressed as a convolution: A convolution becomes a product for the Laplace transform such that: If this is true, one finds a linear scaling with the axial position z, Here, the factor ln( ẽ(∆z, s))/∆z is expected to become almost independent of ∆z for ∆z large enough.It then follows that if Equation ( 8) holds, all cumulants, computed via Equation ( 4), individually scale linearly with z.This linear scaling can be tested and provides insight into the correlation of dispersion along the axial direction in the column.If all cumulants scale linearly, this indicates that mixing in different parts of the reactor is uncorrelated.Conversely, if linear scaling is not found it shows that, for this length of reactor, mixing in different parts of the reactor is correlated.
The linear scaling of cumulants also gives that for large z the standard deviation scales as σ ∝ √ z, k 3 /σ 3 ∝ 1/ √ z and k 4 /σ 4 ∝ 1/z.Thus, the skewness and excess kurtosis tend to zero for large z.This is in accordance with the central limit theory that states that a Gaussian distribution is recovered when summing statistically independent stochastic variables.

Interpretation Models for Dispersion 2.4.1. 1D Convection-Diffusion Model
The simplest interpretation model of axial dispersion is the solution of a 1D convectiondiffusion equation, where dispersion is modelled as a diffusion term with axial/longitudinal dispersion coefficient D L .
∂c ∂t where c(z, t) is the concentration at position z at time t.Note that in this equation we did not introduce a porosity, and the velocity, v, should be interpreted as an interstitial velocity.For a 1D convection-diffusion equation on a semi-infinite domain z ∈ [0, ∞), where at z = 0 a no-flux boundary conditions is imposed after the pulse has been applied, the Laplace transform of the residence distribution is: In Appendix A this expression is derived, and also E(z, t) is provided.The Laplace transform is shown here to point out it is, as expected, linearly dependent on z, not only for large z, but for any z.In this formula, we use a Péclet number based on the axial position.This is also referred to as the Bond number and compares the axial position with the length scale D L /v.By applying Equation ( 6) the first four cumulants are found as: The residence time equals distance divided by velocity: τ = z/v.The second cumulant is the variance, or standard deviation squared.It can be expressed as: The standard deviation in the residence time distribution is related to the standard deviation in space, which is √ 2D L τ for diffusive processes, divided by the velocity.The linear dependence of σ 2 on z obtained from the 1D convection-diffusion equation can be used to determine the dispersion coefficient D L experimentally (or from a simulation).In practice, however, there are always entrance and exit effects.Therefore, near the entrance (and exit) region, there will likely be a nonlinear dependence of σ 2 on z.The procedure we use to obtain D L is to, first, identify the region where the variance grows linear with z and fit a straight line there.The axial dispersion coefficient is then obtained by identifying (2D L /v 3 ) with the slope of this line.
For the skewness and excess kurtosis we see that it decays with increasing z as z and κ 4 /σ 4 = 30 Pe −1 z ), such that the residence time distribution peak becomes more and more Gaussian for increasing z.In the other limit, i.e., z → 0 the skewness and kurtosis diverge.These deviations from the limiting Gaussian shape remain significant even for relatively large Pe z .This is illustrated in Figure A1a, which shows the RTD corresponding to a 1D convection-diffusion equation.
Note that, conventionally, the RTD is evaluated at the outlet of the column, while we here also consider intermediate axial positions.In the 1D interpretation model, different choices can be made for the boundary condition at the outlet.The exact dependence of the spreading characterized by σ 2 (at the outlet) on the Péclet number Pe L is dependent on this choice.See, for example, Table 3 in [17].For axial positions away from the exit, however, this influence of the outlet boundary condition decays on a length scale D L /v.In Appendix A, this is further analyzed.

2D Convection-Diffusion Model
A more extended interpretation model is a 2D convection-diffusion equation.Next to the axial (or longitudinal) dispersion coefficient, D L , a transverse dispersion coefficient, D T , is used.According to this model the concentration, c(z, r, t), evolves as: In this equation, ε is the porosity and v the interstitial axial velocity.We consider two models: the first where both are constant and the second where both have a radial profile, i.e., ε(r) and v(r).The equation will be solved numerically, where we use the porosity and velocity profiles extracted from the CFD simulation results.
In the case that v is constant, the 1D and 2D models provide the same RTD.In the case of profiles, v(r) and ε(r), the transverse dispersion influences the longitudinal dispersion because the transverse dispersion makes sure that the profiles are sampled.This is similar to Taylor dispersion in channel flow.The limiting 1D dispersion coefficient can be computed from the velocity and porosity profile as: In this equation ε and v va are the radially averaged porosity and volume averaged interstitial velocity, respectively.The volume averaged interstitial velocity is defined as v va = ε v / ε , so it equals the superficial velocity divided by the average porosity.
The numerical factor α depends on the shape of the dimensionless radial velocity and porosity profiles.Here, the particle diameter d p is used as a length scale to make the radial coordinate dimensionless.In Appendix B, the formulae needed to compute α are provided.When performing a tracer experiment with a point injection, rather then a uniform concentration over the inlet cross-section, interpretation of the RTD is also complicated by the presence of a radial concentration profile.Integration over the radial direction of Equation (13) gives: Here, two different types of averages appear, indicated by subscripts 'va' and 'ca', respectively.The subscript 'ca' is used to denote cup-averaging, which for concentration equals the molar flow rate divided by the volumetric flow rate at an axial position z.So, in fact, when dispersion is modeled as a diffusive term it is actually not a real average.In this case, c ca = v c − D L ∂c/∂z va / v va .For concentration-step inflow by definition of the accumulative RTD, c ca = c ca,in F. Using Equation ( 3), an interesting result is that for the residence time distribution, by integrating time and space, where c va and c ca are the stationary state values of the volume-averaged and cupaveraged concentrations.For a uniform inlet, the stationary-state concentration in the fluid will be uniform in the radial and axial directions and the integrand equals one.For pointinjection there can be a deviation between the stationary volume-averaged concentration, which has a z-dependence, and the cup-averaged one, which is independent of z.
It is also good to note that, when there are no radial profiles for the porosity and velocity, one can integrate the 2D equations, Equation ( 13), over the radial direction and obtain the closed-form 1D Equation (10).It follows in this case that results for axial dispersion deduced from the 1D and 2D convection-dispersion equations are the same.In conclusion, the integrand in Equation ( 16) can only deviate from one when there is both a non-uniform injection and a non-uniform velocity profile.
Next to the interplay between radial and axial dispersion, the 2D interpretation model can also be used to quantify transverse dispersion in a more direct way.By multiplying Equation ( 13) by r 2 and integrating over the radial direction, we obtain the following relation: We can define r 2 c = r 2 c ca / c ca , where c ca is independent of z in steady state.This root-mean-square average of r can be interpreted as an average where the axial molar flux is used instead of the concentration.
In the simplified case where radial profiles are modeled as uniform, then for steady state (with This shows that r 2 c initially grows linearly with slope 4D T /v.At larger z values the concentration at the wall, c(R), will build up and the growth of r 2 c will level off, which is especially important in slender packed beds.Note that for a uniform concentration in a circular cross-section, with diameter d bed , r 2 c = d 2 bed /8.The relation between r 2 c and D T is not fully closed because of the appearance of c(R) in the formula.For example, the value of c(R) depends on D L via back-diffusion.An analysis of the 2D convection-diffusion equation gives that this effect is small when the dimensionless number D L D T /(v 2 d 2 bed ) is small.For mechanical dispersion, one expects D L to be of order v d p with a prefactor of order 1 and D T has a similar scaling but with a smaller prefactor.Therefore, the dimensionless number scales as (d p /d bed ) 2 .
In the case where a radial velocity and/or porosity profile is present, the formulas become more complicated.As a general approach, we will obtain values of r 2 c from the particle-resolved CFD simulations and fit this to numerical results of the 2D convectiondiffusion model, Equation (13).This is performed for cases with and without radial profiles of velocity and porosity.

Post-Processing
The core datasets are porosity, velocity and concentration fields for three packings obtained from particle-resolved CFD simulations.To obtain profiles as function of the axial coordinate, z, slabs of one cell width were evaluated.Volume-averaged quantities were computed by averaging over the computational cells in a slab that have their centers inside the fluid phase (and not the solid phase).For radial profiles, binning was performed based on the distance of the cell's center to the center line of the tube.This provides an average over the angular direction.
For computing RTDs, cup-averaging is needed.For cup-averaged quantities, the molar fluxes in the axial direction were used as weights in the averaging procedure.From the CFD results, the convective part of the molar fluxes were computed as local concentration times and local axial velocity.In principle, there is also a diffusive contribution to the molar flux, which is related to the concentration gradient and molecular diffusivity via Fick's law.Note, however, that we used a low diffusion coefficient (typical for a liquid).Therefore, the diffusive contribution was negligible compared to the convective one, as indicated by the large Péclet number based on the particle size: Pe p,m = vd p /D m = Re p Sc ≈ 50 × 10 3 .Note that the resolution used in the simulations was 20 cells per particle diameter.Therefore, even the numerical cell-based Péclet number was very large.Diffusion observed in the particle-resolved CFD results will be due to numerical false diffusion rather than being caused by the modeled very small molecular diffusion.Therefore, we did not attempt to fit our results, e.g., with literature correlations, using the particle-based molecular Péclet number, Pe p,m .

Porosity, Velocity and Concentration Profiles
Figure 1 shows the porosity profiles for the three packings considered in this study.The profiles in the middle row are angularly averaged.This provides a good impression of both the order in radial an axial directions as induced by the tube walls.The layering is most significant for the 5 mm particle.Due to the low ratio between tube diameter and particle diameter, the ordering induced by the walls persists radially for a significant part of the tube also in the case of 3 mm spheres.The ordering decays more quickly for the spherocylinders than for the spheres.The bottom of the packing as created by DEM corresponds to the left side in the plots.This explains that the particles are layered more on the left than on the right side.The radial profiles shown in the bottom row of graphs were obtained by averaging the radial profiles over the axial direction.This was performed over the 'bulk' region that we quite arbitrarily defined as starting 3 d p from the bottom of the packing up to 3 d p from its top.The averaged porosities in this region for the three cases are: 0.428, 0.463 and 0.492, respectively.From left to right, the results for the 3 particle types are shown: spheres with diameter of 3 mm, 5 mm and spherocylinders with 3 mm and 9 mm shaft.The rows from top to bottom show a lateral cross-section through the center, the angularly averaged porosity profiles and the radial porosity profile laterally averaged over the bulk region.
The stationary axial-velocity profiles for the three cases are shown in Figure 2. The tracer component is introduced in the center with a jet.This local high velocity jet decays very quickly due to interaction with the packing.Comparison of local and averaged velocities show that, on the local scale, there is quite some heterogeneity.This is induced by the heterogeneity in the local packing structure, inducing significant dispersion.The angularly averaged velocity fields show a clear radial profile.There is preferential flow close to the tube walls and a decaying oscillatory as a function of the radial distance from the wall.This radial profile of axial velocities is induced by the radial porosity profiles of Figure 1.In the bottom rows of Figure 2, both the profiles of superficial and interstitial velocities are shown.The fact that there is a significant effect on the interstitial velocity illustrates that in regions of higher porosity there is also a lower flow resistance, giving rise to channeling. Figure 3 shows the concentration profiles at steady state.On the left side, concentration enters the system at a significantly higher velocity than the background velocity.Because the excess fluid velocity decays quickly, what is of importance is the molar flow rate of the tracer.After angularly averaging the concentration profile, an initial peak becomes clear, and broadens when moving downstream.This broadening will be used to characterize the transverse dispersion.3 the dimensionless inlet concentrations at the inlet are 7.1, 3.6 and 4.4, respectively.For the well-mixed state in the packed region, the volume averaged dimensionless concentration equals 1/ε.The middle row depicts angularly volume-averaged concentration profiles, while the bottom row shows angularly cup-averaged concentration profiles.The black regions near the center-lines correspond to angular regions that do not contain fluid, such that a concentration cannot be defined.

Residence Time Distributions and Axial Dispersion
In Figure 4, the cumulative residence time distributions are shown for the three packings in the top, while the first three cumulants computed from these distributions are shown in the bottom.The cumulative RTD, F(z, t), is computed by monitoring the cup-averaged concentration at axial position z, (normalized by the steady state value) as a function of time, i.e., the breakthrough curve.It can be clearly seen that the breakthrough curve widens with increasing z.To quantify the z-dependence of the RTD we computed the first three cumulants of E(z, t) from the results of F(z, t) by means of Equations ( 3) and ( 4).According to the reasoning in Section 2.3 and also the expressions of the 1D convection-diffusion interpretation model, Equation (12), one expects a linear dependence with z for cumulants when entrance effects have decayed.This is observed in the bulk region for the residence time τ(z) and the variance σ 2 (z) but not for the skewness.The results are presented in a dimensionless manner using the particle diameter, d p , and the superficial velocity.The superficial instead of the average interstitial velocity is used because, for incompressible flow as considered here, the superficial velocity is independent of axial position, so it is a less ambiguous quantity independent of local porosities.There is a clear linear regime for the average residence time as a function of z.
The 1D interpretation model provides a strict prediction for the proportionality constant, namely, τ = z/v int = εz/v sup .The slopes in the τ(z)-curves, shown in Figure 4, are therefore expected to be the average porosities.However, for the three cases, the slopes found by the fit shown in Figure 4 are 0.496, 0.484 and 0.530, which is significantly higher than the bulk porosities of 0.428, 0.463 and 0.492.This difference might be partly explained by Equation ( 16), because after the entrance region, the computed volume average concentration is about 5% to 7% higher than the cup-averaged concentration.
For the variance, σ 2 , in the regime where the dependence is linear, the 1D convectiondiffusion interpretation model predicts a slope 2D L z/v 3 = 2D L zε 3 /v 3  sup .The dimensionless slopes provided in Figure 4 for σ 2 according to this expression then correspond to D L ε 3 /(v sup d p ).The values are 0.262, 0.259 and 0.388, respectively, for the three packings.The discrepancy between the bulk porosity of the packing and the 'porosity' obtained from the average residence time gives an uncertainty for the value of ε that should be used in the formula.The fact that it appears to the power of 3 in the formula increases the uncertainty significantly.The ranges of values found are D L /(v d p ) in the interval [1.07, 1.43] for 3 mm spheres, in [1.11, 1.21] for 5 mm spheres and in [1.38, 1.60] for the packing of spherocylinders.Here, the smaller values correspond to the porosities as obtained from the proportionality constant for the residence time and the larger values from the void fraction.Often dispersion coefficients are presented as a Péclet number where particle diameter and interstitial velocities (v = v sup /ε) are used.The range of the Péclet numbers, Pe p,ax = v d p /D L , is in the intervals [0.70, 0.94], [0.83, 0.90] and [0.62, 0.72], for the three cases.This is consistent with the values reported by Delgado [18] for the parameters of the simulations (Re p ≈ 130, Sc = 1000, Pe p,m ≈ 1.3 × 10 5 , based on interstitial velocities).
For the third cumulant, κ 3 , the lines in Figure 4 are not convincing straight lines.This might provide a clue about the fact that there is still some correlation along the axial direction.The slope of a straight line as predicted by the 1D model, Equation (12), would be ε 3 × (D L /v d p ) 2 ≈ 0.15 for the two spherical cases and 0.3 for the spherocylinders.This is of the correct magnitude for the curves in the figure.

Transverse Dispersion
As discussed in Section 2.4.2, a way to deduce a transverse dispersion coefficient is by analyzing r 2 c in steady state as a function of axial position z. Figure 5 shows the results of r 2 c computed from the particle-resolved CFD, for the three particle packings considered in this paper.For all three cases, there is a steep increase in the radial spread with axial position at the entrance.The characteristic axial distance over which this fast spread occurs is about three times the particle diameter, and after this the value of r 2 c is 1.5 to 3 times d 2 p .This initial quick spread is also very clear from the contour plots in Figure 3.It is caused by the impact of the high velocity jet, which contains the tracer, with the first particles it encounters.After this initial region in the bed, the high velocity of the jet is dissipated, and the tracer moves with the velocity field that is established.This is the region where we tried to fit the results of the 2D convection-diffusion model.The data used for the model fitting are indicated by the symbols in Figure 5. Two versions of the 2D model were considered, the first with uniform profiles for velocities and porosities.These fits are indicated as 'no prof.' in Figure 5.The second version is the numerical solution of the 2D convection-diffusion model with radial profiles for porosities and velocities (as shown in Figures 1 and 2).The results of the model that fit with this version are indicated as 'with prof.'.The value for the axial dispersion coefficient, D L , was fixed at the value obtained from the analysis of the residence time distribution.
Two parameters were fitted, namely, D T and a shift of the origin along the axial direction.The purpose of the shift is to match the fast spreading caused by the initial impact of the jet.For the 3 mm case the fitted curves convincingly follow the CFD values, and result in a value of D T = 0.175 v sup d p = 0.075 v d p in the case where the 2D model without radial porosity and velocity profiles is used.The value obtained from the model that included profiles is D T = 0.143 v sup d p ≈ 0.061 v d p .The pre-factors 0.075 and 0.061 are consistent with the value of 1/12 provided by Delgado [18] for the case at high molecular particle-based P'eclèt.The fluctuation in the CFD results is larger in magnitude than the difference in the fitted curves.The fits provide no evidence that the result from the 2D model with profiles is better.For the case of 3 mm spherical particles, there is a clear growth of r 2 c ; however, fluctuations are also significant.Curves fitted from the 2D models give a slightly smaller pre-factor compared to the 3 mm case.However, due to the large level of fluctuations, no conclusion can be drawn from the difference between the 3 mm and 5 mm cases for transverse dispersion.
As illustrated in the third plot of Figure 5, in the case of the spherocylinders the initial spread due to the colliding jet is so large that subsequently no dispersion can be measured.The lines in the graph show the limiting values of r 2 c for the 2D interpretation model with and without radial profiles of porosity and axial velocity.The limiting values are quite different.This is due to the fact that the average is a cup-average.Therefore, even if the concentration is homogeneously spread, the value of r 2 c depends on the velocity profile.
As discussed in Section 2.4.2, there is a coupling between transverse and longitudinal dispersion.The radial velocity profile gives rise to extra longitudinal dispersion because it is sampled in the transverse direction.This is expressed in a formula similar to Taylor dispersion, Equation (14).Here the pre-factor depends on the profiles and can be computed by evaluating Equation (A13).This expression was numerically evaluated (using spline fits of the profiles) and provided values for α of: 8.3 × 10 −3 , 4.2 × 10 −3 and 4.3 × 10 −2 , for the 3 mm, 5 mm spheres and spherocylinders, respectively.For the 3 mm spheres, using the value D T = 0.061 v d p gives by means of Equation ( 14) D Taylor = 0.14 v d p .This is small compared to D L = 1.4 v d p , which corresponds to Pe p,ax = 0.70 in Section 3.2.

Discussion and Conclusions
In this paper, longitudinal and transverse dispersion have been analyzed from data generated by particle-resolved CFD simulations for two slender packed beds of spheres and a packed bed of spherocylinders.For each bed, only one Reynolds number was considered.This is too limited to draw general conclusions.
The datasets were mostly meant to illustrate the proposed analysis method for full field data.The goal of the research is to further develop methods to characterize dispersion from full field data, such as particle-resolved CFD and MRI measurements.The final goal of the research program is to compare particle-resolved simulations and full field measurements, and in this way enhance the insight into dispersion phenomena in packed beds.Here we will also discuss some developments needed to reach this long term goal.
As concerning the analysis method to determine the axial dispersion coefficient from the z-dependence of the variance of σ 2 as computed from F(z, t), this seems quite accurate.Inside the bed there was a clear region where σ 2 grew linearly with z.It should be noted that in the exit region there was a rapid change in the slope of σ 2 (z).This suggests that determining it from the RTD measured at the exit of a packed bed is quite tricky, and obtaining a good value is sensitive to boundary conditions used in the interpretation model.
It should be remarked that there seems to be some inaccuracy in determining the mean residence time.The interstitial velocity measured from it is not fully consistent with the measured porosity.This is partly explained by dead volume due to the use of point injection.To avoid this complication it is recommended, when interested in axial dispersion, to supply the tracer at the entrance uniformly over the cross-section.
Only for the case with the smaller 3 mm spheres, we were able to determine a transverse dispersion coefficient.Supplying a tracer by means of an injection with a high velocity relative to the background gives an initial quick spread due to the collision of the jet with the packing.For slender beds, this then limits the amount of radial space still available for transverse dispersion.It might be preferred to introduce the tracer with a velocity close to the background value (but more concentrated).Clearly in simulations there is no issue in doing this.In experiments there might be issues.
As far as the particle-resolved simulations are concerned, clearly a larger parameter space needs to be investigated.At a fixed Schmidt number, a range of Reynolds numbers need to be investigated.The question of the influence of channeling in slender beds on dispersion is interesting for chemical engineering practices.However, to investigate that systematically, it would be good to have a non-confined reference case, e.g., a random particle-packing with periodic boundaries in the lateral directions, and next vary the tubeto-particle diameter ratio systematically.Once the reference case for spherical particles is well established, packings of non-spherical particles and/or polydisperse mixtures of particles can be studied and compared systematically.
In the current paper, we considered one Schmidt number characteristic for liquids.Experimentally, it is found that, for high Schmidt numbers, the axial dispersion coefficient only depends on particle Reynolds number but not on the Schmidt number.See, for example, Figure 8 in Delgado [18].For Re p ≈ 130 we find Pe p,ax = v d p /D L ≈ 0.7 for spheres, while that experimental graph suggests a value of about 1.4.Our larger dispersion can tentatively be explained by the channeling in the slender beds.
There is a mismatch between experimental data and theoretical predictions that we have not addressed yet.In the present simulations, a very low molecular diffusivity was used.It was taken so low that one might expect that molecular diffusion wsould be dominated by false numerical diffusion in the CFD simulations.We, without further consideration, assumed that the mechanism of dispersion was so-called mechanical dispersion and the diffusion contributions can be neglected.In mechanical dispersion, the tracer is assumed to flow along the streamlines of the velocity field.The spread in residence times is caused by the tortuosity in the random bed.When non-dimensionalizing the longitudinal dispersion with the molecular diffusivity, for mechanical dispersion, D L /D m ∝ Sc f (Re p ), which reduces for Stokes flow to D L /D m ∝ Sc Re p = Pe p,m .Here, D L does not depend on D m .However, at high Péclet numbers, according to theory, "holdup dispersion" is expected to dominate, as was convincingly argued by Koch and Brady [19] for Stokes flow.If there are circulation regions (or porous particles) present in the packed bed, the tracer will diffuse in these stagnant regions.The smaller the diffusion coefficient, the longer the residence time in the stagnant regions.This long tail in the RTD will give rise to high axial dispersion coefficients.In the case of stagnant regions, Koch and Brady [19] predict D L /D m ∝ Pe 2 p,m .Even in the absence of stagnant regions, but with no-slip boundary conditions on the particle surfaces, the theoretical development of Koch and Brady [19] predicts: D L /D m ∝ Pe p,m logPe p,m .This is sometimes referred to as boundary layer dispersion.
For large Péclet numbers, one might expect this dependency to show up in the axial dispersion, but apparently experimentally it does not.This is in need of an explanation.Either there is an issue with the theoretical development, or it is related to the measurements.One hypothesis we have is that in common experiments the phenomenon has a very long timescale.It provides contributions to the tail of the RTD.This tail decays slowly and, even though the amount of substance with these long residence times is possibly small, mathematically it would contribute to σ 2 when integrating E(t) from 0 to ∞.However, experimentally, the tail is not measured accurately because of an unfavorable signal-tonoise ratio and/or too short measurement times.Note that, for very long beds with residence times longer than the hypothesized long timescale, the diffusive phenomenon would contribute to the width of the front (via the central limit theorem) and would be measurable.
CFD simulations can provide extra insight into such a mismatch of theory and experiment, especially because there are fewer sources of noise and better control of parameters.However, there are limiting inaccuracies, e.g., due to finite resolution in simulations too.For high Schmidt numbers, as considered here, very thin boundary layers will be present.Although our simulations presented here are particle-resolved, the grid is rather coarse with a resolution of 20, and does not resolve the mass boundary layers concerned.To investigate the contributions of holdup and boundary layer dispersion, simulations need to be performed that resolve these boundary layers.One possible approach is the adaptive grid refinement as proposed by Panda et al. [20].This should then be combined with a systematic variation of the Schmidt number.The analysis methods based on cumulants (or Laplace transform) of the RTD, E(z, t), as proposed in the current paper should be of use to shed light on the inconsistency between experiments and theory.As our hypothesis is true, the long time-tail contribution will prevent the linear scaling with axial position for moderate values of z.
It can be questioned whether, in the case of significant holdup dispersion, it actually makes sense to capture the long timescale behavior in a dispersion coefficient.The axial dispersion coefficient, D L , is a parameter in an interpretation model, namely the 1D convection-dispersion model.Lumping the long timescale behavior in a parameter limits the applicability of the interpretation model.It might be much preferred to explicitly model the long time dynamics.Further, in many applications of packed beds, particles are porous.Therefore, when lumping all dispersion phenomena in a dispersion coefficient, holdup dispersion will be significant.However, if something happens inside the particle, such as a catalytic reaction, this again will also influence the dispersion of the substances.In such cases one would like to have an interpretation model, say a reactor model in the case of a catalytically active packed bed, that treats the fluid and solid phases separately.Such reasoning advocates the use of interpretation models more tailored to the application that also have tailored parameters.It matters for the value of a model parameter, such as a dispersion coefficient, if a phenomenon (such as holdup dispersion) is lumped into the parameter or is modeled explicitly.Instead of posing general dispersion coefficients, it might make more sense to fit parameters of custom interpretation models using simulation data and experiments.An interesting limit, although not directly relevant for our application, is that of D L large.With the diffusive flux equal to zero at the outlet, the 1D convection-diffusion model in this limit should become equivalent to a fully mixed tank.In this limit E(t) = L/v exp(−v t/L) independent of z and thus E(z, s) → 1/(1 + (s L/v)).

Appendix B. General Expression for Taylor Dispersion
Here we provide a general derivation for the contribution to the effective 1D axial dispersion coefficient due to the interplay of transverse dispersion and a radial velocity profile, i.e., Taylor dispersion.The starting point is Equation (13) where both the superficial velocity, v, and the porosity profile, ε, are assumed to have a radial dependence.When integrating this equation over the radial direction the 1D equation is where the angular brackets denote averaging over the radial direction.To get a closed form 1D equation ε vc needs to be expressed in terms of c.When we consider the limit D T → ∞ the mixing in the radial direction is so fast that the concentration is uniform.Let's call this concentration, c∞ .In this limit ε v c = ε v c∞ .For large D T we will make the decomposition c = c∞ + c, where c is small.When substituting this expression in Equation ( 13) and maintaining the leading order terms (in c) this gives The time derivative was eliminated by substituting its expression from Equation (A9) in the limit D T → ∞.The resulting final equation is solved for c with the constraint ε c = 0.With this expression, The extra term does not contribute because of the constraint ε c = 0.However, if this constraint is not enforced and c is determined up to an integration constant, this integration constant does not contribute (because its proportionality constant is ε v − ε v = 0. Finally, where the numerical constant α follows from the dimensionless problem: Here, we choose d p as dimensionless length scale such that r = d p ζ.The reason to choose this length scale (instead of e.g., D bed ) is that the velocity and porosity profiles deviate from the bulk values near the wall on a characteristic length of several d p .

Figure 1 .
Figure 1.Porosity maps for the packings in the cylindrical tube of 21 mm generated using DEM.From left to right, the results for the 3 particle types are shown: spheres with diameter of 3 mm, 5 mm and spherocylinders with 3 mm and 9 mm shaft.The rows from top to bottom show a lateral cross-section through the center, the angularly averaged porosity profiles and the radial porosity profile laterally averaged over the bulk region.

Figure 2 .
Figure 2. Axial velocity maps for the three packings.The top row shows a cross-sectional map of the axial velocities.The middle graph shows angularly averaged superficial velocity profiles.For both rows, the velocities are made dimensionless through division by the average superficial velocity.The same scale as indicated by the color bar is used for all color plots.The bottom row depicts radial velocity profiles obtained by averaging the superficial velocity profiles in the axial direction in the bulk region.The red lines show the dimensionless radial interstitial velocity profiles.

Figure 3 .
Figure3.Steady state tracer concentration maps for the three packings.For all graphs, the concentrations are normalized using the cup-averaged concentration, which at steady state is independent of the axial position.The top row shows longitudinal cross-sections through the center of the tube.From the inlet concentrations and flow rates provided in Table3the dimensionless inlet concentrations at the inlet are 7.1, 3.6 and 4.4, respectively.For the well-mixed state in the packed region, the volume averaged dimensionless concentration equals 1/ε.The middle row depicts angularly volume-averaged concentration profiles, while the bottom row shows angularly cup-averaged concentration profiles.The black regions near the center-lines correspond to angular regions that do not contain fluid, such that a concentration cannot be defined.

Figure 4 .
Figure 4.The top row shows accumulative RTDs for the three packings at different axial positions.The are obtained by computing the cup-averaged concentration and normalizing it with the steadystate value.The second row shows the average residence time, variance and skewness as computed by Equations (4) and (3).

Figure 5 .
Figure 5.The spread of the tracer in the radial direction is characterized by r 2 c .The three graphs show this quantity as a function of z computed from the steady state results of the particle-resolved simulation for the three cases: 3 mm and 5 mm spheres and the spherocylinders.The curves are fitted to results obtained from a 2D convection-dispersion model.The symbols indicate the part of the CFD data that is used for the fit.

Figure A1 .
Figure A1.The residence time distribution and its cumulants for 1D convection-diffusion. (a) A normalized depiction of the RTD corresponding to 1D convection-diffusion in an semi-infinite domain.The different curves correspond to different axial positions characterized by the longitudinal Péclet number, Pe z = v z/D L .The inset shows the decay of the tails of the RTDs.(b) For 1D convectiondispersion cumulants of E(z, t) on a finite domain with length L cumulants scale linearly with z until the point where the exit boundary condition becomes relevant.

Table 1 .
Summary of previous CFD modeling studies of dispersion in packed beds.

Table 3 .
Parameters used in the IBM-DNS simulations for spheres with d = 3 mm, spheres with d = 5 mm and spherocylinders with d = 3 mm, h = 9 mm.