Simulation of bidisperse colloidal centrifugal sedimentation using a mixture viscosity model Open Access

Gravitational sedimentation is one of the most important solid–liquid separation methods and is widely used in industries, such as water treatment, minerals processing, and pharmaceuticals, due to its simplicity and flexibility.^1–3 The importance of sedimentation on stability is also critical to new applications of colloidal science, such as nanofluids for energy applications.⁴ Monodisperse systems have been studied intensively,^5–8 but there are several aspects of polydisperse settling that remain relatively poorly understood. For example, it is believed⁹ that the hindrance effect in polydisperse suspensions of low total volume fractions is more complex than for monodisperse particles, because the settling velocity of each species is affected to different degrees by the counter flow of the displaced fluid.

Bidisperse sedimentation can be regarded as a model case¹⁰ for understanding interactions in polydisperse systems while maintaining a relatively high degree of dispersion control. Investigations into bidisperse dynamics are also crucial for systems where size segregation and sorting take place. For example, recent work has considered particle control and migration behavior in microchannels of dense suspensions with acoustic handing,¹¹ where flow stratification is significantly influenced by bidisperse sedimentation theories. A suspension balance model (SBM) for viscous Stokes flow has also been recently extended to bidisperse systems, to investigate the effect of particle size and detail structuring within the suspensions.¹² Improving our understanding of bidisperse suspension sedimentation is, therefore, of considerable importance for particle segregation and flows in many industrial processes,^13,14 and a number of studies have been reported on the analysis of the sedimentation in bidisperse^15–17 or polydisperse suspensions.^18–20 Nevertheless, the effects of particle size and concentration ratios are not fully established, while most techniques employed for simulating such systems are generally very computationally intensive.

For well dispersed spherical systems, the Stokes settling velocity of individual particles is well known,²¹ where a hindrance effect will increase with concentration due to the combined effects of viscosity enhancement by the nearby particles, as empirically proposed by Richardson and Zaki.⁵ They derived a relatively simple “hindered settling” equation to account for this effect that gives an empirical power-law function. Later, Bachelor²² introduced a linear settling law for bidisperse and polydisperse sedimentation that was restricted to dilute suspensions, which has been extended to a wider range by authors such as Bachelor and Wen⁸ and by Al-Naafa and Selim.²³ 

Over the past decades, there have been several experimental and modeling approaches to quantifying the critical parameters in bidisperse sedimentation. For example, Thies-Weesie et al.²⁴ examined large diameter ratio silica particles, where the settling velocity of the larger particles was reduced in correlation to the addition of the smaller species. In this case, differences were attributed to surface irregularities on the large spheres with the adsorption of the small particles, which increased the hydrodynamic friction. Non-colloidal suspensions were also studied by Cheung et al.,²⁵ who found that the interparticle interactions of bidisperse systems are non-negligible, especially in condensed suspensions. Hu and Guo²⁶ simulated the interaction between particle clusters and single particles, showing the dynamics of clusters is significantly influenced by particle interactions. Cao et al.¹⁵ conducted centrifugal sedimentation of polymethylmethacrylate (pMMA) colloidal suspensions, where they proposed a new sedimentation coefficient for bidisperse systems to distinguish velocities of settling occurring independently and collectively. The formula describes the relationship between the sedimentation coefficient and porosity, linking the diameter ratio and the mixing ratio.

There are also other complex simulations using approaches such as the Masliyah–Lockett–Bassoon (MLB) model,^19,27 force coupling method (FCM),²⁸ effective-medium model,²⁹ or effective voidage model³⁰ to predict the settling velocities of two interfaces, which require significant computation time. Furthermore, the settling of different shaped particles, in particular ellipses,^31,32 has been simulated to establish the dynamics of systems with wider application. Other non-spherical particles such as rods^33,34 are also drawing attention in recent years and are more oriented toward industrial applications such as pharmaceuticals. Very recently, Oshima et al.³⁵ have simulated the sedimentation under a centrifugal field of bidisperse SiO₂ and CeO₂ nanoparticles, by modifying an advection-diffusion model, as used by Antonopoulou et al.⁷ for monodisperse silica. Particle size effects were incorporated by considering an averaged composite Gaussian function (and so each sized wasn't explicitly simulated). Nevertheless, several important observations could be determined on the nature of nanoparticle separation of different size classes, with good correlation between experiments and modeling, especially for particles ≥100 nm.

To simplify sedimentation simulations, suspension mixture viscosity models³⁶ are one of the most popular methods to predict the influence of concentration on settling, with the advantage of reducing the velocity to zero at the maximum volume fraction. Bidisperse mixture viscosity models have also been developed, such as those by Qi and Tanner,^37,38 although as of yet these have not been incorporated into sedimentation predictions. In their model, a critical factor to calculate the relative viscosity is the maximum volume fraction of bidisperse suspensions, under the principle that the smaller particles can partially fill the voids between the larger fractions. Additionally, particle–particle interactions are usually neglected for convenience,^39–41 where the maximum volume fraction is based on a hard-sphere assumption (which means there is no restriction on the minimum separation distance between two particles). However, the assumption is not suitable for highly charged colloidal particles,⁷ due to strong electrostatic repulsion and the electrical double layer, which can, thus, significantly impact effective volume fractions increasing specific viscosity. To describe colloidal suspensions, Metin⁴² demonstrated the use of an effective particle size taking into account the minimum interparticle distance, to calculate an effective maximum volume fraction. In the work of Antonopoulou et al.,⁷ interparticle interactions were captured through the modeling of Derjaguin–Landau–Verwey–Overbeek (DLVO) theory, where the secondary predicted minimum interparticle distance was used to estimate the effective particle diameter, and such methods may be extended to bidisperse and polydisperse systems.

Thus, there is still a critical need for validated bidisperse colloidal settling models that are computationally inexpensive, which can confirm particle fraction behavior and track both species. There is an important need to investigate the influence of different size and concentration ratios, as well as to incorporate particle–particle interactions. Although complex hydrodynamic direct numerical simulations may reflect every detail in the sedimentation process, they are very computationally time-consuming and it is difficult to track zonal behavior at a full sample scale. Therefore, as an alternative, a one-dimensional (1-D) advection-diffusion bidisperse sedimentation simulation is proposed here, based on a modified mixture viscosity model, to track the concentration and settling rates of disperse colloidal particle fractions. Specifically, a continuum-based approach has been used with numerical implementation using an adaptive finite-difference solver. Particle phase concentrations are tracked to predict the bidisperse velocities of interfaces for dilute (less than 3 vol. %) colloidal particle suspensions with a variable size ratio under both centrifuge and normal gravity conditions. In particular, spherical silica mixtures of 500 and 800 nm (size ratio of 1.6) and 100 and 500 nm (size ratio of 5) are considered. A mathematical bidisperse mixture viscosity model has also been applied, where hindered settling effects are explicitly considered,⁴³ and results are compared to the empirical Richardson–Zaki (R–Z) function. Furthermore, an effective maximum volume fraction⁴² method is also used, where interparticle repulsions are taken into account. Simulations of concentration profiles, fractional interfaces vs time, and averaged settling rates are also compared to experimental data for the same systems (measured by analytical centrifuge) for validation. The model is shown to be a close fit to the data for both size ratio dispersions at three mixing ratios, with differences <10%. Slightly higher levels of variation are detected for the 500:800 nm system, owing to the smaller size ratio and potentially greater resulting effect of uncounted secondary hydrodynamic factors.

The colloidal silicas used were Angstrom Sphere^® silica powder, with nominal quoted mean sizes of 100, 500, and 800 nm (Fiber Optic Center, Inc., USA). The particles were initially characterized with SEM to confirm that their sizes were close to the manufacturer's estimates, as used in the modeling (see the previous paper by the current authors⁴⁴). The particle densities (ρ_p) were measured as 2.20 g/cm³ for 100 nm silica particles and 1.92 g/cm³ for 500 and 800 nm particles, using a Pycnomatic ATC gas pycnometer (Thermo Electron, USA). The density results were close to expectations for silica particles.⁴⁵ Suspensions were prepared by adding silica powder to a 1 × 10⁻⁴ M potassium chloride (KCl) background electrolyte solution, using KCl crystalline powder (Fluka Chemie GmbH, Germany) and ultrapure Milli-Q^TM water, with a resistivity of 18.2 MΩ cm at 298 K (Millipore, USA). The suspensions were placed in an ultrasonic bath (XUBA3, Grant) for 15 min prior to characterization without heating. Then, samples were further dispersed using an ultrasonic probe (Sonic Dismembrator, Fisher Scientific) at 80% amplitude for 5 min. The purpose of the ultrasonic bath and probe was to fully homogenize the suspensions before any measurements were taken. Bidisperse nanoparticle suspensions of 100 and 500 nm particles (labelled 100:500) as well as 500 and 800 nm particles (labelled 500:800) were prepared in various mixing ratios (on a volume basis) as shown in Table I.

Sedimentation of both monodisperses and bidisperse suspensions was monitored using a LUMiSizer^® analytical centrifuge (LUM GmbH, Germany). The LUMiSizer⁴⁶ enhances sedimentation by exposing sample dispersions to a centrifugal acceleration greater than Earth gravity. The centrifuge rotation rate can be altered from 0 to 4000 rpm, and the effective distance (radial) ranges from 105 to 130 mm. For the experiments conducted, a sample height of 21 mm was used (with a corresponding sample volume of 430 μl), while the rotation rate depended on the bidisperse size fraction (see Table I). A light source pulsed near-infrared (865 nm) light through the side of each sample cell at user-specified times. The light intensity was normalized prior to each run. A 25 mm 2048 element CCD-line detected the intensity of transmitted light across the length of the sample, yielding transmission profiles. The transmission profiles were recorded every 10 s. The intensity of transmitted light is related to the particle volume fraction according to the Lambert-Beer law,⁴⁷ with high transmission intensity corresponding to a low local volume fraction. For a complete description of the technique, please refer to a previous publication by the current authors.⁴⁴ 

To verify the sedimentation under Earth gravity conditions, a Turbiscan^® (Formulation, France) Earth gravity optical analyzer was used. Here, 500:800 nm particle dispersions were placed in 20 ml sedimentation cells. As with the LUMiSizer, the sedimentation fronts were analyzed through a pulsed optical source (near infrared, 880 nm) scanned every 10 min automatically. The Turbiscan utilizes detectors at 180° (transmission) and 45° (backscatter) allowing analysis of opaque dispersions. The whole measuring time was more than 48 h due to the slow settling time of the colloidal particles under Earth gravity.

To verify predicted concentration profiles, a LUMiReader X-Ray^® (LUM GmbH, Germany) was also used. This device utilizes the same cells as the LUMiSizer, where x-rays rather than light is used to track the settling fronts. Due to the greater penetration of x-rays, information on the relative concentration of dispersions is also available by measuring the per distance attenuation. This technique cannot, however, directly process samples under centrifuge. Therefore, 500:800 nm mixtures at total volume fractions of 0.03 were first accelerated at 1000 rpm in the LUMiSizer for different time periods, whereupon they were placed within the LUMiReader X-Ray and the x-ray attenuation was profiled for different settling times.

When the suspension is monodisperse, such as in the case when the larger particle fraction has completely settled out and only the smaller particles remain suspended, the viscosity model reverts back to the monodisperse modified Quemada model.⁷ It is recognized that the DLVO model used here becomes less accurate for particles that are smaller than 500 nm in diameter.^58,59 For particles as small as 100 nm, as an example, the calculation of separation distance is several times larger than the diameter of particles in a very diluted salt solution (such as 10⁻⁴ M KCl). Therefore, the standard DLVO secondary separation distance might be several times larger than the particle size, which is not physical for a particle sediment bed under centrifugal compression. However, it is worth noting that the simulations mainly focused on the prediction of the settling rate and the volume fraction change of the bulk fluid, with the overprediction of the minimum separation distance only slightly influencing the settling rate prediction. Thus, in this study, we used the calculated secondary minima separation distance of the 500:800 system in the simulation, but for 100:500 cases, we used two times the Debye length to replace the minima distance calculated from the DLVO model.

The numerical model was implemented in MATLAB (MathWorks. Inc., version 2016a). A very fine spatial mesh is required near the settling front to accurately capture the solution behavior and prevent numerical instabilities.⁶⁰ The NAG Library (mark 25) was chosen in this study, which includes a collection of general-propose subroutines in MATLAB. The d03pp subroutine, a finite difference-based solver, is used that integrates a system of advection-diffusion equations with the scope to use an adaptive spatial mesh. The numerical implementation is similar to that described in the study of Antonopoulou et al.⁷ who verified the numerical stability and efficiency using various numbers of spatial mesh points.

The sedimentation of bidisperse mixtures of 500 and 800 nm particles (500:800) with a total initial volume fraction ϕ₀ = 0.02, was simulated under Earth gravity conditions for a 1:1 particle ratio. No Earth gravity simulations of the 100:500 nm systems were attempted due to the longer timescales involved. An example volume fraction profile is displayed in Fig. 2 after a time of 11 h, while a second example after 16 h is shown within SM [Fig. S2(a)], where the differences in interfacial position between two species make them both clearly visible. It is noted that the individual species concentration is shown, rather than the total concentration.

The simulated volume fractions for both species rise abruptly in the sediment bed region (near $X / X max=0$⁠), where Fig. S2(b) gives an example of the concentration profiles in the near-bed zone. The overall total volume fraction within the bed approaches the modified bidisperse maximum volume fraction, although it only occupies a very small number of nodes, due to the low dispersion volume fraction. In the supernatant zone above the upper sediment interface, the volume fractions rapidly fall toward zero. In the suspension zones, the situation is more complicated than for monodisperse systems. The suspension consists of two zones: suspension 1 (containing only 500 nm particles) and suspension 2 (containing a mixture of 500:800 nm particles). In suspension 1, the mixture viscosity was reduced to the monodisperse Quemada viscosity model. In suspension 2, the particle size interaction effect was captured by the full bidisperse mixture viscosity model, which led to a small effective viscosity change (and also that the total volume fraction in suspension 2 was double that in suspension 1 and so hindered settling effects were greater). Therefore, the volume fractions in the two suspension areas were slightly different due to the hindrance function changes from bidisperse to monodisperse, which also affected the settling rate. However, the volume fractions remained constant in each zone, as would be expected without high levels of polydispersity,^62,63 and only varied when one size of particle left that zone. There is a clear zonal interface between the small particles in the upper region and the two mixtures.

The two mean interfaces positions (between Suspensions 2 and 1, and between Suspension 1 and the supernatant) are also marked in Fig. 2 and were used to validate the model against experimental sedimentation fronts (see Sec. III C). The volume fraction values used to define the interfaces were half the initial volume fractions of each species (ϕ_i = 0.005). The faction profiles around the interfaces also show a slight “S-shape” curve under Earth gravity, which is due to the diffusive effect of the particles' Brownian motion. Here, the Earth gravity force is not enough to completely overcome diffusion and there remains a relatively wide transition region between zones. This effect is most clear in suspension 1 near the supernatant interface (containing only the 500 nm particles) as would be expected from their smaller size.

Examples of the simulated volume fraction profiles of 500:800 nm centrifugal bidisperse sedimentation using the modified viscosity model are shown in Fig. 3. Three total volume fractions are displayed (ϕ₀ = 0.01, 0.02, and 0.03) with suspensions subject to a centrifugal force consistent with a rotation rate of 1000 rpm (allowing direct validation of the front settling rates against experimental data from the LUMiSizer) with a spatially varying range of centrifugal acceleration from 121.8 g at the upper meniscus to 145.3 g at the base. A mixing ratio of 1:2 (500:800) was chosen for the example to make the distinction in size fractions clear. Profiles are presented at four time periods for the 500 nm species and three times for the faster settling 800 nm particles, spanning most of the total settling period. Similar results for other mixing ratios (2:1 and 1:1) are shown within the SM in Figs. S3 and S4 for the same total volume fractions, while examples of the consolidated bed zones are also given in the SM figures.

It is noted that the reason that ϕ₀ = 0.03 was chosen as the highest total volume fraction simulated was due to the experimental limitations of the LUMiSizer in measuring bidisperse systems at greater concentrations, as discussed in a previous publication by the current authors.⁴⁴ Above this concentration, the transmitted light was totally extinguished in either particle zone, making it impossible to detect an interface between the faster and slower settling species. In contrast, the proposed model would be able to simulate any concentration below the effective maximum packing fraction.

It is observed that, unlike under Earth gravity conditions, there is a small but distinct suspension concentration decrease with time for both particle species (which is most notable for the larger 800 nm particles). In Fig. 3(a), for the large particles, the concentration decreased by ∼13.2% (from 0.0067 to 0.0058 in 360 s), while for the smaller 500 nm species particles, the decrease was ∼12.9% (from 0.00333 to 0.0029 in 960 s). These decreases are similar to previous simulations of centrifugal sedimentation of 500 nm monodisperse particles.⁷ Due to the radially varying acceleration in the setup, particles in suspension closer to the sediment bed have a higher sedimentation velocity than those above. As a result, they are evacuated into the sediment bed at a faster rate than they are replaced from above, which results in a decreasing volume fraction across the suspension. As the height of the sediment bed increases in time, the centrifugal acceleration at the top of the bed is slightly reduced with time, and so the rate of decrease in the volume fraction also declines. Oshima et al.³⁵ in both simulations and calibrated experiments reported similar decreases in concentration for monodisperse and bidisperse colloidal particles under centrifugal acceleration. They also evidenced increased diffusion between interfaces in centrifugal cases, due to the use of an effective axial dispersion coefficient (although the particle sizes were much smaller than those used here).

The relative increase in gravitational acceleration over time also results in an acceleration of the separation velocity, which can be seen from Fig. 3(a) by the spatial distance between each profile front increasing in time (again considered at a relative concentration of each species of 0.5 × ϕ_0i). The arrows in Fig. 3(a) illustrate the distances between profiles, which are seen to increase over time under centrifugal acceleration. It is also noted that the centrifugal acceleration minimizes secondary particle movement due to Brownian motion (since the overall time is much shorter than for Earth gravity) resulting in an almost completely vertical concentration change at the suspension interfaces (and so there is no characteristic “S-type” curve as for Earth gravity). Concentration profiles from the bed zones [SM, Figs. S3 and S4(d)–S4(f)] are similar to the Earth gravity case, with a very distinct and rapid rise toward the maximum packing fraction. There is more evidence of segregation under centrifugation, as would be expected, with a mixed zone containing both species and an upper zone containing only 500 nm particles. There was also some instability in the simulation for the bed fractions in some cases, again as the beds themselves only occupy the first few nodes.

An attempt was made to establish whether the reduction in particle concentrations in the near-bed region could be verified experimentally. Given in Fig. 4 is an example of an attenuation profile measurement (for the lower 0.3 relative height) from the LUMiReader X-Ray for a 1:1 system of 500:800 nm particles (total concentration, ϕ₀ = 0.03). The full profile is given within the SM (Fig. S5). Here, samples were initially accelerated using the analytical centrifuge for set time periods. The profile changes are complicated by the fact that, over time, a drop in attenuation occurs as the particle phases separate, as well as a larger drop at the upper supernatant boundary. However, a focus on the near-bed region does show a consistent reduction in attenuation over time (see Fig. 4). There is a more considerable drop between 400 and 550 s, where the larger particles are completely sedimented, but beyond this point, the attenuation continues to reduce, even though the 500 nm particles are fully suspended in the lower 8 mm section throughout the whole measurement time. Therefore, qualitatively, the centrifugal disruption to concentration correlates well experimentally. Further study is required to enable conversion of attenuation measurements to volume fractions directly via linear calibration measurements.^64,65 Due to the low total volume fractions used, accurate calibration was difficult to achieve in the present system with low noise, and future work will assess the particle depletion effect in more concentrated systems.

Mixtures of 100 and 500 nm particles (100:500) were also simulated, as shown in Fig. 5, for a particle volume ratio of 2:1 and a total volume fraction, ϕ₀, of 0.03, under centrifugation at 4000 rpm. Similar data for mixing ratios of 1:1 are presented within the SM (Fig. S6). Previous experimental studies of the same systems indicated that high centrifugal speeds of >3000 rpm are required to gain representative data with low noise levels, due to the additional Brownian motion of the smaller 100 nm particles.⁴⁴ One of the resulting limitations of such high centrifuge speeds is the large disparity in timescale of the two species, due to the greater size ratio. Indeed, here the larger 500 nm particles in suspension settled within just 77 s, while the smaller 100 nm particles settled in a similar time to the smaller particles in the 500:800 simulations (approximately 1000 s).

Concentration decreases of both species over time are also observed in Fig. 5, in a similar trend to the 500:800 nm case. For the small (100 nm) particles, the concentration decreased by ∼10.9% (in 960 s), while the larger (500 nm) particle concentration decreased ∼12.7% (in 60 s). Again, the reduced concentrations also correlated with an increase in the sedimentation rate over time (with distance differences for each time step again shown within the arrow). Additionally, the influence of particle Brownian motion on diffusion at zonal interfaces was minimal. It is noted that the calculated Peclet number for the 100 nm particles under 4000 rpm is 0.164, which means that Brownian motion is important to be captured in simulation for these smaller particles, although the centrifugal force still dominates.

In Fig. 6, experimental and model interface height profiles under Earth gravity are compared for a 500:800 nm suspension of initial total volume fraction ϕ₀ = 0.02, and a mixing ratio of 1:1. The experimental data were collected from the Turbiscan using a backscatter sensor. The model data were obtained in the corresponding conditions under the same Earth gravity. Here, the “lower” interface represents the distinction between the lower zone containing a mixture of 500 and 800 nm particles and the zone containing only 500 nm particles. The “upper” interface is the distinction between the 500 nm suspension and the particle-free supernatant. While each particle size can be tracked independently in the simulation, it is assumed experimentally that the lower observable interface approximates the sedimentation rate of the larger 800 nm particles.

Comparisons between the experimental and simulation data sets are very close, inferring that experimentally the technique can track both particle species accurately, and most importantly, the underlying physics are correctly captured in the model. However, there are some small deviations in the height vs time profiles. The settling velocity deviation percentage between experiment and the model was averaged, showing 4% for the smaller (500 nm) species and 12.4% for larger (800 nm) species, with an overprediction in the separation rate by the model in both cases. There is a greater difference at the initial period between two data which may be caused by the particle adsorption at the liquid–gas meniscus. The fact that the simulation overpredicts the settling rates is interesting, as it indicates that the mixture viscosity model may be underrepresenting hindered settling effects, to some degree.

Nonetheless, some of the deviations found with the larger species may be due to experimental uncertainty and limitation. At the early time periods, it was found that the experimental data did not extrapolate linearly from a zero-bed height for the 800 nm front (which can be observed in Fig. 6) although the total settling times of the larger particles were very similar between the experiments and simulation. This small deviation implies that there may have been some potential error or uncertainty in the experimental measurements. In fact, while the lower interface was able to be measured from a difference in the backscatter signal, it was not always evident which was the most appropriate strength value to use as the interface (owing to the diffuse nature of the interface). Other effects, such as secondary flows formed by sample mixing on preparation, may have also affected separation at early time periods. Additionally, the small deviation may have been caused by setup differences in the simulation. The height of the cell in the simulation was smaller than that of the experimental height, due to the computational cost (with a 2 cm total height being used for the simulation, against the experimental sample cell of 4 cm). Despite these differences, the correlation is close enough to give high confidence in the model, allowing front height validation to be extended to the centrifugal cases.

Figure 7 presents the comparison of experimental and simulated front height vs time data for the 500:800 bidisperse suspensions, at a mixing ratio of 1:2, under a 1000 rpm at three total volume fractions (ϕ₀ = 0.01–0.03). Comparison data for mixing ratios of 2:1 and 1:1, for the same three volume fractions, are given within the SM (Figs. S7 and S8). There was a slight difference in the suspension height between the experiments and simulations, due to different initializations for the convenience of the numerical mesh (experimentally, total effective heights were about 21 mm, while they were 20 mm in the simulations). Here, the normalized height is displayed, as calculated from the ratio of interface position to the total cell height. Thus, simulation results can be directly compared with the experimental data.

Generally, the simulation results match well with the experiments, although there were again deviations, with the model overpredicting both interfaces (and so consistent with the gravitational results for the same system). Interestingly in this case, unlike in the gravitational system, the differences between the simulation and experimental data were similar for both larger and smaller particle interfaces, averaging at ∼8%–12% overprediction in all cases, with no clear effect of concentration evident in the differences. It is noted that the experimental data for the lower (800 nm) interface were difficult to resolve experimentally for the highest volume fraction case, with a high degree of uncertainty, due to the difficulty measuring the lower and upper zones in the cases of low transmission. This problem also led to difficulties determining the lower interface near the base of the same in the other cases. Indeed, for the 1:1 and 2:1 particle mixtures with a larger proportion of smaller particles in the upper zone, it was not possible to extract lower interface data experimentally at all for the higher concentrations (see SM, Figs. S7 and S8).

Despite these limitations, the centrifugal data do suggest a similar level of accuracy for the monodisperse and bidisperse mixture viscosity models, and additional differences evident in the gravitational case for the larger particle interface may, therefore, be related more to experimental uncertainty. Also, if one considers the relative error of the upper interface, it appears the effect of centrifugation may accentuate any limitations of the model in accurately describing the underlying physics (with an error increasing from ∼4% to ∼10% on average). Essentially, it appears the experimental data display a slightly higher level of hindrance to settling than the mixture viscosity model predicts.

The reason for the discrepancy between Earth gravity and centrifugal conditions may be from a number of sources. It may be that the minimum interparticle distance from the electric double layer is in fact greater than the calculated value (leading to a lower maximum effective volume fraction and increased relative hindrance effects). Alternatively, it may be due to the influence of dynamic effects that are not currently incorporated into the model. For example, the faster particles mean the velocities in the surrounding fluid are higher, which leads to more secondary flows. Under centrifugal conditions, the low Re_p assumption for the Stokes settling velocity may not be strictly valid. In particular, previous studies on bidisperse colloidal sedimentation have observed effects such as dispersion anisotropy caused by particle wakes, or backflows, that enhance drag from the larger particle fraction as it settles through the smaller particle zone.^66–68 Numerical simulations⁶⁹ illustrate the complex nature by which spheres may be affected from the motion of fluid that they disturb, or twinning effects, such as smaller particles being pulled by larger particle wakes. While most studies would suggest these begin to occur in monodisperse systems for Re_p > 1 (which is greater than for any of the particles simulated), we believe that dispersion anisotropy effects are more pronounced in bidisperse cases, where one particle phase settles through another.⁴⁴ The fact that the comparisons were much closer under Earth gravity (with considerably smaller Re_p) would also infer secondary hydrodynamic effects may be the cause of the deviation. It does also appear that the mixture fraction influences the inconsistencies as well, where the differences between the experimental and simulated upper interfaces are larger with an increasing proportion of smaller particles that the larger phase settles through (e.g., when one observes Fig. 7 in comparison with Figs. S7 and S8).

An example comparison of experimental data and simulation results for the 100:500 nm dispersion is shown in Fig. 8, for a mixing ratio of 2:1 and a total volume fraction, ϕ₀, of 0.03. Data for a mixture ratio of 1:1 at total volume fractions of 0.01–0.03 are also presented within the SM (Fig. S9). It is obvious that the 500 nm particles settle much faster than in the 500:800 case, due to the larger centrifugal speed of 4000 rpm. This higher centrifugation, as well as the greater size ratio, meant that there was quite a restricted time where both particle fractions were present in the suspension together, which potentially reduced the overall influence of bidisperse interactions. Additionally, there was a lack of resolution in the upper 100 nm interface for some systems, where there was evidence of adhesion to the meniscus (as it is known that attractive forces may potentially adhere particles to the cell walls⁷⁰). Collectively, this led to some small jumps in the interface data for the upper 100 nm interface [e.g., Fig S9(a)]. However, for most systems, meaningful interface velocities could still be extracted.

As an addition, there were two stages of sedimentation considered separately for the 100:500 cases, because of the much larger size ratio. First, both species exist in the suspension for a small time period; second, only the 100 nm species are left after the fast sedimentation of 500 nm particles is complete. For the 100 nm particles, this two-step process meant that two separate sedimentation velocities were able to be measured both experimentally and from the simulation, equating to those before and after the 500 nm particles settled using linear fitting, as shown in Table II.

Comparing the simulations and experiments, the calculated deviation from the averaged linear settling rate for the larger 500 nm particles was only 0.37%, while for the smaller 100 nm particles settling together with the 500 nm fraction the deviation was a larger 24.9% reducing to 10.5% for the section where the 100 nm species settled alone. However, it must be noted that experimentally this equated to only the first few data points, and so there is likely to be a relatively large degree of uncertainty associated with the experimental data for the initial mixed period. Interestingly, the simulated velocity of the 100 nm particles increased after the 500 nm particles settled, which was the same as the experimental data but to a lesser degree (see Table II). In the simulation, the velocity change is caused by the effective viscosity decreasing due to the lack of the 500 nm species, as well as the increase in the centrifugal force as the settling front progresses. It may be that the deviations between velocities for the 100 nm particles are because the double layer is under-represented for these small particles, and so the effective volume fraction may be slightly larger. Interestingly, the simulation also does not overpredict the settling rate of the smaller species, as in the case of the 500:800 nm system. This result may further suggest that additional hydrodynamic interactions between the species in the 500:800 nm mixtures may be the reason for enhanced hindered settling, as the greater size ratio for the 100:500 nm systems means that the species only interact together for a very small period.

The average bidisperse settling velocities for the 500:800 nm mixtures from the simulations are plotted in Fig. 9 for three initial volume fractions ranging from 0.01 to 0.03 and three mixing ratios (1:2, 1:1, and 2:1). The data are compared to a number of experimental measurements over a similar concentration region and are also plotted against Richardson–Zaki (R–Z) predictions,⁵ using an exponent of n = 5.1.

There are two key trends that are evident from this comparison. First, both the smaller and larger particles compare well with the R–Z prediction in almost all systems, apart from the larger species at the highest 0.03 total volume fraction simulated. This slight deviation may suggest that the R–Z settling exponent may be a slightly lower value for the larger particles (which could be expected, given that the 800 nm species are toward the upper limit for colloidal systems, where the double layer is not significant). Also, both the simulations and the R–Z predictions deviate from the experimental data to varying degrees, depending on the mixture ratio. Here, the predictions for the larger size fraction compared to the experimental data are weakest for the 2:1 mixture ratio (with the highest proportion of smaller particles), while they are weakest for the smaller species with the greatest proportion of larger particles. As evident in the height vs time predictions, the R-Z prediction is consistent with the simulations in overpredicting the settling rates in both cases. Therefore, these results would agree with the hypothesis that additional hydrodynamic effects are present as the larger particle fraction settles through the smaller particles, where neither the mixture viscosity model in the simulations nor the empirical R–Z model takes these secondary effects into explicit consideration. It is stressed, however, that the deviations to the experimental data are minor (<10% in all cases) and it would be extremely computationally intensive to try to incorporate hydrodynamic mixing within the simulation, as the tracking of individual particles is probably required. However, the current work is assessing whether the bidisperse mixture viscosity model may be improved to implicitly improve the correlations.

The averaged bidisperse settling velocities of the 100:500 case from the simulations are plotted in Fig. 10, again in comparison with the R–Z model and experimental data. Here, only a mixing ratio of 1:1 was used, due to the much greater computational time required to simulate these systems of smaller nanoparticles. In this case, as shown from the single interface vs time plot (Fig. 6), the comparison between the simulations and experimental data is generally closer than for the 500:800 systems. It is also clear that the prediction is again very closely aligned to the empirical R–Z prediction. It is assumed that the better performance of the simulation in the 100:500 nm cases is due to the large size ratio, where each particle fraction is not in contact for a long period of time (as discussed). Hence, the average velocities from the slower 100 nm particle phase are largely taken from a monodisperse suspension. For the larger 500 nm phase, there is also no clear reduction in the experimental settling rate, within error, that would infer additional hydrodynamic effects. It may be such effects could become more prevalent at other mixture ratios, but it is likely the larger size ratio itself, and the small contact time under high centrifuge reduces secondary interactions.

In this study, a 1-D advection-diffusion sedimentation model for bidisperse colloidal suspensions was proposed, which uses a bidisperse mixture viscosity model to consider the effect of hindered settling and particle interactions. The model also modified the hard sphere assumption to account for colloidal interactions by using an effective maximum volume fraction, which was computed based on the DLVO theory. Simulations were validated using the experimental data and the empirical Richardson–Zaki model, using two different size ratios (of 500:800 nm and 100:500 nm particles) at different mixing ratios and total concentrations. Results highlighted that interparticle interaction is an important aspect in determining the sedimentation velocities of colloidal suspensions, especially for bidisperse systems. The inclusion of a minimum interparticle distance is a significant addition to sedimentation models and will be especially applicable at high volume fractions. The use of a 1-D advection-diffusion model allowed the important physics of particle separation to be captured, without requiring a full computational fluid dynamics approach. Also, the adaptive spatial mesh was shown to be numerically stable while giving sufficient resolution around zonal interfaces to simulate the influence of diffusion on concentration.

The proposed model was validated under both Earth gravity and various centrifugal force conditions. It was found that the predicted volume fraction of both species reduced in time under centrifugation (due to the increase in force at the bottom of the simulated dispersion cell), whereas they remained almost constant at Earth gravity, as concentrations were not great enough to force physical interaction between the size classes. The bidisperse simulation was also used to track the suspension interfaces of each size fraction and, therefore, gain an average estimated settling rate. For the 500:800 nm mixture cases, deviations in average settling rates between simulations and experimental data were <10% in all cases. However, there was a clear trend of the simulations slightly overpredicting settling rates of the particle phase with the lowest mixture fraction. These results suggested that hindrance effects in bidisperse systems are higher than in the mixture viscosity model prediction, which may be due to secondary hydrodynamic factors that are not explicitly accounted for. Settling rate validation for the 100:500 nm cases was generally closer than for the 500:800 cases, as the centrifugal force was larger, resulting in the interaction time between the species being smaller and reducing secondary hydrodynamic interactions. The proposed model could also be extended to a true polydisperse system, with appropriate modification to the viscosity to consider the interactions of one species with all others. Importantly, due to the importance of defining the effective maximum volume fraction of multiple species, simulations of particle packing may offer the best route toward polydisperse extension. Therefore, the model may be used to predict concentration changes and settling rates of many systems undergoing centrifugal sedimentation and separation.

See the supplementary material for example calculation of maximum volume fraction and relative viscosity [Eqs. (S1–S4)]; independence grid verification (Fig. S1); example of Earth gravity sedimentation at 16 h, including near-bed zone (Fig. S2); simulated bidisperse concentration profiles of 500:800 nm mixtures under a 1000 rpm centrifuge and a mixing ratio of 1:1 (Fig. S3) and 2:1 (Fig. S4), both for three total volume fractions, including near-bed zones; normalized x-ray sedimentation profile for a 500:800 nm mixture, a mixing ratio of 1:1, and a total volume fraction of 0.03, at five time periods (Fig. S5); simulated concentration profiles of 100:500 nm mixtures under a 4000 rpm centrifuge and mixing ratio of 1:1 and three total volume fractions (Fig. S6); height vs time profiles for 500:800 nm mixtures at three total volume fractions, and volume ratios of 1:1 (Fig. S7) and 2:1 (Fig. S8); And height vs time profiles for 100:500 nm mixtures at three total volume fractions and a volume ratio of 1:1 (Fig. S9).

The authors would like to thank the Engineering and Physical Sciences Research Council (EPSRC), UK, for funding and access to the x-ray sedimentation system, as part of the MULTIphase Fluid flOw in nucleaR systeMs (MULTIFORM) National Nuclear User Facility (Grant No. EP/V034898/1). T.C.S. was also supported by EPSRC through Grant No. EP/L01615X/1. Dr. Ben Douglas is thanked for experimental support.

The authors have no conflicts to disclose.

Hangyu Chen: Conceptualization (lead); Formal analysis (lead); Methodology (equal); Visualization (lead); Writing – original draft (lead). Thomas C. Sykes: Methodology (equal); Software (lead); Writing – review & editing (equal). Oguzhan Kivan: Investigation (supporting); Validation (supporting); Writing – review & editing (equal). Xiaodong Jia: Supervision (equal); Writing – review & editing (equal). Michael Fairweather: Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). Timothy N. Hunter: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are openly available in University of Leeds Data Repository at https://doi.org/10.5518/1421.⁷¹