Numerical and experimental analysis of the sedimentation of spherical colloidal suspensions under centrifugal force

Suspensions of colloidal particles can be found in many industrial applications, such as consumer products, surface coatings, and printer inks.^1,2 Examples of such consumer products include fabric enhancers, cosmetics, and shampoos.³ In many of these applications, the dispersed phase has a higher density than the continuous phase. Therefore, sedimentation of the dispersed phase may occur due to gravity. Sedimentation leads to the separation of the continuous and dispersed phases (phase separation), resulting in the formation of a particle sediment bed. Similarly, phase separation can occur when the dispersed phase is less dense than the continuous phase, with the dispersed phase rising to the top in a process known as creaming. A crucial requirement of colloidal suspensions in many applications is that noticeable phase separation must not occur during the product’s shelf life.⁴ Some products are formulated as a space-spanning network, or gel, in an attempt to prolong their shelf life. However, even gels can be mechanically unstable.⁵ Other products must remain as suspensions, which motivates the study of sedimentation of colloidal suspensions in this work.

This work considers monodisperse suspensions of spherical colloidal particles without particle aggregation. For such a suspension, in which the dispersed phase is more dense than the continuous phase, three components will form due to gravity.⁶ A clear region forms above the suspension (the supernatant) and a sediment bed with a high volume fraction forms below the suspension. The volume fraction of the suspension component remains constant at the initial volume fraction due to the constant gravitational acceleration across the suspension. There is a rapid change in volume fraction between the supernatant and suspension, called the supernatant-suspension interface. Polydispersity in the particle size, and Brownian motion of colloidal sized particles, may cause the thickness of this interface to increase in time.⁷ Sedimentation continues until all of the dispersed phase is contained within the sediment bed, with only clear fluid above. When this is the case, complete phase separation has occurred.

The first mathematical treatment of sedimentation which related batch sedimentation rates to the microscopic properties of the suspension was provided by Kynch,⁸ who proposed a one-dimensional monodisperse sedimentation model based on the particle flux through a given height of the sample. One-dimensional models are appropriate for describing sedimentation since the particle volume fraction is one-dimensional in height.⁶ Kynch’s model was extended by Davis and Russel⁹ to include diffusion due to Brownian motion for colloidal particles. The addition of the diffusive term to take account of Brownian motion allows the sediment to be included in the model, rather than relying on a time-dependent boundary condition at the top of the sediment. The resulting particle volume fraction around the supernatant-suspension interface is a balance between diffusive front spreading and sharpening due to hindered settling effects.¹⁰ By including the sediment within the model, compression of the sediment can be taken into account, using a compressive yield stress.¹¹ Such one-dimensional models have also been successfully extended to account for the effects of polydispersity in particle size,¹² and density stratification in the continuous phase.¹³ 

Particles in suspension do not sediment alone, and as a result their motion is hindered by surrounding particles, where the effect of hindrance increases as the particle volume fraction increases. The effect of hindrance is taken into account by a hindrance function, usually^9,12 the Richardson-Zaki empirical relationship.¹⁴ Other power-law expressions have been proposed based on the Richardson-Zaki relationship¹⁵ and implemented in sedimentation models.¹⁶ However, for monodisperse suspensions, in which the maximum volume fraction is not one, Richardson-Zaki power-law expressions can become inaccurate by allowing unphysical volume fractions.¹⁷ Alternatively, particle hindrance can be taken into account by using a relationship based on a suspension viscosity model.¹⁸ Using a suspension viscosity model has the advantage of reducing the sedimentation velocity to zero at the maximum volume fraction. Such a formulation has been successfully implemented in a two-dimensional sedimentation model by Rao et al.⁶ 

Much of the sedimentation literature neglects interactions between the particles in suspension, though hydrodynamic interactions for sedimentation at moderate Reynolds number have been considered.¹⁵ A common assumption made is the hard-sphere assumption, an idealisation which assumes no interparticle interactions other than infinite repulsion at contact.^19,20 This means that there is no restriction on the minimum separation distance between two colloidal particles, so the maximum volume fraction is determined by random sphere packing theory in this case. Whilst this assumption may hold for larger non-colloidal particles,²¹ it may not be suitable for highly charged colloidal particles.

Genovese¹⁹ has validated the use of an effective volume fraction to describe interparticle interactions and thus determine the shear rheology of colloidal suspensions. An effective volume fraction is defined based on an effective particle size which is the sum of the particle’s physical size and an exclusion volume around the particle due to interparticle interactions. An effective diameter has also been used to describe depletion interaction.²² Metin²³ has demonstrated the use of an effective particle size for describing colloidal suspensions, but taking into account the interparticle interactions using an effective maximum volume fraction, which has the advantage of remaining constant for a given suspension. Metin’s approach can be applied to account for interparticle interactions in a sedimentation model, as demonstrated in this work.

Experimental data can be used to characterise the stability of suspensions and validate sedimentation models. Many sedimentation experiments use custom charge-coupled device (CCD) camera arrangements with image processing software.^1,5 The aim is usually to track the supernatant-suspension interface and obtain its velocity, known as the separation velocity, v(t). The separation velocity characterises the rate of phase separation. Other techniques, such as magnetic resonance imaging,⁴ nuclear magnetic resonance,^6,12 and acoustics²⁴ can determine the local particle volume fraction. Such experiments tend to use a static sample, so are subject to Earth gravity. Centrifugation may be used to accelerate sedimentation, by effectively replacing the gravitational acceleration with a centrifugal acceleration. Centrifugation is useful since sedimentation experiments are often very lengthy, especially for consumer products which are formulated to have long shelf lives. The time required for complete sedimentation can be reduced by over an order of magnitude using centrifugation. Instruments are available to monitor transmission through the sample in time in a centrifuge.²⁵ The separation velocity under centrifugal force can be determined using these transmission profiles, and the corresponding Earth gravity separation velocity can be determined. However, recent studies have cast doubt on the accuracy of recovering Earth gravity settling velocities from centrifuge experiments.²⁶ We hypothesise that the difficultly in recovering Earth gravity sedimentation rates is due to the variation in centrifugal acceleration across a sample in a centrifuge since the centrifugal acceleration depends on the radial distance from the centre of rotation, which is an effect typically not directly taken into account in Earth gravity conversions.

In this work, a mathematical model for the sedimentation of monodisperse colloidal particles has been developed. Hindered settling effects are taken into account in a one-dimensional sedimentation model⁹ using the suspension viscosity,⁶ rather than a hindrance function like the Richardson-Zaki relationship. By introducing an effective maximum volume fraction, a method which has been successfully validated for describing the shear rheology of colloidal suspensions in the past,²³ interparticle interactions are taken into account in the model. The numerical implementation of this model using an adaptive finite difference solver is described. An extension to the model for sedimentation under centrifugal acceleration is proposed which allows model results to be directly compared to unconverted experimental data for nearly monodisperse colloidal silica suspensions collected under centrifugal force at 705 rpm, corresponding to a centrifugal acceleration of between 59.6g and 71.8g across the sample suspension. An investigation into the effects of centrifugation on sedimenting suspensions is described, which leads to the validation of a proposed method for successfully converting centrifugal force sedimentation data to Earth gravity conditions and an improved understanding of the behaviour of colloidal suspensions under centrifugation.

Spherical silica nanoparticles (AngstromSphere, labeled 500 nm) were purchased as a dry powder from Fiber Optic Center (New Bedford, USA). The particle density was measured as 1.92 g cm⁻³ using a Pycnomatic ATC gas pycnometer (Thermo Electron, USA). A known mass of the silica powder was added to a calibrated volume, with the displaced volume determined by monitoring the pressure change of helium gas, and the mass measured using an analytical balance (0.1 mg precision).

Suspensions were prepared by adding silica powder to a 0.01M Potassium Chloride (KCl) solution. This solution was prepared using KCl crystalline powder (Fluka Chemie GmbH, Germany) and ultrapure water, with resistivity 18.2 MΩ cm at 298.15 K (Milli-Q, Millipore, USA). Before each sedimentation experiment, the suspension was magnetically stirred for 5 min, then placed in an ultrasonic bath for 5 min, and finally magnetically stirred for a further 25 min. The purpose of magnetic stirring was to homogenise the suspension, thus ensuring it was well-dispersed. The ultrasonic bath was used to break up any particle aggregates, to ensure the suspension was monodisperse.

A Malvern Zetasizer Nano ZS ZEN3600 (Malvern Instruments Ltd., UK), with a 632.8 nm, 4 mW laser was used to measure the size of the particles in suspension using dynamic light scattering. Approximately 1.2 ml of suspension was added to a 12 mm square polystyrene cuvette (DTS0012, Malvern), using a 1 ml syringe. For the silica dispersed phase, the refractive index used was 1.44, and the absorption 0.001. For the continuous phase, the refractive index used was 1.33 and the dynamic viscosity was 8.9 × 10⁻⁴ Pa s.²⁷ For each sample, 3 measurements were made, with 14 sub runs (division of each size measurement) each. The average number-mean particle size was approximately 480 nm. This value varied by approximately 10 nm across the different samples measured. A typical particle size distribution for a sample, showing three separate measurements with different line types, is shown in Fig. 1. There is a single peak at the primary particle size, with a relatively tight distribution around the peak. Further, the polydispersity index (PDI), a measure for the width of the particle size distribution, was calculated using the Zetasizer. The measured PDI of 0.069, with a standard deviation of 0.008, confirms that the suspension is nearly monodisperse.²⁸ 

The Zetasizer was also used to measure the zeta potential of the suspension by determining the electrophoretic mobility of the dispersed phase and then applying the Henry equation. 1.0 ml of the suspension was added to a folded capillary cell (DTS1070, Malvern), using a 1 ml syringe. Care was taken to ensure bubbles did not form in the sample. Five zeta potential measurements were made. The mean average zeta potential was −43.5 mV, with a standard deviation of 2.3 mV.

Phase separation was monitored using a LUMiSizer analytical centrifuge (LUM GmbH, Germany). The LUMiSizer accelerates sedimentation by exposing a sample suspension to a centrifugal acceleration greater than the Earth’s gravitational acceleration. Up to 12 samples can be analysed simultaneously. For each sample, approximately 430 μl of suspension was extracted with a 500 μl autopipette and transferred into a 2 mm rectangular polyamide cell (LUMiSizer cell type 3), resulting in a sample height of approximately 22 mm. Each sample was sealed and placed horizontally in the LUMiSizer, along the radial axis, with the top of the sample nearest to the centre of rotation. The base of the LUMiSizer sample was determined to be at a radial distance of 129.3 mm from the centre of rotation, by analysing light transmission profiles of ultrapure water. For each suspension investigated, three repeats were made. Samples were analysed at centrifuge rotation rates of 705 rpm and 1000 rpm, at a temperature of 25 °C (298.15 K). In the LUMiSizer, a light source pulses near-infrared (865 nm) light through the side of each sample cell at user-specified times. The light intensity was normalised prior to each run. A 25 mm 2048 element CCD-line detects the intensity of transmitted light across the length of the sample, yielding transmission profiles. A profile was collected every 10 s, until the sample had fully sedimented. Each transmission profile was automatically smoothed by the instrument’s software using a 9 point moving average. 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. Since the suspensions are nearly monodisperse, there is a large and rapid change in transmission intensity near the supernatant-suspension interface in each transmission profile. The interface height for each profile was determined in the “front tracking” module of the LUMiSizer software (SEPView 6.2). A transmission threshold defines the transmission intensity corresponding to the interface. The interface height was found to depend weakly on the transmission threshold chosen since the suspension was nearly monodisperse. Therefore, 50% transmission intensity was used to represent the interface position. The final output from the LUMiSizer is a time series of interface heights obtained under centrifugal acceleration for each sample investigated. The experimental conditions for the sedimentation experiments are summarised in Table I.

A sedimentation model for monodisperse, colloidal suspensions is proposed, based on the model of Davis and Russel.⁹ The particle volume fraction, ϕ at a height x and time t, is described by the advection-diffusion equation

where u(ϕ) is the hindered settling velocity of a particle in suspension (positive in the upward direction) and D(ϕ) is the diffusion coefficient. The continuous phase is assumed to be stagnant, except for displacement due to sedimentation. This allows the application of the zero-flux boundary conditions

where x = 0 at the base of the sample and x = x_max at the meniscus. The suspension is assumed to be initially well-dispersed, so the initial condition is

where ϕ₀ is constant. The hindered settling velocity of a particle in suspension differs from its terminal velocity since surrounding particles obstruct and restrict the backflow of the surrounding continuous phase. The hindered settling velocity is given by

where f_h is a function describing the hindrance and u_T is the particle terminal velocity. Since sedimentation problems are in Stokes’ regime, the terminal velocity is given by

where d is the dispersed phase diameter, ρ_d is the dispersed phase density, and g is the Earth gravity acceleration (a modification for centrifugal acceleration is described below). ρ_c and μ_c are the density and viscosity of the continuous phase, respectively. For ultrapure water, ρ_c = 997.1 kg m⁻³ and μ_c = 8.9 × 10⁻⁴ Pa s.²⁷ A commonly used hindrance function is the Richardson-Zaki relationship,¹⁴ 

where p is a constant generally determined experimentally. A significant disadvantage of using Eq. (6) in sedimentation models for colloidal particles is that it does not reach zero at the maximum volume fraction and so allows the volume fraction to become unphysically large. Instead of Eq. (6), this work uses the closed system version of the Acrivos hindered settling function,⁶ 

By using Eq. (7), the mixture viscosity of the suspension, μ is used in the model, rather than the continuous phase viscosity, μ_c, thus taking into account hindered settling effects. The mixture viscosity model can be chosen such that it diverges as the maximum volume fraction is approached, so the hindered settling velocity becomes zero there. Further, the use of Eq. (7) negates the need to specify p in Eq. (6) which is a parameter that takes many different values in the literature.

Due to interparticle interactions, the maximum volume fraction of the dispersed phase is not necessarily well described by random sphere packing theory for colloidal particles. Electrostatic repulsion provides an energy barrier which gives rise to a non-zero minimum separation distance between particles. An effective maximum volume fraction, $ϕmaxeff$ is defined based on the effective diameter of the particles, with the minimum separation distance taken into account. The Quemada model with the effective volume fraction¹⁹ is used for the mixture viscosity model in Eq. (7),

The Quemada model has been shown to model the mixture viscosity of silica suspensions well.³⁰ Equation (8) allows the hindered settling velocity to become zero at the effective maximum volume fraction and takes into account interparticle interactions. This eliminates the hard-sphere assumption.

In a similar manner to advection, diffusion is hindered by surrounding particles. Therefore, the mixture viscosity replaces the continuous phase viscosity in the Stokes-Einstein equation to yield the diffusion coefficient¹² 

where k_B is Boltzmann’s constant and T is the absolute temperature. Equation (9) is an approximation of the diffusion coefficient proposed by Kholodenko and Douglas,³¹ in the low volume fraction limit.

The effective maximum volume fraction, $ϕmaxeff$ is given by²³ 

where $s¯$ is the minimum separation distance between the surfaces of two particles and ϕ_max is the low shear maximum volume fraction for randomly packed hard-spheres. Here, ϕ_max = 0.639 is assumed, as determined by Qi and Tanner.³² Note that the sediment bed compresses towards the maximum volume fraction as increased pressure is applied. Hence, the maximum volume fraction is unlikely to be fully realized.

The minimum separation distance can be calculated using the DLVO (Derjaguin-Landau-Verwey-Overbeek) theory.^33,34 DLVO theory provides a quantitative framework to determine the interaction potential of two particles based on attractive (van der Waals) and repulsive (electrostatic) intermolecular forces. The gravitational potential is assumed to be negligible in the calculation of the interaction potential. For suspensions with a large interaction potential relative to the thermal energy, random collisions between particles due to Brownian motion are unlikely to result in aggregation,³⁵ giving rise to a non-zero separation distance between the particles. In this case, the minimum separation distance may correspond to the location of the second minimum in the interaction potential, for the interacting particles to achieve the minimum energy state. The separation distance corresponding to the secondary minimum is taken to be the minimum separation distance, $s¯$ in this work.

To identify the minimum separation distance, DLVO theory is used to determine the interaction potential of the particles, from which the secondary minimum separation distance can be identified. The total interaction potential, U is given by³⁰ 

where U_VDW is the attractive (van der Waals) component and U_EDL is the repulsive (electrostatic) component of the interaction potential.

The attractive (van der Waals) component of the interaction potential for two spherical particles of diameter d is given by^36,37

where A_H is the Hamaker constant for the colloidal particle interacting through the given medium. For the silica-water-silica system studied, the Hamaker constant is taken to be A_H = 2.5 × 10⁻²⁰ J. This choice is discussed in Sec. V A.

The repulsive (electrostatic) component of the interaction potential is given by^37,38

where k_B is the Boltzmann constant, T is the absolute temperature, N_A is the Avagadro constant, I is the ionic strength of the electrolyte (mol m⁻³), γ is the reduced surface potential, and κ is the inverse Debye length. Equation (13) is constructed on the basis of Derjaguin’s approximation, which has been shown to be valid for spherical colloidal particles.³⁹ Derjaguin’s approximation requires s ≪ 0.5d, which is satisfied here since DLVO theory is only used to identify the minimum separation distance of colloidal particles. γ is given by³⁸ 

where e is the elementary charge. Rather than the surface potential, ζ is the zeta potential, which is closely related to the surface potential and easily measurable.³⁰ For monovalent electrolytes (like KCl), the inverse Debye length is given by⁴⁰ 

where ε₀ = 8.85 × 10⁻¹² F m⁻¹ is the permittivity of free space³⁰ and ε_r is the relative permittivity of the medium. In the case of ultrapure water at 298.15 K, ε_r = 78.304.⁴¹ 

Numerically, the location of the secondary minimum is found by using a linear difference scheme with picometer resolution to identify the second sign change in the gradient of the DLVO-computed interaction potential (the first sign change corresponds to the primary maximum). The secondary minimum is identified as the node between the differences exhibiting the sign change.

Instead of DLVO, Genovese¹⁹ uses the Debye length κ⁻¹ to approximate the minimum separation distance, taking $s¯≈2κ−1$⁠. This approximation may lead to an underprediction of the minimum separation distance, as κ⁻¹ is dependent only on the electrolyte concentration.

Experimental validation of the numerical value of the effective maximum volume fraction was attempted. Whilst it is theoretically possible to determine the effective maximum volume fraction directly by experiment, the size of colloidal particles and the resulting small height of the sediment gives rise to a large experimental error. Therefore, the DLVO-based method described in this section is the recommended way of determining the effective maximum volume fraction in sedimentation modeling.

The sedimentation model can be adapted for centrifugal force conditions. This allows direct model comparisons to LUMiSizer experiments. To do this, the gravitation acceleration, g in Eq. (5), is replaced by the centrifugal acceleration, a, to give

which is used in Eq. (4), in place of Eq. (5), to define the hindered settling velocity. The centrifugal acceleration, a, is given by

where n is the number of centrifuge revolutions per second and r is the radial position in the centrifuge, with r = 0 at the centre of rotation. Due to the position of the sample in the centrifuge, r is given by

where r_max is the radial position at the base of the sample. Equation (18) shows that the centrifugal acceleration varies with height in the sample. Therefore, by Eq. (16), the terminal velocity of a given particle decreases with height in the sample. Due to the LUMiSizer dimensions, the centrifugal acceleration varies by approximately 19% across the sample cell in the centrifuge used for the reported experiments.

To summarise the model, the form of Eq. (1) proposed is

where Γ is g in Earth gravity conditions and (2πn)² (r_max − x) in centrifugal force conditions. The effective maximum volume fraction is given by

where $s¯$ is taken to be the separation distance corresponding to the secondary minimum in the DLVO-computed interaction potential between two colloidal particles.

The model is implemented in MATLAB⁴² (version 2016a) using a finite difference method, motivated by the code of Beeckmans et al.⁴³ Due to the rapid change in volume fraction at the supernatant-suspension interface, there is a travelling front in the solution. A very fine spatial mesh is therefore required near the front and in the sediment to accurately capture the solution behaviour and prevent numerical instabilities.⁴⁴ For regions where the volume fraction is relatively uniform, a coarser spatial mesh is sufficient. An adaptive spatial mesh is therefore used to follow the time-dependent nature of the solution, which increases efficiency by reducing the number of spatial mesh points required.

The NAG Library is a collection of general-purpose subroutines, which are callable in MATLAB via the NAG Toolbox for MATLAB (Mark 25).⁴⁵ The d03pp subroutine, a finite difference based solver which integrates a system of advection-diffusion equations with scope to use an adaptive spatial mesh, is chosen to solve the initial boundary value problem. This subroutine reduces Eq. (1) to a coupled system of stiff ordinary differential equations using the method of lines, replacing the spatial derivatives by finite difference approximations.⁴⁶ Time integration is then carried out using a backward differentiation formula,⁴⁷ with a variable time step, determined automatically using a user-specified local error tolerance.⁴⁸ This is taken to be 5 × 10⁻⁵ based on the maximum norm.

The subroutine generates a new spatial mesh every ten time steps, based on the current solution profile, as described by a monitor function, m(x, t). The routine equally distributes the integral of the monitor function over the domain.⁴⁵ The monitor function must therefore be large in the sediment and near the front. In this work,

is used as the monitor function, where K(x, t) is the smooth function curvature estimator

given and discretised by Sfakianakis,⁴⁹ where ϕ′ denotes the first spatial derivative of ϕ. This monitor function uses the first spatial derivative of ϕ in the sediment (since curvature is low there) and the curvature everywhere else. The exponent in Eq. (21) is used to improve numerical stability. A maximum adjacent cell ratio of 1.1 is set, allowing this monitor function to be used.

The effect of spatial mesh on the solution was investigated by comparing solutions computed using various numbers of spatial mesh points, in both Earth gravity and centrifugal force conditions. Based on the maximum resolution of the LUMiSizer, 2400 and 3600 spatial mesh points were deemed sufficient for simulations in Earth gravity and centrifugal force conditions, respectively.

The local error tolerance was varied with a fixed spatial mesh to understand its effect on the solution. Varying it between 1 × 10⁻³ and 1 × 10⁻⁷ changed the solution by less than 0.01%, suggesting that the choice of 5 × 10⁻⁵ does not affect the solution obtained.

All numerical simulations described are made to be comparable with experimental data. The sample height in each case, x_max, is taken to correspond to experimental data. The sample height is not completely consistent across samples due to experimental variation in filling the experimental cells. Refer to the supplementary material for the MATLAB code used for the numerical simulations.

The experimental results obtained are a time series of interface heights, measured in centrifugal force conditions. Each time series needs to be converted to Earth gravity conditions, so that the separation velocity, v(t) in Earth gravity conditions can be determined, which is often the quantity of interest. A key challenge of using centrifugation in suspension stability analysis is to reliably convert data collected under centrifugal force to Earth gravity conditions.

Several such conversions have been proposed in the literature. For example, Lerche and Sobisch²⁵ propose a method based on the initial separation velocity and the relative centrifugal force (RCF) at the meniscus, x_max. The RCF is defined as the ratio of centrifugal acceleration to the gravitational acceleration. A significant disadvantage of their method is that light transmission profiles near the interface are generally subject to noise, which makes the initial separation velocity difficult to determine reliably. Furthermore, Tehrani-Bagha²⁶ demonstrated that extrapolation from experiments using multiple RCF values does not accurately recover the separation velocity at Earth gravity.

A conversion to Earth gravity conditions, based on the space-averaged RCF experienced by the interface, is proposed here. Consider an interface height h_c achieved at a time t_c under centrifugal force. The space-averaged RCF experienced by the interface until the time t_c is the mean of the RCF at the meniscus and the RCF at the radial position corresponding to the interface, and so it is given by

for a given n (number of centrifuge revolutions per second), with one term arising from the meniscus and the other from the interface. Since a depends linearly on r, an interface height h_c (achieved at time t_c in centrifugal force conditions) is achieved at a time

in Earth gravity conditions. Both the interface position and suspension concentration affect the time, t_c which ensures these effects are taken into account by the proposed conversion. For each interface height obtained under centrifugal force, the corresponding time can be converted to an equivalent time in Earth gravity conditions using Eq. (24). Performing this conversion for each available data point, a times series of interface heights in Earth gravity conditions can be recovered from data collected under centrifugal force. As explained in Sec. I, the separation velocity in Earth gravity conditions, v(t) is expected to be constant. Therefore, the constant Earth gravity separation velocity can be determined by fitting a straight line to a graph of interface height against time (in Earth gravity conditions) using linear regression, thereby recovering the Earth gravity result which is often desired.

Note that the proposed conversion can be used to generate Earth gravity settling velocities from data that can be directly extracted from an analytical centrifuge. The validity of this conversion process will be assessed in Sec. V D.

Figure 2 shows the DLVO-computed interaction potential for silica dispersed in 0.01M KCl, consistent with the experimental conditions. The location of the secondary minimum is marked, corresponding to a separation distance of 18.7 nm. This separation distance is taken as the minimum separation distance for the evaluation of the effective maximum volume fraction. The minimum separation distance is well-determined by the location of the secondary minimum in the interaction potential due to energy considerations since closer contact would increase the energy of the interaction. Nevertheless, Fig. 2 shows that this secondary minimum in the interaction potential is fairly broad and shallow, so some variation in the minimum separation distance may be expected. However, it is likely that any variation will average out over the number of particles and therefore that the minimum separation distance is well estimated by the separation distance corresponding to the secondary minimum in the interaction potential, as proposed.

Using the calculated minimum separation distance for the 0.01M KCl suspension, the effective maximum volume fraction is determined to be $ϕmaxeff=0.570$ using Eq. (20). The effective maximum volume fraction is approximately 10% less than the classical maximum volume fraction in this case. This variation leads to a significant difference in predicted separation rates using the effective maximum volume fraction as opposed to the theoretical one, as demonstrated in Sec. V E.

As mentioned in Sec. III B, the value of the Hamaker constant for the silica-water-silica system used was 2.5 × 10⁻²⁰ J. The Hamaker constant for a given system is notoriously difficult to obtain experimentally.⁵⁰ However, many different values for the Hamaker constant of silica-water-silica systems are indicated in the literature, for example, 0.7 × 10⁻²⁰ J and 3.8 × 10⁻²⁰ J.^30,51 As a result, the Hamaker constant used in this work was chosen instead on the basis of a qualitative analysis of the experimental silica suspensions. When silica was dispersed in 0.1M KCl, aggregation was observed within several minutes of preparing the suspension. With a primary maximum in the interaction potential of 10k_BT, rapid aggregation is expected, but with a primary maximum of 50k_BT, relatively slow aggregation is expected instead. A primary maximum of 40k_BT was therefore chosen to appropriately describe the experimentally observed aggregation behaviour of the 0.1M KCl suspension. To achieve this condition, the value of the Hamaker constant was calculated as 2.5 × 10⁻²⁰ J. Figure 2 shows the corresponding interaction potential for silica dispersed in both 0.01M KCl and 0.1M KCl, using the calculated value of the Hamaker constant. The 0.1M KCl suspension has the chosen primary maximum of 40k_BT, whereas the primary maximum of the 0.01M KCl suspension is much larger, indicating almost no aggregation during the experimental period. The choice for the value of the Hamaker constant of 2.5 × 10⁻²⁰ J is therefore physically realistic. Note also that if the actual value of the Hamaker constant were less than what is used here, then the minimum separation distance would increase, thereby increasing the difference between the effective and standard maximum volume fractions, adding further support to the proposal of using an effective maximum volume fraction in sedimentation models.

The sedimentation of a suspension with ϕ₀ = 0.09 was simulated in Earth gravity conditions. A volume fraction profile is shown in Fig. 3 after a time of 30 h, which is approximately half the time required for complete sedimentation. The volume fraction is large in the sediment (near x/x_max = 0), where it approaches the effective maximum volume fraction. In the supernatant, the volume fraction is zero, as the particles have sedimented towards the base of the sample. In the suspension component, the volume fraction is constant. A constant volume fraction in the suspension component is expected since the hindrance function remains constant there, so the settling velocity remains constant. This behaviour is also noted in the literature.^6,12

In Fig. 3, the data are plotted as points, located at the spatial positions where the solution is calculated for the given time step. The spatial mesh points cluster in the sediment and around the travelling front (the supernatant-suspension interface). Fewer spatial mesh points are used in the regions where the volume fraction is relatively uniform, enabling the simulation to be numerically stable whilst improving numerical efficiency.

The travelling front represents the supernatant-suspension interface, with the interface position marked in Fig. 3. The volume fraction which defines the interface is half the initial volume fraction, which is consistent with the transmission threshold used for the experimental data. The volume fraction profile around the interface takes the form of a narrow S-shape curve, rather than the interface being defined by a single discrete value, due to the diffusive effect of Brownian motion for colloidal particles. In contrast, there is a sharp change in volume fraction at the top of the sediment. The sediment height grows in time, and the volume fraction of the sediment bed increases towards the effective maximum volume fraction. The volume fraction in the simulation is bounded above by the effective maximum volume fraction according to Eq. (8), with further compression negligible due to the height of the primary maximum in the interparticle interaction potential (Fig. 2).

Figure 4 shows the model volume fraction profiles at various times for a suspension simulated with ϕ₀ = 0.03. Figure 4(a) displays the results under Earth gravity conditions. The volume fraction of the suspension component remains both constant and uniform, maintaining the initial volume fraction, ϕ = ϕ₀. This is because the suspension is monodisperse, so the hindered settling velocity for all particles is constant at a given volume fraction according to Eq. (4).

Figure 4(b) displays the corresponding volume fraction profiles subject to a centrifugal force consistent with a rotation rate of 705 rpm in the LUMiSizer, with a spatially varying range of centrifugal acceleration from 59.6g at the meniscus to 71.8g at the base (assuming a 22 mm sample height). The percentage difference in the centrifugal acceleration between the base and meniscus is 18.6%. In contrast to the Earth gravity simulation, the volume fraction of the suspension component decreases with time. However, the volume fraction remains uniform across the suspension component. The volume fraction decreases from ϕ = 0.03 initially to ϕ = 0.026 at t = 39 min, a reduction of approximately 13%. Indeed, the centrifugal acceleration decreases with height in the sample, so the terminal velocity of a given particle is higher near the base of the sample. This effect results in acceleration of the separation velocity in time, which can be seen from Fig. 4(b) by the spatial distance between each profile front at 0.5ϕ₀ increasing in time, notwithstanding the constant time period between each profile. The double-headed arrows in Fig. 4 illustrate the normalised distances between the profiles, with the normalised distance between subsequent profiles increasing in time under centrifugal acceleration [Fig. 4(b)], but remaining constant in Earth gravity conditions [Fig. 4(a)].

Due to the radially varying acceleration in the experimental set-up, particles in the suspension closer to the sediment bed have a higher sedimentation velocity than those above. As a result, particles are evacuated into the sediment bed faster 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 decreases. Therefore, particles are evacuated into the bed less efficiently as time goes on, so the rate of decrease of the volume fraction of the suspension decreases in time.

However, the volume fraction of the suspension remains uniform. This is because the particles entering the sediment bed make way for those settling behind them. Locally this will result in a reduction in volume fraction and a subsequent decrease in hindrance. From the formulation of the mixture viscosity in Eq. (7), the hindrance experienced by the particles immediately above is reduced, which, by Eq. (4), leads to an increase in the hindered settling velocity of the particles. Therefore, throughout the suspension component, the increase in acceleration towards the sediment bed is balanced by the reduction in the hindrance, resulting in a uniform volume fraction in the suspension component, despite its reducing value.

Note that the profile of the travelling front is more abrupt in the centrifugal force simulation. This is because there is less time for diffusion to act compared to sedimentation at Earth gravity. The Earth gravity conversions in the literature (see Sec. IV) do not account for this nor do they directly address the variation in the centrifugal acceleration across the sample, which may introduce a loss of accuracy into the derived Earth gravity separation velocities. Centrifugation has altered the sedimentation behaviour since the centrifugal acceleration varies across the suspension. As a result, linear conversion to Earth gravity using a single RCF value may be an unrealistic expectation. The conversion proposed in Sec. IV does not assume linearity and is therefore expected to be more reliable.

In Fig. 5, experimental and model normalised interface height (h/x_max) profiles are compared for two suspensions of initial volume fraction ϕ₀ = 0.03 and ϕ₀ = 0.05. For each suspension, profiles under both Earth gravity and centrifugal force (705 rpm) are shown. The experimental data were collected under centrifugal force and converted to Earth gravity conditions using the method proposed in Sec. IV. The model data were obtained in the appropriate condition, using the appropriate acceleration. Notice that under centrifugal force conditions, complete sedimentation is achieved in less than 1 h, whereas approximately 50 h is required in Earth gravity conditions.

Figures 5(a) and 5(b) show the profiles in centrifugal force conditions, for the ϕ₀ = 0.03 and ϕ₀ = 0.05 suspension, respectively. For both suspensions, the model and experimental profiles agree very well, suggesting that the model is capturing the behaviour of the experimental system correctly. Indeed, the percentage difference in the time required for complete sedimentation between the model and experiment is 3.6% for the ϕ₀ = 0.03 suspension, and 1.0% for the ϕ₀ = 0.05 suspension, with an overprediction in the separation rate by the model in both cases. This small overprediction may be due to uncertainty in the location of the base of the cell in the experimental centrifuge (prescribed in the model) or the experimental silica not being truly monodisperse. The effect of polydispersity would be more apparent in lower volume fractions due to the lower number of particles, which may explain why the overprediction is lower for the ϕ₀ = 0.05 suspension than the ϕ₀ = 0.03 suspension. For both suspensions, the profiles are curved due to the increase in centrifugal acceleration as the interface height decreases. The model correctly captures the curve in the experimental data, providing confidence in the centrifugal force extension to the model.

Having demonstrated the model’s ability to capture the unconverted experimental data in centrifugal force conditions, the proposed Earth gravity conversion can be analysed. Figures 5(c) and 5(d) show the experimental data converted to Earth gravity conditions using the proposed conversion. The model results are from the corresponding simulations in Earth gravity conditions, so no conversion is required. The unconverted model follows the converted experimental data very closely, and with similar small overpredictions as seen in the unconverted centrifugal force graphs. The percentage difference in the time required for complete sedimentation between the model and experiment is now 2.3% for the ϕ₀ = 0.03 suspension and less than 0.1% for the ϕ₀ = 0.05 suspension. Further, the curve in the experimental data profile has been almost completely removed, yielding a nominally constant separation velocity in Earth gravity conditions, as would be expected for low volume fraction, monodisperse suspensions. These observations suggest that the proposed Earth gravity conversion is able to convert the experimental data accurately. Therefore, the converted experimental data can be confidently used to further assess the accuracy of the model’s predictions.

In Sec. V D, it was shown that the proposed Earth gravity conversion can correctly recover Earth gravity settling data from experiments conducted in centrifugal force conditions. From the linear interface height profiles of the converted experimental data, separation velocities can be extracted to compare those predicted by the model in Earth gravity conditions.

Figure 6 shows separation velocities for suspensions with an initial volume fraction from ϕ₀ = 0.01 to ϕ₀ = 0.09. For each volume fraction studied experimentally, data from different experiments are shown explicitly, with three data points plotted, except ϕ₀ = 0.04, where two are plotted (due to experimental irregularities). The results of the proposed model (see Sec. III D) are shown (⁠$ϕmaxeff$ model). In addition, the results of a model with the hard-sphere assumption reintroduced, by replacing $ϕmaxeff$ with ϕ_max in Eq. (8), are shown (ϕ_max model) for comparison. The $ϕmaxeff$ model follows the experimental data very closely, showing that the individual settling velocities are well predicted by the proposed model. Further, this observation indicates that estimation of hindered settling effects by altering the fluid viscosity in line with a mixture model [Eq. (8)] closely approximates the fluid dynamics in the suspension, resulting in an appropriate reduction in the separation velocity on increasing the initial volume fraction. Moreover, the ϕ_max model deviates from the experimental data for the larger volume fraction suspensions studied, especially compared to the $ϕmaxeff$ model. Initial volume fractions higher than ϕ₀ = 0.09 were not investigated due to practical limitations on the initial volume fraction in the experimental suspensions. Larger deviations between the two models would be expected for higher volume fractions, where the proportion of the suspension accounted for by the effective maximum volume fraction increases. The improved fit to the experimental data of the proposed model compared to the equivalent hard-sphere model suggests that the elimination of the hard-sphere assumption has had a positive effect on the model. Hence, interparticle interactions are shown to have a significant effect on sedimentation behaviour of suspensions of colloidal particles.

As mentioned in Sec. II A, the silica nanoparticles used in the experiments were labeled 500 nm. Dynamic light scattering experiments were used to measure the hydrodynamic particle diameter, with a mean value of 480 nm (see Sec. II A). An approximate variation of ±10 nm in the mean value was observed across different measurements. To understand the effect of varying the particle size according to this tolerance on the sedimentation results, model simulations assuming three different particle sizes (470 nm, 480 nm and 490 nm) are shown in Fig. 7. The variation in particle size has some effect on the separation velocity for each initial volume fraction, due to Eq. (5). However, the variation in particle size only has a minimal effect on the reduction of separation velocity on increasing the initial volume fraction, which is evidenced by the model curves in Fig. 7 being almost parallel. This minimal effect is due to the ratio of the minimum separation distance to the particle size changing only slightly with a 20 nm variation in particle size. The separation velocity is sensitive to the particle size but does not change the conclusion regarding the effectiveness of the hindered settling function proposed, within the measured tolerance. From Fig. 6, it is determined that the particle size of 480 nm is the best fit to the experimental data, corroborating the results of the dynamic light scattering experiments and justifying the use of this particle size in the model.

Another factor which could affect a given experimental suspension’s separation velocity, when converted to Earth gravity conditions, is the centrifuge rotation rate. Figure 8 shows experimental separation velocities collected using two different rotation rates (705 rpm and 1000 rpm) and converted to Earth gravity conditions. For each volume fraction studied experimentally, three data points are plotted for each rotation rate, except ϕ₀ = 0.04 at 705 rpm and ϕ₀ = 0.03 at 1000 rpm, where two are plotted (due to experimental irregularities). For all volume fractions, the experimental separation velocities are similar for both rotation rates, suggesting that the choice of rotation rate does not significantly alter the converted Earth gravity separation velocity, which provides further support to the proposed Earth gravity conversion method. In practice, the rotation rate used should be sufficient to provide a centrifugal acceleration which dominates the Earth gravity acceleration, minimising slumping in the sample when in the centrifuge.

Combining Eq. (4) with Eq. (6) and taking the logarithm yields the equation for a straight line in the variables log(u) and log(1 − ϕ),

This line can be fitted to experimental data using linear regression to obtain the parameters u_T and p. For the 705 rpm experimental data, u_T = 0.4703 mm h⁻¹ and p = 4.70. Note that the value of p is close to the common literature value of 4.65 which is often used for experiments in laminar flow conditions.^12,14 The ability to recover this value provides further confidence in the Earth gravity conversion and experimental procedure.

Using the fitted value of u_T, the particle diameter, d can be extracted from the sedimentation data using Eq. (5), as a verification to the dynamic light scattering experiments for particle size. Only u_T and well-characterised parameters μ_c, ρ_d, ρ_c, and g are required. By substitution in Eq. (5), the particle size is determined to be 481 nm, which is remarkably close to the dynamic light scattering experiment value. Again, the ability to recover this value provides further confidence in the Earth gravity conversion and the choice of particle size used in the model (480 nm).

Further to the experiment-fitted line, the model curve in Earth gravity conditions is plotted in Fig. 8 for comparison. This model line lies almost entirely on top of the experiment-fitted line, which serves to provide further verification to the proposed model, and in particular the use of an effective maximum volume fraction to take into account interparticle interactions.

A one-dimensional sedimentation model for monodisperse, colloidal suspensions is proposed, which uses a suspension viscosity model to take into account the effect of hindered settling. The model avoids the hard-sphere assumption by using an effective maximum volume fraction, which is computed using DLVO theory. The proposed model is shown to fit the experimental data well, indicating that interparticle interactions are an important aspect in determining the sedimentation velocities of colloidal suspensions. The inclusion of interparticle interactions is a significant addition to sedimentation models in the literature and will be especially applicable at high volume fractions where the additional proportion of the suspension accounted for by the effective maximum volume fraction is larger.

The use of an adaptive spatial mesh in a finite difference solver has been shown to be an efficient way to numerically implement the sedimentation model, allowing the simulation to be numerically stable, whilst using a relatively low number of spatial mesh points.

Using the proposed model’s extension to centrifugal force conditions, the effects of centrifugal acceleration on a sedimenting suspension were investigated. It was found that the volume fraction of the suspension component reduced in time, whereas it remained constant at Earth gravity. Nevertheless, the volume fraction of the suspension remains uniform for all times, an effect which can be explained due to hindrance. The separation velocity was shown to increase in time under centrifugal force, due to the decrease in centrifugal acceleration with height in a sample in a centrifuge. On the basis of these observations, a method to convert centrifugal acceleration sedimentation data to Earth gravity conditions, which uses the space-averaged centrifugal acceleration experienced by the supernatant-suspension interface to account for the variation in centrifugal acceleration across the sample, was proposed. The conversion was shown to accurately recover Earth gravity sedimentation data and could be extremely useful in consumer product stability testing since experiments that usually take months can be accurately carried out in hours, with confidence in the derived Earth gravity results.

Further experimental data, particularly in Earth gravity conditions, could be used to provide additional validation for the conversion method. It would be particularly insightful to investigate larger volume fractions using a suitable experimental suspension, where the effect of using an effective maximum volume fraction is likely to be more pronounced. The proposed model could be extended to account for polydispersity in particle size following Abeynaike et al.,¹² which would increase the number of physical systems it is applicable to.

See supplementary material for the MATLAB code (implementing the sedimentation model) described in Sec. III E.

This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) Centre for Doctoral Training in Fluid Dynamics at the University of Leeds under Grant No. EP/L01615X/1. The authors would like to thank Vincenzo Guida for helpful discussion and Procter and Gamble (P&G) for supporting the project.