How does column packing microstructure affect column efficiency in liquid chromatography?

Journal of Chromatography A, 1126 (2006) 58–69

How does column packing microstructure affect column efficiency in liquid chromatography? Mark R. Schure a,∗ , Robert S. Maier b a

Theoretical Separation Science Laboratory, Rohm and Haas Company, 727 Norristown Road, Box 0904, Spring House, PA 19477-0904, United States b Information Technology Laboratory, U.S. Army Engineer Research and Development Center, Vicksburg, MS 39180, United States Available online 27 June 2006

Abstract Full three-dimensional computer simulations of the fluid flow and dispersion characteristics of model nonporous chromatographic packings are reported. Interstitial porosity and packing defects are varied in an attempt to understand the chromatographic consequences of the packing microstructure. The tracer zone dispersion is calculated in the form of plate height as a function of fluid velocity for seven model particle packs where particles are selectively removed from the packs in clusters of varying size and topology. In an attempt to examine the consequences of loose but random packs, the velocities and zone dispersion of seven defect-free packs are simulated over the range 0.36 ≤
≤ 0.50, where
is the interstitial porosity. The results indicate that defect-free loose packings can give good chromatographic efficiency but the efficiency can vary depending on subtle details of the pack. When the defect population increases, the zone dispersion increases accordingly. For a particle pack where 6% of the particles are removed from an
= 0.36 pack, ≈33% of the column efficiency is lost. These results show that it is far more important in column packing to prevent defect sites leading to inhomogeneous packing rather than obtaining the highest density pack with the smallest interstitial void volume. © 2006 Elsevier B.V. All rights reserved. Keywords: Simulation; Fluid mechanics; Zone broadening; Packed column; Microstructure; Column efficiency

1. Introduction Advances in the art and science of high performance liquid chromatography (HPLC) are strongly dependent on the many parameters governing chromatographic efficiency [1,2]. When packing HPLC columns, control of these parameters is often impossible and a heuristic approach is accepted. For example, the practitioner who packs HPLC columns may choose various solvents, additives, and mechanical procedures based on column performance, and not necessarily based upon fundamental chemical and physical considerations. This approach is necessary in part because of the inherent difficulty in probing the column packing microstructure as a function of various process treatments. In addition, the complete mechanical description of the column including fluid velocities and mechanical stresses at all points is lacking. Hence, column preparation is treated phenomenologically and is largely an art rather than a science.

∗

Corresponding author. Tel.: +1 215 641 7854; fax: +1 215 619 1616. E-mail addresses: [email protected] (M.R. Schure), [email protected] (R.S. Maier). 0021-9673/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.chroma.2006.05.066

It has long been known that column packings are structurally nonuniform, with density variations and efficiencies being spatially dependent [3–18]. Various techniques have been utilized to probe these but mostly all are based on measuring local column dispersion with some form of spectrometry or electrochemical probe placed in a specific region of the column. It has been suggested in many studies that axial variations in packing density lead to minor variations in efficiency while radial variations in packing density lead to major losses of potential efficiency [3– 18]. We have recently developed a simulation-based method capable of performing zone broadening studies of HPLC column packings where the spatial positions of individual packing beads are specified [19,20]. In this paper we attempt to look at the effects of packing heterogeneities and interpret packing defects quantitatively. Specifically, we calculate the zone broadening of model packings produced with differing amounts of packing defects. In addition, we vary the topology of the defect and examine the consequence of the defect shape on efficiency. The interstitial porosity is also varied for randomly uniform packings to compare these results to the defect results and determine the loss of efficiency from loose but uniform packings. The particle

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

59

Table 1 Simulation conditions of the defect-free random packs Interstitial porositya

Random pack no. 1 2A 2B 2C 2D 3 4 a b

Nondimensional permeabilityb

Ω

0.3599 [36] 0.3998 [40] 0.3993 [40] 0.3997 [40] 0.3998 [40] 0.4491 [45] 0.4990 [50]

× 10−4

0.01 0.01 0.10 0.50 0.90 0.50 0.50

6.440 9.769 × 10−4 9.694 × 10−4 1.015 × 10−3 1.000 × 10−3 1.640 × 10−3 2.655 × 10−3

Q6 0.1797 0.1736 0.1866 0.1755 0.1706 0.1883 0.1857

Number in brackets is the nominal porosity as a percentage. Nondimensionalized as k/d 2 , where k is the permeability and d is the sphere diameter.

packs synthesized for these studies are not intended as models of real packings in columns. Rather these packs are to be utilized for the purpose of sensitivity analysis regarding the defects and the loss of plates due to defects. We will show that the interstitial porosity has a minor impact on column efficiency when the packing is macroscopically uniform but still random. However, it is shown that the incorporation of small microstructural defects is a major contribution to the loss of efficiency in chromatographic columns. The asymmetry factor is shown to be a possible diagnostic for packs with percolating defects when asymmetry factors somewhat less than 1 are encountered. These simulations are performed with periodic boundary conditions, which eliminate wall effects that are known to be problematic in chromatographic columns. In addition, the particles are monodisperse. These two assumptions are used here to restrict the complexity of the present study. These aspects will be studied in future papers separately so as to isolate the individual contributions to chromatographic efficiency. 2. Simulations 2.1. Particle packing The method used here for assembling particle packs has been described previously [19,21] although we provide more detail below as some of the results are shown to depend critically on the packing algorithm. The packing algorithm places spherical particles with various packing densities in a rectangular region with periodic boundary conditions (PBC’s) [22,23]. The packing porosities of the defect-free uniform random packs are given in Table 1 along with the calculated permeabilities. The most dense pack (
= 0.36) has been used in previous papers [19,20]. Four packs are utilized for the 40% interstitial porosity studies to examine the effect of microstructural packing variations. We

will refer to all of these packs in terms of the pack number or the nominal interstitial porosity percentages of 36%, 40%, 45%, and 50% as given in Table 1. In the case of the four 40% packs we will refer to these as 2A–2D as noted in Table 1. Monodisperse spheres are randomly packed using a hardsphere Monte Carlo algorithm. Spheres of uniform diameter, d, are arranged initially in a dilute cubic array. Each sphere is moved in a random direction and the move is accepted if it does not result in a collision with another sphere. After a series of such Monte Carlo steps, the minimum distance δ between sphere centers is evaluated and the coordinate system is scaled by the factor d/[d + Ω(δ − d)] while maintaining constant sphere diameter, and where 0 < Ω ≤ 1. By compressing the coordinate system, the packing density is increased. The rate of compression is regulated by Ω. For Ω = 1 the packing is compressed as fast as possible without causing particle intersection. For Ω = 0, the pack is not compressed. Thus, larger values of Ω correspond to more rapid compression. These packing parameters are summarized in Table 1. The interstitial porosities of the particle packs with defect sites are given in Table 2. The process used to create the defective packs starts with the
= 0.36 pack and selectively removes particles in clusters or in line defects. This way empty space is added to the pack, so as to simulate an aggregate of particles with regions that do not pack due to various degrees or types of aggregation, particle jamming, and/or bridging. This approach is at least partly motivated by visualizations of the microstructure of packs discussed previously [24]. No attempt is made here to remove particles that are not anchored solidly on other particles, i.e. particles that “jiggle.” For clusters of particles removed, the number of clusters and the number of particles per cluster can be varied. The sites of cluster removal are chosen under program control to be randomly placed within the computational box. For line defects,

Table 2 Simulation conditions of the defective packs Defective pack no.

Interstitial porosity

No. of clusters removed

D1 D2 D3 D4 D5 D6 D7

0.3626 0.3664 0.3787 0.3914 0.3974 0.3662 0.3830

4 10 20 40 40 1 1

a

Cluster is a line defect parallel with the z-axis.

No. of particles per cluster 4 4 6 5 6 39a 144a

Particles removed (%)

Nondimensional permeability, k/d 2

0.4000 1.000 3.000 5.000 6.000 0.9750 3.600

6.603 × 10−4 6.831 × 10−4 7.542 × 10−4 8.688 × 10−4 9.158 × 10−4 8.483 × 10−4 3.704 × 10−3

60

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

Fig. 1. The seven defect packs used in the study labeled D1–D7 where the spheres shown are the vacant defect sites in the pack. The periodic boundary box edges are highlighted.

control of the number of particles per line cluster and the starting and ending coordinates are specified. We have found that line defects more prevalent and/or larger than one line defect with one particle width per line defect give chromatographic results which yield only a few plates near the minimum plate height. This has set an upper limit to these types of percolating defects, i.e. where the defect allows fluid flow to take a quick route around or through the packing. The defect cluster size is chosen to be 4, 5, or 6 particle sites, as shown in Table 2. These sizes were chosen to simulate effects where both particle jamming and aggregation are present. However, since there is virtually no information on defect length scales in real columns, these parameters are chosen in an ad hoc manner. Because the different packs have different characteristic lengths of defects, the time scale for tracer lateral reequilibration is unique to each particle pack. A study with a larger range of defect sizes, specifically with smaller cluster sizes, is forthcoming from our laboratories. Fig. 1 shows the site of particle removal for all seven of the packs with cluster defects. The visualizations here are composed using a ray-tracing technique [25]. Note that the line defect pack D6 appears to lack full nearest-neighbor connectivity. Nonetheless, fluid and tracer solute can still flow easier through this route as judged by the percent of particles removed (0.975%) and the permeability (8.483 × 10−4 ). This is to be compared with the D4 pack where 5.00% of the particles have been removed and the

permeability (8.688 × 10−4 ) is comparable to the D6 pack. The D7 pack, which has a much larger permeability (3.704 × 10−3 ) as compared to the other packs, has a percolating defect through the pack in the axial flow direction. All particle packs used in this paper have been checked to ensure that particles do not overlap. 2.2. Packing characterization The synthesized packed beds have been analyzed for the amount of orientational order present in the pack using the Q6 metric [26–31]. The structural characterization of the packing material is also analyzed by radial distribution function (RDF) analysis [22,23]. The Q6 metric is computed as follows. The distance from the sphere center of interest to the first coordination shell maximum is chosen by analysis of the RDF of the pack, as suggested previously [28–30]. Spheres with distances less than or equal to this value are connected by bonds and the φ and θ values of these bonds are determined in spherical coordinates. These are then used to obtain Qlm (r) = Ylm (θ(r), φ(r))

(1)

for each bond where r is the radial bond vector and Ylm are spherical harmonic functions [27]. The spherical harmonics are

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

where ti is the arrival time of the ith particle, N  the number of particles, and M1 is the first moment and

averaged over all bonds to get Qlm = Qlm (r)

(2)

and the rotationally invariant quantity Ql is calculated as:
Ql =

l 4π  |Qlm |2 2l + 1

61



N 1  M2 =  (ti − M1 )2 N

(5)

i=1

1/2

and (3)

m=−l



N 1  M3 =  (ti − M1 )3 N

(6)

i=1

We only focus on the Q6 metric as it has a number of desirable properties [26–31]. For purely crystalline packs Q6 = 0.35355, 0.51059, 0.48476, and 0.57452 for simple cubic, body-centered cubic, hexagonal close packed, and face-centered cubic geometry, respectively [28]. The RDF, also denoted as g(r) or g2 (r), is calculated as a function of the nondimensional distance r/d, where r is the dimensional distance and d is the packing particle diameter. The RDF g(r) gives the probability of finding a pair of spheres a distance r apart, relative to the probability expected for a completely random distribution at the same density [22]. The method of calculation is reviewed in refs. [22,23]. We use 300 bins to form the histogram inherent in this method. 2.3. Velocity and dispersion calculations All flow fields utilize grids of 256 × 256 × 1024 (x, y, and z, respectively) and the velocities are computed using the Lattice Boltzmann technique [21,32] that effectively utilizes parallel computing [21]. Four thousand particles are used in the packs which are free of defects, and are packed with the same aspect ratio as used in the flow field grids. All velocities are rescaled from the primary calculated flow field which was computed with a Reynolds number of 0.05. PBC’s are utilized for the flow fields so that long columns are mimicked by passing the solute through the pack multiple times in all coordinates and vapor– fluid interfaces are eliminated. Ten thousand random walkers are used for each dispersion calculation. This was determined to be adequate by comparing model simulations of varying random walker number and examining the resulting moments of the arrival time distribution. The diffusion coefficient D of the random walkers is 8 × 10−10 m2 s−1 and the effective column length is 1.0 cm. Although we show nondimensional results throughout this paper we use a dimensional particle diameter d of 10 ␮m for the purpose of establishing length scales and diffusion coefficients. Other details of these types of simulations are provided in our previous publications [19,20] using these techniques. 2.4. Plate height and asymmetry calculations The first, second, and third moments are easily obtained by direct summation of individual particle arrival times: 

N 1  M1 =  ti N i=1

(4)

for the second and third moments, respectively. The peak skew is M3 S= (7) (M2 )1.5 A common method of reporting peak asymmetry, As , is to report the ratio As = B/A, where B = tB − tmax is the time difference between the peak maximum tmax and tB . The time value where the peak amplitude is 10% of the amplitude at tmax is tB noting that tB > tmax . The quantity A is similar except that A = tmax − tA and tA is the time before the peak maximum where the peak amplitude is 10% of the amplitude at tmax . Rather than smoothing the arrival time histogram and estimating tA , tmax , and tB directly from the peak shape, we use the following procedure which is far more immune to noise. This is essential because the reconstructed peak shape has a finite level of noise and its shape and noise level are dependent on the histogram width used for reconstruction and the number of arrival times, N . We utilize the Exponentially Modified Gaussian (EMG) peak model [33–36] given the parameters σ and τ for that model and estimate A and B using a noiseless peak shape calculation, as described below. Since M2 = σ 2 + τ 2 and M3 = 2τ 3 , as given previously [33–36], the moments obtained through Eqs. (4)–(6) can be used to estimate σ and τ by:    M3 2/3 σ = M2 − (8) 2 and τ=



M3 2

1/3 (9)

The resulting σ and τ values are used to construct an artificial EMG peak of 6 million data points with the subsequent direct estimation of tA , tmax , and tB from these data points. For cases where the simulation peaks are fronting, M3 is negative (the skew is also negative) and the reciprocal relationship As = A/B is utilized in this procedure by substituting the absolute value of the third moment into Eqs. (8) and (9). This is done because the EMG model can only produce tailing peaks but symmetry can be applied to allow a fronting peak to be produced by essentially reversing the independent variable axis. This procedure allows noiseless estimation of As at the expense of invoking a peak shape model into the estimation process. Although the validity of any peak shape model must always be in question, the relative magnitudes of the As values are of importance. Establishing a reliable noise-immune procedure for

62

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

extracting this asymmetry is more important in this work than its absolute estimation accuracy. The plate height H and the reduced plate height h are given as: H=

LM2 M12

and

h=

H d

(10)

where L is the column length. The loss in plates due to loose packings and/or defects can be expressed as a percentage and is calculated as follows. The number of plates in a perfectly packed column is denoted as N∞ . By combining the relationship N = L/H with Eq. (10) to give N = L/ hd, ratios can be written as: N h∞ = N∞ h

(11)

which is also expressed as the percentage loss in efficiency with respect to a “wall-less” column as the quantity (1 − h∞ / h) × 100, where h∞ = 1. The quantity h∞ is typically equated to the plate height minimum hmin which is obtained through model curve fits and described below. Wall effects [13] appear to add a serious amount of zone broadening in packed column chromatography and may explain why many commercial columns have minimum h values of 2. We take a simple approach here and equate the additional zone broadening due to wall effects, hwall , as hwall ≈ 1. Hence, another useful quantity is the percentage loss of efficiency for a column when wall effects are included. In this case h∞ ≈ 2 and this percentage loss of efficiency is given as (1 − [2/(1 + h)]) × 100. These two metrics are used below to quantitate the loss in plates compared to an ideal random packed column. 2.5. Models and curve fitting All results from the zone broadening simulations are obtained by tabulating h as a function of the nondimensional average velocity ν (also called the P´eclet number) at various values of ν noting that ν = vd/D, where v is the average fluid velocity. The h versus ν curves can be fit to a number of models. We have described this procedure previously [19] for a number of the models fit to the 36% porosity pack data. For the study in this paper we use the Knox four-parameter model [37,38] which has the functional form: h = aνn +

b + cν ν

(12)

and has been shown previously [19] to fit zone broadening data extremely well. Another useful model is the Giddings fiveparameter model [39]: h=

2λ 2γ + + cν ν 1 + ων−n

(13)

This model was tried but suffers from convergence problems during parameter estimation due to a lack of linear independence in the λ and ω parameters. This was seen previously in our own studies [19] and others’ studies [40,41]. The parameter estima-

tion appears to be more problematic with these defect loaded particle packs than with the 36% pack used without defects. Other models of interest include the five-parameter Knox model [37]: h=

  1 1 −1 b + cν + +  n ν a dν

(14)

where d  here is an additional parameter that describes the diffusion length needed to smooth out flow irregularities [37]. Eq. (14) also suffers from lack of independence during parameter estimation between the a and d  parameters. The reliability of curve fitting with Eq. (14) appears to increase if a is substituted for 1/a and the reciprocal taken after parameter estimation. The parameter estimation by nonlinear-least squares procedures was described previously [19,20]. The curve fitting is useful for examining the parameter dependencies in Eqs. (12) and (14). A number of physical interpretations have been given to these model parameters [37] and one can probe these parameters with the model packs utilized in this paper. The parameter estimation is also useful as an interpolation basis for determining the plate height minimum, hmin , and the velocity at which hmin is found, νmin |hmin , by local region searching. We use h = hmin for comparing efficiencies between packs described in Eq. (11). 3. Results and discussion 3.1. Interstitial porosity study The conditions used for this study are given in Table 1 and the resulting seven h versus ν curves are shown in Figs. 2 and 3. Two elements of these curves stand out. First, with the exception of the 40% pack results, the more porous packs do not add that much additional zone broadening as compared to the 36% pack. This can be quantitatively viewed in Table 3 where the two methods of looking at plate loss described previously are calculated for the curves at their plate height minima; the additional zone broadening results in only a few percent difference. Second, the four 40% packs show a fairly wide range of plate height curves; these are shown in Fig. 3. The only difference in the packs shown in Fig. 3 is in the preparation of the pack; the porosity is the same for all of these, however, these packs were prepared with different rates of compression, Ω. The distance between nearest neighbor particle surfaces in a crystal lattice will scale [42] as φ−1/3 where the volume fraction of particles φ is equal to 1 −
. This scaling shows that the distance dependence between particle surfaces (and particle– particle center distances) and interstitial porosity is a weak function. These distances control the diffusion length of solute in a chromatographic system and the smaller the distance, the less zone broadening occurs due to Taylor or ‘nonequilibirum” dispersion mechanisms [38,39]. More accurate formulas exist for the mean nearest neighbor and mean nearest surface distance in a random high density pack [31]. The mean nearest neighbor distance scales as a complicated function of
and the pack RDF [31]. For the packs used here the

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

Fig. 2. The dimensionless plate height, h, plotted against the dimensionless average velocity ν or P´eclet number for the four interstitial porosity datasets given in Table 1. The 40% pack is that of the 2C pack given in Table 1. The lines are curve fits from Eq. (12).

Fig. 3. The dimensionless plate height, h, plotted against the dimensionless average velocity ν or P´eclet number for the four 40% porosity datasets given in Table 1 labeled as 2A–2D. The lines are curve fits from Eq. (12). Table 3 Plate loss at the velocity of maximum efficiency νmin |hmin

hmin

Loss in platesa (%)

Loss in platesb (%)

Random packs 1 5.109 2A 5.307 2B 5.226 2C 4.628 2D 4.842 3 4.880 4 4.957

0.9803 0.9156 0.9221 1.094 1.019 1.025 1.025

−2.01c −9.22c −8.45c 8.63 1.88 2.43 2.44

−0.995c −4.40c −4.05c 4.51 0.950 1.23 1.24

Defective packs D1 4.933 D2 4.829 D3 4.091 D4 4.113 D5 3.774

1.004 1.055 1.426 1.776 1.977

0.416 5.18 29.9 43.7 49.4

0.208 2.66 17.6 27.9 32.8

Pack no.

a

With respect to a “wall-less” column. Assuming that wall effects give hwall = 1. c Negative numbers represent the percentage gain in plates compared to a column with hmin = 1. b

63

mean distance between 10 ␮m particle surfaces is 1.144 ␮m, 1.231 ␮m, 1.329 ␮m, and 1.418 ␮m for the 36%, 40%, 45%, and 50% packs, respectively. These numbers are computed for a distance cutoff of r/d = 1.1 noting that the 40% result is the average of the four 40% packs. This shows that the distance between particle surfaces is a weak function of
. This partially explains why the plate height only shows a mild increase with the increase in interstitial porosity. The results of Q6 analysis for the defect-free packs are shown in Table 1. These results suggest [28–30] that the packs are in a noncrystalline, nonequilibrium fluid state typical of granular media [43,44]. There appears to be little difference between these packs as judged by Q6 . This metric gives a bulk value for the pack and regions of crystallinity if they exist can be balanced out with other regions of random packing. We are currently extending this method to analyze the density distribution and spatial distribution of order in the sphere pack noting that correlation functions of local bond order parameters have been given previously [26,29]. It will be shown that the RDF is more useful for examining subtleties between the packs below. The variability in the h versus ν plots in Fig. 3 shows that packing variations for random packs can lead to significant departures from ideal chromatographic behavior. To examine how different these four 40% porosity packs are we view the RDF’s for both the range of porosities and for the four 40% packs. The RDF’s of the packs labeled as 1, 2C, 3, and 4 in Table 1 are shown in Fig. 4 and the four 40% pack RDF’s are shown in Fig. 5. When viewed near the first coordination shell in Fig. 4A where g(r/d) ≈ 1, the packs appear in the correct order of diffuseness with the 36% pack having sharper features than the 50% pack and the other two packs filling in the space between these. The 36% pack result for g(r/d) appears to be exactly the same as the pair correlation function, h(r), given previously in the literature for a random dense sphere packing [45] with
= 0.36 noting that h(r) = g(r) − 1. Finite values of g(r/d) are noted to occur for g(r/d) < 1. This is due to the binning process averaging used to calculate g(r/d) and is a numerical artefact due to a small number of total bins. However, when one inspects the RDF given in Fig. 4B, it can be seen that the 40% pack (pack 2C in Table 1) does not follow the order on the right of the first coordination shell peak maximum; rather the 40% pack is too densely packed at the first coordination shell. This would result from a fraction of nearest neighbor particles all touching and having a looser packing at farther distances. This might occur if aggregate clusters of particles allowed for flow around the aggregates. It appears that packing subtleties in the first coordination shell control the zone broadening in a more important way than zone broadening just being a function of the packing density. The differences in RDF’s for the four 40% packs are shown in Fig. 5. The more efficient chromatographic simulations were derived from packs 2A and 2B. These packs exhibit first coordination shells that are more dispersed on average than the other two packs which show a loss in plate height compared to the 36% porosity pack. However, both 2A and 2B have plate height minima less than unity. This could result from small amounts of

64

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

Fig. 4. (A and B) The radial distribution function, g(r/d), as a function of the dimensionless ratio of distance to particle diameter for the 36%, 40% (dataset 2C), 45%, and 50% porosity data.

crystallinity in the packing. It was previously shown [20] that a face-centered cubic sphere pack gives a plate height minima of h ≈ 0.6. As noted years ago by Lubachevsky and Stillinger [46], their packing algorithm can generate denser, more ordered packs when slower growth rate parameters are utilized. For our packing algorithm using small values of the rate of compression, Ω, the packs may also be retaining some of the crystallinity of the initial starting geometry. On the other hand, the packs with Ω = 0.5 and 0.9 show less efficient performance as judged by Fig. 3. These packs are not purely random either as the RDF analysis demonstrates a localized clumping of particles which have interparticle spacing close to the interparticle contact limit. These results suggest that at the constant density of
= 40%, the microstructural details of how the particles are packed can modulate the chromatographic efficiency over a fairly significant range. In Fig. 3, the variation in h takes place over a range of ≈10% between packs 2C and

Fig. 5. The radial distribution function, g(r/d), as a function of the dimensionless ratio of distance to particle diameter for the four 40% porosity packs.

2D. This range is ≈20% for the 2A–2D packs when possible crystalline effects are included and the plates quantitated for “wall-less” efficiency differences. For most columns, 0.36 ≤
≤ 0.42 [47–49]. The limit of
= 0.36 is often referred to as “dense random packing” [30,50,51] and
= 0.40 as “loose random packing” [30,50,51]. Typically most columns are nearer to the “loose random packing” of
= 0.40 as viewed from the bulk value and not from any microstructural consideration. The results given in this paper suggest that even with the less dense packings the loss in plates is not critical. Since the anomaly in g(r/d) for the 40% pack is short-ranged and probably occurs spatially throughout the packing, the loss in plates here is still not severe. We will see that this is not the case once defect sites are introduced and have specific orientations. 3.2. Defect study The h versus ν results for the packs with defects are shown in Fig. 6 and the percentage of plates lost due to defects is given in Table 3. The data from Table 3 are presented graphically in Fig. 7. As can be seen from these data, a small percentage of defects results in almost no appreciable loss of plates. However, as the number of defect sites exceeds 1%, a significant amount of plate loss is observed. We do not plot the data from the D6 and D7 simulations as these packs have less than four plates at the minimum plate height. Once the defect percolates through the pack, the pack is useless for chromatographic separations. Although the loss of plates is shown to be monotonically increasing in Fig. 7, it is obviously not linear. Even though 6% of the packing sites are evacuated for the D5 pack, which originally had a 36% porosity, the porosity after particle removal is ≈40% and the loss in plates is 25–50% depending on the method of loss quantitation. Hence, defects due to packing anomalies are more

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

65

3.3. Estimation of model parameters

Fig. 6. The dimensionless plate height, h, plotted against the dimensionless average velocity ν or P´eclet number for the five defective pack datasets given in Table 2. The lines are curve fits from Eq. (12).

important than nearest neighbor particle–particle distance variations given in the defect-free random packing studies. Therefore, radial inhomogeneities resulting from macroscopic defects whose origins are in aggregation, particle jamming, particle bridging, and other granular mechanical mechanisms should be avoided if possible in column packing. The RDF analysis of the defect packs D1–D7 does not show significant differences with the RDF of the native 36% pack shown in Fig. 4. This is because RDF’s do not have sensitivity towards capturing individual events. RDF’s show average properties and the defects apparently get averaged out as they occur at relatively small concentrations throughout the packs. This is consistent with a study [52] where the RDF of a sphere pack was studied with one-third of the spheres removed. In this case the RDF was scaled in probability but retained the same shape. It was also shown [52] that the coordination number distribution was radically altered by the removal of these spheres.

Fig. 7. The percent loss in plates from Table 3 as a function of the percentage of defects in the pack.

The results of curve fitting the h versus ν simulation data to the four-parameter model (Eq. (12)) are given in Table 4. The results of curve fitting to the five-parameter model (Eq. (14)) are given in Table 5. The curve fits, although not shown, have included data at ν = 150, 200, 250, 300, 350, and 400. We show simulation results where ν ≤ 100 in our figures because this range is far more relevant to analytical separations. For the purpose of interpreting the mechanical meaning of the parameters, we use this extended range to include the convective region where h exhibits concave behavior as a function of ν [37–39]. A number of trends are observed from this data. First, the four-parameter model fits the data slightly better, as viewed from the relative standard deviation of the fits given in these tables. This was known previously [37,38]; however, the fourparameter model does not have a clear physical description of the coefficients [37]. Second, nearly all of the b coefficients are very close to 2, which is expected for nonporous particles with infinitely small tracer molecules suggesting a near-unity tortuosity [38,39]. Chromatographic theory [38,39] gives the dimensional plate height, H, due to axial diffusion as 2D/v. After making H nondimensional by dividing H by d, the reduced plate height h due to axial diffusion is equal to 2/ν and b should be equal to 2. The c values are very small in both Tables 4 and 5. This is expected for retentionless, small molecule tracers and is consistent with a large number of experimental data [37]. The c term is a combination of stationary and mobile phase mass transport resistance terms [38,39,49], noting that the stationary phase term is zero for nonporous, retentionless solutes and the mobile phase term is small compared to unity. These tables have many negative c values which we assume are an artefact of curve fitting with small magnitude c coefficients. We have discovered that the other parameters are not adequately estimated when c is not included in the model. This suggests that c is small but essential to describe the data from these simulation conditions. The a parameter values are shown in Table 4 and appear to directly correlate with the plate height h. As the packs contain more defects, a gets larger. For the defect-free random packs, the a parameter from Eq. (12) is in the range 0.1838 ≤ a ≤ 0.2625 in an inverse correlation to efficiency. It appears that the a parameter in this model simply scales the plate height so that curves with higher hmin values have higher a parameters. Hence, the physical meaning is complex and dependent on all of the factors which lead to higher minimum plate heights. However, this is the same meaning as the a parameter in the simple van Deemter equation [38,39] where eddy dispersion and other effects were proposed to control zone dispersion as a packing-specific and velocity independent parameter. The a parameter values for the five-parameter model are shown in Table 5. These values do not correlate consistently to a structural model nor the hmin values for the various packs. However, lower hmin values appear to have for the most part higher a values in the five-parameter model. It is difficult to assign any specific physical meaning to the a parameter in the five-parameter model based on the datasets here. Giddings’ model of packed bed chromatography incorporated a as

66

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

Table 4 Curve fit coefficients for the Knox four-parameter model Pack no.

a

b

c

n

RSD (%)

Random packs 1 2A 2B 2C 2D 3 4

0.2104 0.1838 0.1895 0.2625 0.2271 0.2265 0.2309

1.986 1.999 2.012 1.963 1.992 1.964 1.949

−9.826 × 10−3 −1.855 × 10−2 −2.863 × 10−2 −1.041 × 10−2 −2.134 × 10−2 −3.999 × 10−3 3.486 × 10−3

0.6837 0.7452 0.7784 0.6571 0.7237 0.6573 0.6115

1.34 1.29 2.94 1.57 2.23 1.68 1.28

Defective packs D1 D2 D3 D4 D5

0.2278 0.2457 0.4690 0.8396 1.023

2.009 1.966 1.944 1.804 1.712

−3.130 × 10−2 −6.493 × 10−3 2.308 × 10−3 2.007 × 10−2 1.910 × 10−2

0.7475 0.6456 0.4946 0.2839 0.2628

1.34 2.11 1.67 1.55 1.52

a velocity-independent microstructural parameter that described flow-based contributions to zone broadening in so-called transchannel, transparticle, interchannel and transcolumn configurations [39]. The essential lumping of all of these characteristics into one parameter may make a unique physical assignment of the a parameter most difficult. The additional problem in assigning physical meaning to a is that the d  parameter also has microstructural meaning; however, d  is related to diffusive exchange of solute across streamlines [37–39]. From the data in Table 5 it appears that d  is correlated closely to hmin and more efficient packs have smaller d  values. The d  values are in a smaller range than a in Table 5; for nondefective packs we see that 0.1621 ≤ d  ≤ 0.2266. However, when viewing the d  values we see that d  scales monotonically with the defect concentration where 0.1878 ≤ d  ≤ 0.7763. The most interesting parameter in Tables 4 and 5 is the exponent, n. This parameter appears to be well-correlated with the number and type of defect in the defective packs; the larger the percentage of defects, the lower is the n value. The functional dependence of the exponent n appears to be similar for the two models. The n values approximately track the random packs given in packs 1–4 with the efficiencies given in Figs. 2 and 3 and Table 3. The higher efficiency results give a higher value of n. It has been suggested [37] that n is related to “the regularity

of packing.” The n parameter correlates well with efficiency; higher efficiency packs have higher n values. That trend is inverse to that of d  . Perhaps a more comprehensive testing of the model parameters is to more precisely vary the defect size from one-particle vacancies up to the level given in Table 2. However, there is a lower limit here where smaller vacancies will not cause appreciably different h versus ν curves. Considering the simplicity of these four- and five-parameter models, it may be most difficult to ever establish that a model of this simplicity can accurately give physical information about a particle pack that contains a complex RDF, is spatially heterogeneous, and may include long length-scale phenomena. 3.4. Pack permeability and fluid velocity distribution The permeability data shown in Tables 1 and 2 are obtained from the Lattice Boltzmann calculation of the velocity field. These permeabilities correlate well with the interstitial porosity for the defect-free packs. The four 40% packs have permeabilities that differ by a few percent. Although this is expected, it is nonetheless interesting that shifting the nearest-neighbor statistics of the pack while keeping
constant can result in differing permeabilities. This type of variation shows that the fluid permeates through the column differently among the four

Table 5 Curve fit coefficients for the Knox five-parameter model Pack no.

a

b

c

d

n

RSD (%)

Random packs 1 2A 2B 2C 2D 3 4

85.94 13.57 36.12 17.45 17.81 45.44 17.15

1.987 1.986 2.000 1.944 1.998 1.965 1.975

−1.119 × 10−2 6.039 × 10−3 −1.257 × 10−3 4.676 × 10−3 3.113 × 10−3 1.792 × 10−3 1.013 × 10−2

0.1982 0.1621 0.1771 0.2266 0.1995 0.2190 0.2116

0.7164 0.7088 0.7079 0.6757 0.7010 0.6563 0.6336

1.33 1.34 2.47 1.83 2.42 1.84 1.38

Defective packs D1 D2 D3 D4 D5

10.31 33.16 68.81 63.74 38.54

2.015 1.983 1.990 2.047 1.978

6.267 × 10−3 −2.466 × 10−3 −1.768 × 10−3 3.136 × 10−3 6.089 × 10−3

0.1878 0.2251 0.4244 0.6084 0.7763

0.7159 0.6761 0.5152 0.4657 0.4222

1.50 2.13 2.53 3.99 3.77

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

Fig. 8. The flow velocity probability density function for the defect-free packs at different interstitial porosities.

40% packs and this indeed results in different zone broadening results. The velocity density function is shown for the defect-free packs in Fig. 8 and for the defect-filled packs in Fig. 9. In both figures the fluid velocity is made nondimensional by dividing the velocity v by the average velocity v. The seven different defectfree packs have slightly different velocity density functions as shown in Fig. 8. There are subtle differences between the four 40% packs as well and this is expected but the differences are small and are most pronounced at low flow velocities. The 45% and 50% packs have larger permeabilities as compared to the 36% and 40% packs. This is also shown in Fig. 8 where the 45% and 50% packs have less probability of slower flow regions, as compared to the smaller
packs. The zones of higher plate height correlate with the curves of narrowest and most peaked velocity density curves. Hence, the 2C curve shows more flow at a smaller velocity than the other packs. The two line defect packs, D6 and D7, show much narrower velocity densities than the other defect-filled packs in Fig. 9. This is due to most of the flow going through the defect as this

Fig. 9. The flow velocity probability density function for the defect packs.

67

Fig. 10. The skew and asymmetry values for datasets D1 and D6 as a function of the nondimensional flow velocity.

path has the lowest resistance to flow. As was the case for the velocity densities shown in Fig. 8, the zones of higher plate height correlate with the curves of narrowest and most peaked velocity density curves in Fig. 9. The flow velocity differences between D1 and D5 are somewhat subtle but the zone broadening results show dramatically different efficiencies. 3.5. Zone asymmetry The asymmetry values for datasets D1 and D6 are given in Fig. 10. We use D1 as an example of a well-packed column with minimal defects and D6 as an example of a defective column with a serious defect. For D1 it is seen that the asymmetry stays relatively constant and is greater than 1, indicative of a small amount of tailing, for the velocities greater than the velocity where the plate height is minimum. This is contrasted with the D6 dataset where the asymmetry shows a peak fronting effect that increases with velocity. The dataset D1 has an interstitial porosity of 0.3625 and D6 has an interstitial porosity of 0.3662 so the comparison is mostly one of defect topology. A well-packed commercial column gives asymmetry values in line with the D1 simulation, typically 1.05–1.12 and it is relatively constant across the velocity range as also given by the D1 simulation. The D6 pack, which models a badly packed column where multipath fluid flow occurs, shows asymmetry values in the range of 0.7 ≤ As ≤ 0.85. This is an extreme case but shows that fronting, under linear chromatographic conditions, is probably indicative of a poorly packed column where some solute “sneaks” through a defect. This suggests that if the asymmetry or skew decreases dramatically with velocity in a real column, a long length-scale packing defect may exist. An interesting consequence of plotting the data in both As and skew is that these descriptors appear to be almost linearly related in this case. The data almost collapse to one curve for the D1 result and one for the D5 result. The point that these two descriptors of zone asymmetry are linearly different, i.e. As ≈ c1 + c2 S over the range of ν does not appear in the chromatographic literature on peak shape. The extent with which this relationship is universal over a wide range of data is unknown.

68

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

3.6. Structural stability and metastability of packs A question of tantamount importance is whether one can learn more about dispersion mechanisms in column chromatography by utilizing particle packs made with a physically correct description of the packing process, including accurate interparticle forces. The particle packs used in this paper were not intended as being accurate models of real packed columns. Rather, these packs are more homogeneous and lack the load-bearing structures, or force chains, formed by particle contacts in a real packing. One difficulty in modeling a real packing process is that the slurry concentrations are sufficiently high that long length scale collective particle motion must be incorporated with fluid mechanics. Attempts to do this with the lattice Boltzmann algorithm are known [53,54]. For a model of slurry packing to be accurate, surface interactions must be accounted for and under the concentrated conditions of the particle slurry, pair-wise interactions may not be accurate. Three-bodied potentials [55] may be needed with some form of a lubrication approximation to recover the interactions because the particle density is high. In addition, we have not attempted to follow a program of mechanical stability regarding the size of the defects [56]. This is especially critical as it appears from the literature on the physics of granular media that defects can be minimized only through a timed process of mechanical vibration [43,44,56]. Structural inhomogeneities commonly found in monolith columns are more prevalent than those found in packed beds, as revealed by microscopy. These may ultimately be the limiting aspect to monolith performance and may be where monolithic technology can be improved. Many of the concepts used for interpreting variation in column efficiency were described over 40 years ago [38] when Giddings discussed transchannel and transparticle flow. These ideas are descriptively simple yet powerful: defects add a faster way for solute to flow through the packed bed and as such provide a zone broadening pathway which lowers efficiency. It is convenient to use simulations to search for correlation between the parameters in the curve-fitting equations and the fundamental attributes of the pack. However, the most parsimonious models will inevitably remain highly approximate characterizations of defects on both short and long length scales. Simulation, when used empirically, allows a deeper quantitative understanding of the particle packing consequences of chromatographic columns and advances in this field can be guided by these results. Nomenclature a, b, c, d  , n parameters in the four- and five-parameter Knox equation a reciprocal of a used in parameter estimation A tmax − tA As peak asymmetry B tB − tmax c1 , c2 coefficients that approximately equate skew and asymmetry

d particle diameter D diffusion coefficient g(r), g2 (r) radial distribution function (RDF) h nondimensional plate height, h = H/d hmin nondimensional plate height minimum hwall wall effects contribution to plate height h∞ reduced plate height in a perfectly packed column h(r) pair-wise correlation function, h(r) = g(r) − 1 H dimensional plate height L length of column M1 , M2 , M3 first, second, and third moments, respectively N number of plates N number of tracers in simulation N∞ number of plates in a perfectly packed column Qlm (r) spherical harmonic functions of the radial neighbor bond vector Qlm (r)averaged Qlm over all bonds r distance between spheres r radial neighbor bond vector S peak skew tA time of the 10% height location of the peak, tA < tmax tB time of the 10% height location of the peak, tB > tmax ti arrival time of each tracer in dispersion simulation tmax time of the peak maximum v, v velocity and average velocity of fluid x, y coordinates perpendicular to axial flow coordinate Ylm (θ(r), φ(r)) spherical harmonic functions z axial flow coordinate Greek letters γ, λ, ω, c, n parameters in the five-parameter Giddings equation δ minimum distance between spheres
interstitial porosity θ(r) spherical coordinate angle with neighbor bond ν reduced velocity νmin |hmin reduced velocity at the plate height minimum σ standard deviation of peak in time units τ exponential time constant for tailing peaks φ volume fraction of particle pack, φ = 1 −
φ(r) spherical coordinate angle with neighbor bond Ω rate of compression in packing algorithm

Acknowledgments Georges Guiochon has made extremely important contributions to the understanding of column inhomogeneities in chromatography. His experimental work on fluid flow in the wall region of packed columns is highly acknowledged and respected throughout the fluid mechanics community. It is a pleasure and honor to acknowledge Georges and his contributions as he turns 75 years young. MRS gratefully acknowledges the National Scientific Foundation under grant CHE-0213387 and the continued support of the Advanced Biosciences division of the Rohm and Haas Company.

M.R. Schure, R.S. Maier / J. Chromatogr. A 1126 (2006) 58–69

References [1] L.R. Snyder, J.J. Kirkland, Introduction to Modern Liquid Chromatography, John Wiley and Sons, New York, 1979. [2] P.R. Brown, R.A. Hartwick (Eds.), High Performance Liquid Chromatography, John Wiley and Sons, New York, 1989. [3] J.H. Knox, J.F. Parcher, Anal. Chem. 41 (1969) 1599. [4] J.H. Knox, G.R. Laird, P.A. Raven, J. Chromatogr. 12 (1976) 129. [5] E. Baur, E.W. Kristensen, R.M. Wightman, Anal. Chem. 60 (1988) 2334. [6] T. Farkas, M.J. Sepaniak, G. Guiochon, J. Chromatogr. A 740 (1996) 169. [7] T. Yun, G. Guiochon, J. Chromatogr. A 760 (1997) 17. [8] T. Farkas, M.J. Sepaniak, G. Guiochon, AIChE J. 43 (1997) 1964. [9] T. Farkas, G. Guiochon, Anal. Chem. 69 (1997) 4592. [10] M.S. Smith, G. Guiochon, J. Chromatogr. A 827 (1998) 241. [11] K. Miyabe, G. Guiochon, J. Chromatogr. A 830 (1999) 263. [12] B.S. Broyles, R.A. Shalliker, G. Guiochon, J. Chromatogr. A 867 (2000) 71. [13] R.A. Shalliker, B.S. Broyles, G. Guiochon, J. Chromatogr. A 888 (2000) 1. [14] R.A. Shalliker, B.S. Broyles, G. Guiochon, Anal. Chem. 72 (2000) 323. [15] S.G. Harding, H. Baumann, J. Chromatogr. A 905 (2001) 19. [16] A.J. Sederman, P. Alexander, L.F. Gladden, Phys. Fluids 117 (2001) 255. [17] U. Tallarek, T.W.J. Scheenen, H. Van As, J. Phys. Chem. B 105 (2001) 8591. [18] J.C. Park, K. Raghavan, S.J. Gibbs, J. Chromatogr. A 945 (2002) 65. [19] M.R. Schure, R.S. Maier, D.M. Kroll, H.T. Davis, Anal. Chem. 74 (2002) 6006. [20] M.R. Schure, R.S. Maier, D.M. Kroll, H.T. Davis, J. Chromatogr. A 1031 (2004) 79. [21] R.S. Maier, D.M. Kroll, Y.E. Kutsovsky, H.T. Davis, R.S. Bernard, Phys. Fluids 10 (1) (1998) 60. [22] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, New York, 1987. [23] D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, San Diego, CA, 1996. [24] J. Dingenen, Analusis Mag. 26 (7) (1998) 18. [25] POV-Ray—The Persistence of Vision Raytracer, www.povray.org/. [26] P.J. Steinhardt, D.R. Nelson, M. Ronchetti, Phys. Rev. B 28 (1983) 784. [27] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes, second ed., Cambridge University Press, 1992.

69

[28] M.D. Rintoul, S. Torquato, J. Chem. Phys. 105 (1996) 9258. [29] T.M. Truskett, S. Torquato, P.G. Debenedetti, Phys. Rev. E 62 (2000) 993. [30] S. Torquato, T.M. Truskett, P.G. Debenedetti, Phys. Rev. Lett. 84 (2000) 2064. [31] S. Torquato, Random Heterogeneous Materials, Springer-Verlag, New York, 2002. [32] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, Oxford, 2001. [33] I.G. McWilliam, H.G. Bolton, Anal. Chem. 41 (1969) 1755. [34] E. Grushka, Anal. Chem. 44 (1972) 1733. [35] J.P. Foley, J.G. Dorsey, J. Chromatogr. Sci. 22 (1984) 40. [36] D. Hanggi, P.W. Carr, Anal. Chem. 57 (1985) 2395. [37] J.H. Knox, J. Chromatogr. A 960 (2002) 7. [38] J.C. Giddings, Unified Separation Science, John Wiley and Sons, New York, 1991. [39] J.C. Giddings, Dynamics of Chromatography, Marcel Dekker, New York, 1965. [40] E. Katz, K.L. Ogan, R.P.W. Scott, J. Chromatogr. 270 (1983) 51. [41] R.T. Kennedy, J.W. Jorgenson, Anal. Chem. 61 (1989) 1128. [42] W.H. Boersma, J. Laven, N.H. Stein, AIChE J. 36 (1990) 321. [43] D. Bideau, A. Hansen (Eds.), Disorder and Granular Media, North-Holland Publishing, Amsterdam, 1993. [44] A. Coniglio, A. Fierro, H.J. Herrmann, M. Nicodemi M. (Eds.), Unifying Concepts in Granular Media and Glasses, Elsevier, Amsterdam, 2004. [45] A.S. Clarke, J.D. Wiley, Phys. Rev. B 35 (1987) 7350. [46] B.D. Lubachevsky, F.H. Stillinger, J. Stat. Phys. 60 (1990) 561. [47] B.J. Stanley, C.R. Foster, G. Guiochon, J. Chromatogr. A 761 (1997) 41. [48] P. DePhillips, A.M. Lenhoff, J. Chromatogr. A 883 (2000) 39. [49] G. Guiochon, S. Golshan-Shirazi, A.M. Katti, Fundamentals of Preparative and Nonlinear Chromatography, Academic Press, New York, 1994. [50] G.D. Scott, Nature 188 (1960) 908. [51] R.M. German, Particle Packing Characteristics, Metal Powder Industries Federation, Princeton, NJ, 1989. [52] J.D. Bernal, S.V. King, Discuss. Faraday Soc. 43 (1967) 60. [53] A.J.C. Ladd, R. Verberg, J. Stat. Phys. 104 (5/6) (2001) 1191. [54] N.-Q. Nguyen, A.J.C. Ladd, Phys. Rev. E 66 (2002) 046708. [55] C. Russ, H.H. von Gr¨unberg, M. Dijkstra, R. van Roij, Phys. Rev. E 66 (2002) 011402. [56] X.Z. An, R.Y. Yang, K.J. Dong, R.P. Zou, A.B. Yu, Phys. Rev. Lett. 95 (2005) 205502.