Multidimensional Separation by Magnetic Seeded Filtration: Theoretical Study

Abstract

Magnetic seeded filtration (MSF) is a multidimensional solid–liquid separation process capable of fractionating a multimaterial suspension based on particle size and surface properties. It relies on the selective hetero-agglomeration between nonmagnetic target and magnetic seed particles followed by a magnetic separation. Experimental investigations of multimaterial suspensions are challenging and limited. Therefore, a Monte Carlo model for the simulation of hetero-agglomeration processes is developed, validated, and compared to a discrete population balance model. The numerical investigation of both charge-based and hydrophobicity-based separation in an 11-material system, using synthetic agglomeration kernels based on real-world observations, yields results consistent with prior experimental studies and expectations: Although a multidimensional separation is indeed possible, unwanted hetero-agglomeration between target particles results in a reduced selectivity. This effect is more pronounced when separation is based on a dissimilarity rather than a similarity in the separation criterion and emphasizes the advantages of hydrophobicity-based systems. For the first time, 2D grade efficiency functions $T(\phi ,d)$ are presented for MSF. However, it is shown that these functions strongly depend on the initial state of the suspension, which casts doubt on their general definition for agglomeration-based processes and underlines the importance of a simulation tool like the developed MC model.

1. Introduction

The start-up and scale-up of industrial processes involving solid particles has always been a major challenge in the field of process engineering [1,2]. Unlike fluids, solid particles are characterized by distributed properties that make predicting their behavior challenging. The traditional approach of assuming that particles are only distributed in size, while, e.g., shape and surface properties are constant, proves to be both inadequate and even neglectful. Thus, a multidimensional description of particle property distributions is required and is currently a focus of research [3,4]. Naturally, this concept permeates all areas of particle technology, including solid–liquid separation: For one, separation processes are required to operate on multiple particle properties simultaneously [5,6]. Additionally, the description of the separation result—traditionally given via the size-dependent grade efficiency $T\left(d\right)$—must be further developed to account for the multidimensional nature of particle systems [7,8].

A previous study [9] proposed magnetic seeded filtration (MSF) as a multidimensional process that separates both by particle size and surface properties. Here, magnetic seed particles are added to a suspension containing the nonmagnetic target particles, and a selective hetero-agglomeration is induced. This is the process-defining step and is dependent on both the particle surface and size. These hetero-agglomerates are then removed from the suspension by a simple magnetic separation. During the previous experimental investigations [9], a polydisperse suspension with two nonmagnetic materials of different zeta potential was separated into two fractions of more or less defined particle size and material. Thus, in general, MSF is able to achieve a selective and multidimensional separation. However, the expected selectivity based on component-specific separation experiments was not achieved, which was shown to be due to an undesired hetero-agglomeration between the target particles. Another preceding experimental study [10] investigated the separation based on the hydrophobic properties of the surface: By addition of hydrophobic magnetic material, hydrophobic microplastic particles were selectively separated from hydrophilic cellulose particles. The loss of selectivity observed in the charge-based system was not found in the hydrophobicity-based separation.

Experimental work investigating a multidimensional separation currently faces several downsides and limitations. As discussed in [9], analytical methods are intricate, even in suspensions with only two nonmagnetic materials, necessitating the combination of multiple offline analyses to obtain material-specific information. Thus far, separation was quantified continuously with respect to size and discretely, i.e., material-wise with respect to surface charge. A complete multidimensional description of the separation requires a continuous 2D grade efficiency function $T(\phi ,d)$, where $\phi$ is a surface property of interest. This necessitates either the possession and analytical capability of a material continuously distributed with respect to the surface property or a mixture of indefinitely more discrete materials. Both solutions are presently unfeasible given the current experimental capabilities. Furthermore, due to the reliance on offline analyses before and after separation, the statements that can be made are confined to integral quantities, such as the material-specific separation efficiency. However, the pivotal hetero-agglomeration processes, which define the overall outcome, remain concealed.

Simulations play a crucial role in elucidating processes and phenomena that are challenging to measure directly. Population balance equations (PBEs) are the standard way to describe the kinetics of particle agglomeration processes [11,12]. The underlying integro-differential equations can be solved numerically by discretizing the agglomerate property space, resulting in so-called discrete population balance equations (dPBEs). In previous work, a numerical dPBE solver for systems of up to three materials was developed and applied to MSF [13,14]. However, computation times of dPBE models scale exponentially with the number of materials, making the desired transition to a more continuous multidimensional separation challenging. Stochastic or Monte Carlo (MC) models take an entirely different approach by simulating the evolution of a discrete particle population without the need for a mass balance and differential equations [11,12]. They are easy to implement and can handle an arbitrary number of materials, although they are considered computationally expensive [12].

The present work aims to address the research gaps identified above by pursuing the following objectives. First, an MC model is developed and validated to simulate agglomeration processes in arbitrary multimaterial systems. Coupled with a magnetic separation criterion, the continuous 2D grade efficiency curve $T(\phi ,d)$ for MSF is estimated and visualized for the first time. The performance and scaling behavior of the MC model are compared to the existing dPBE model and discussed in the context of multidimensional separation. Subsequently, the MC model is applied in two theoretical case studies: One based on charge and one based on hydrophobicity. Both cases are conceptually different as the separation is based either on dissimilarity (charge) or similarity (hydrophobicity) in the separation criterion. The use of synthetic and well-defined kinetic models enables an in-depth investigation into the effect of intentional changes in multidimensional particle property distributions on the 2D grade efficiency function and selectivity.

2. Models

2.1. Discrete PBE

The fundamental kinetics of the agglomeration process between two agglomerates i and j, resulting in agglomerate $ij$, are given by

$$\frac{\mathrm{d}{n}_{ij}}{\mathrm{d}t}={\beta }_{ij}{\alpha }_{ij}{N}_i{N}_j,$$

(1)

where n is the number density of the respective class. The kinetic rates or so-called kernels $\beta$ and $\alpha$ contain the physics of the agglomeration process. $\beta$ describes the collision frequency, i.e., the number of collisions per unit time, and $\alpha$ the collision efficiency, i.e., the probability that a collision leads to an agglomeration. In a particle system with M discrete materials, an agglomerate $ij$ is completely defined by its M partial volumes ${\mathbf{x}}_{ij}=\{x_1,x_2,\dots ,x_M\}$, if agglomerate shape is neglected and its porosity is assumed to be 0. Equation (1) alone is not sufficient to describe an entire particle system with varying agglomerate sizes and compositions (i.e., varying ${\mathbf{x}}_{ij}$). The population balance equations (PBEs) are obtained by applying Equation (1) to all possible agglomerate combinations and balancing the respective production and loss terms of n. In this context, $\mathbf{x}$ are called the internal coordinates of the PBE and due to mass conservation, the partial volumes are additive, i.e., ${\mathbf{x}}_{ij}={\mathbf{x}}_i+{\mathbf{x}}_j$. The PBE is as follows:

$$\frac{\partial n(\mathbf{x},t)}{\partial t}=\frac{1}{2}{\int }_0^{\mathbf{x}}k(\mathbf{x}-{\mathbf{x}}^{\prime },{\mathbf{x}}^{\prime })n(\mathbf{x}-{\mathbf{x}}^{\prime },t)n({\mathbf{x}}^{\prime },t)\mathrm{d}{\mathbf{x}}^{\prime }-{\int }_0^{\infty }k(\mathbf{x},{\mathbf{x}}^{\prime })n(\mathbf{x},t)n({\mathbf{x}}^{\prime },t)\mathrm{d}{\mathbf{x}}^{\prime },$$

(2)

where $k(i,j)={\beta }_{ij}{\alpha }_{ij}$. Equation (2) is an integro-differential equation and can only be solved analytically for certain kernels and initial conditions. The numerical solution requires a discretization along its internal coordinates $\mathbf{x}$, which transforms Equation (2) into a set of ordinary differential equations (ODE) that can be solved numerically. The simplest discretization scheme is the so-called uniform grid, which scales the discrete partial volume of each material $m=\{1,2,\dots ,M\}$ linearly according to

$$x_m\left(z\right)=z{x}_m\left(1\right),$$

(3)

where z is the discrete index $z=\{1,2,\dots ,Z\}$, Z is the number of grid points per coordinate, and $x_m\left(1\right)$ is the smallest partial volume of material m. Since each partial volume is discretized individually, the number density for all possible agglomerate classes n is now a $Z^M$ matrix. The number of materials M is called the dimension of the discrete PBE (dPBE). Modeling the selective separation between two nonmagnetic material systems by addition of magnetic particles therefore requires a three-dimensional ($M=3$) dPBE.

In addition to the numerical solution scheme, the kernel values $\beta$ and $\alpha$ are the bottleneck for accurate results. Various mechanistic and (semi-) empirical models for calculating the kernels from process parameters exist and are summarized in [15,16,17]. However, as is discussed in detail in a previous publication [14], the necessary process and material data are often either not sufficiently accurate or not available at all, making purely predictive calculations impossible. As proposed in [14], parametrization of the kernel functions and optimization of the unknown parameters to experimental data allows the training of data-driven models and results in so-called hybrid models. In this work, the orthokinetic model [17,18]

$${\beta }_{ij}=\overline{G}{(r_i+r_j)}^3$$

(4)

is used for the collision frequency, with r being the agglomerate radius and $\overline{G}$ being the mean shear rate that is calibrated with experiments. Note that multiscale modeling approaches like coupling the dPBE with CFD-DEM models might provide additional information on this parameter. The previously developed collision case model [13,14] is used to predict ${\alpha }_{ij}$. It requires the (known) partial volumes ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ and the material-specific agglomeration efficiencies $\mathit{\alpha }$. $\mathit{\alpha }$ is a symmetric $M\times M$ matrix (${\alpha }_{ij}={\alpha }_{ji}$) containing the agglomeration efficiency for each material combination. In essence, ${\alpha }_{ij}$ summarizes the particle–particle interactions (van der Waals, electrostatic, hydrophobic) and can be estimated from experiments. It should be noted that agglomeration occurs outside of the magnetic field; therefore, the magnetic interaction is not considered in the dPBE model.

2.2. Monte Carlo

Instead of tracking and balancing all possible agglomeration events simultaneously (see Equation (2)), Monte Carlo (MC) models follow a discrete set of individual agglomerates and represent a completely different approach to computing the temporal evolution of the number density distribution n. A population matrix X ($M\times N$) is used for the calculations, containing the M partial volumes of N individual agglomerates. At first, X is initialized with $N_0$ agglomerates, which are chosen stochastically according to the initial number density distribution $n(\mathbf{x},0)$. Setting the initial number of agglomerates $N_0$ equal to its true value in a real-world system is infeasible, and, instead, $N_0$ is chosen as a representative sample. This implies that the MC model does not simulate the entire volume of the system, but only a fraction of it, the so-called control volume

$$V_c=\frac{\sum X}{c_v},$$

(5)

where $\sum X$ is the total solid volume of the sample and $c_v$ is the overall volume concentration of the real-world suspension. Since every possible agglomeration event requires the collision of two unique agglomerates, the collision frequency array B ($\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.N(N-1)\times 1$) can be calculated for any given population matrix via Equation (4), while the probability of each collision event is given by $P=B/\sum B$. The probability that a particular collision will result in an agglomeration is described by the collision efficiency $\alpha$, which is calculated according to the collision case model described above and stored in a collision efficiency array A ($\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.N(N-1)\times 1$).

The MC algorithm selects a specific collision l stochastically according to its probability $P\left(l\right)$, generates a random number $\mathrm{rnd}\in [0,1]$, and checks via $\mathrm{rnd}<A\left(l\right)$ whether an agglomeration occurs or not. If so, both collision partners i and j are removed from the population matrix while a new agglomerate $ij$ with ${\mathbf{x}}_{ij}={\mathbf{x}}_i+{\mathbf{x}}_j$ is added. This reduces the number of individual agglomerates by 1 and requires both B and A to be updated. On average, this collision event takes

$$\Delta t=\frac{2V_c}{\overline{B}{N}^2},$$

(6)

where $\overline{B}$ is the mean collision frequency of all possible collision events. $\Delta t$ is called the inter-event time and advances the overall simulation time. All previous steps are repeated until the simulation time reaches the predefined agglomeration time $t>t_{\mathrm{A}}$. Since each agglomeration event reduces the total number of agglomerates by 1, choosing $N_0$ to be too small with respect to $t_{\mathrm{A}}$ results in the formation of a single large agglomerate and a premature termination of the simulation. To counteract this, when the current number of agglomerates N reaches half of the initial number $N_0$, i.e., $N\le N_0/2$, the control volume $V_c$ is doubled and all individual agglomerates are duplicated [19].

2.3. Magnetic Separation

Both the dPBE and MC model yield the time-dependent number density distribution $n(\mathbf{x},t)$. However, from a process engineering standpoint, the overall separation efficiency and selectivity of the MSF process are more relevant. To close this gap, a magnetic separation model derived in a previous publication [14] has to be employed. The probability that a given agglomerate i is separated is estimated by

$${\tau }_i=\left\{\begin{array}{cc}1\hfill & r_i^2{\tilde{x}}_{\mathrm{mag},i}\ge C_{\mathrm{mag}}\hfill \\ 0\hfill & r_i^2{\tilde{x}}_{\mathrm{mag},i}<C_{\mathrm{mag}}\hfill \end{array}\right.,$$

(7)

where $r_i$ is the agglomerate radius and the magnetic volume fraction ${\tilde{x}}_{\mathrm{mag},i}$ is calculated according to

$${\tilde{x}}_m=\frac{x_m}{{\sum }_{j=1}^M{x}_j}.$$

(8)

Equation (7) represents a perfectly sharp separation and $C_{\mathrm{mag}}$ defines the properties of an agglomerate that is just barely collected at the separation matrix and therefore summarizes all relevant parameters like, e.g., the magnetization of the magnetic material. The material-specific separation efficiency $A_m$ describes the fraction of overall separated volume of material m with respect to its initial volume:

$$A_m=\frac{{\int }_0^{\infty }n(\mathbf{x},t_{\mathrm{A}})\tau \left(\mathbf{x}\right)x_m\mathrm{d}\mathbf{x}}{{\int }_0^{\infty }n(\mathbf{x},0)x_m\mathrm{d}\mathbf{x}}$$

(9)

As derived in the previous study [9], the material-specific selectivity $S_m$ is defined by

$$S_m=\frac{A_m{\tilde{v}}_{m,0}}{{\sum }_{j=1}^M{A}_j{\tilde{v}}_{j,0}}$$

(10)

with ${\tilde{v}}_{m,0}$ being the initial volume fraction of material m given by

$${\tilde{v}}_{m,0}=\frac{{\int }_0^{\infty }n(\mathbf{x},0)x_m\mathrm{d}\mathbf{x}}{{\sum }_{j=1}^M{\int }_0^{\infty }n(\mathbf{x},0)x_j\mathrm{d}\mathbf{x}}$$

(11)

The material-specific separation efficiency $A_m$ and selectivity $S_m$ can be obtained from both the dPBE and MC results. However, they represent only a single dependence (of material m) and do not give any information on the size-specific separation. As mentioned in the introduction, the goal of a multidimensional separation is to find a continuous 2D grade efficiency function $T(\phi ,x)$ with respect to surface property $\phi$ and agglomerate size (volume) x. This is challenging for two reasons: First, due to the discrete nature of materials in the real world, $\phi$ will never be truly continuous. At best, the number of materials m can be increased ($m\to \infty$) and the material-specific grade efficiencies $T_m\left(x\right)$ calculated. Second, and most importantly, the dPBE model is not able to provide the required data because it does not retain individual agglomerate information and because partial volumes are constantly flowing between size classes. Vividly, when looking at any agglomerate class i at time $t_{\mathrm{A}}$, only the number density and the partial volumes in this class are known, but not which individual, initially present agglomerates have accumulated. Therefore, a size- and material-specific comparison before and after separation is not possible. MC simulations are different, as the history and evolution of individual particles over the course of an agglomeration process is traceable. Therefore, applying Equation (7) to all agglomerates at time $t_{\mathrm{A}}$ ($X\left(t_{\mathrm{A}}\right)$) directly gives information about the separation of initially present agglomerates at time $t=0$ ($X\left(0\right)$).

An example that illustrates this discussion is shown in Figure 1: Starting with an initial state of 10 discrete particles of $M=4$ different materials, both the MC and dPBE model are able to calculate the state of agglomeration after $t_{\mathrm{A}}$. In this example, three agglomerates A₁, A₂, and A₃ are generated, and applying Equation (7) yields that A₁ and A₃ are separated (indicated by lowered saturation), while A₂ remains. Since the dPBE only has information on the partial volumes $\mathbf{x}$ of these agglomerates, only the overall, material-specific separation efficiency $A_m$ and selectivity $S_m$ can be calculated, while the connection to the initially present, pure particles is lost. MC simulations, on the other hand, are able to keep track of which agglomerates contain which initial particles, and the separation of A₁ and A₃ is directly linked to the initial state.

Sorting the initial agglomerate population into defined classes of overall surface properties and size ($\widehat{\phi }$,$\widehat{d}$), where d is the diameter of the volume-equivalent sphere, allows the calculation of the volume-weighted, 2D grade efficiency function

$$T(\widehat{\phi },\widehat{d})=\frac{V_{\mathrm{sep}}(\widehat{\phi },\widehat{d})}{V_0(\widehat{\phi },\widehat{d})}$$

(12)

by simply dividing the separated $V_{\mathrm{sep}}$ by the initially present volume $V_0$ in each class. This assumes that at $t=0$, only single-material particles or agglomerates are present to avoid the need for defining an “average” surface property $\overline{\phi }$ for multimaterial agglomerates.

3. Validation and Grid Study

3.1. Validation

Both the dPBE and MC model discussed in Section 2 are validated by comparison with an analytical solution in the two-component system ($M=2$). Instead of comparing the full number density distribution n, visualizing the moments

$${\mu }_{ij}\left(t\right)={\int }_0^{\infty }{\int }_0^{\infty }{x}_1^i{x}_2^{j}n(\mathbf{x},t)\mathrm{d}{x}_1\mathrm{d}{x}_2$$

(13)

of the distribution allows for a straightforward assessment of the accuracy of the calculations. In particular, the zeroth moment ${\mu }_{00}$ represents the total number of agglomerates, while the first moments ${\mu }_{10}$ and ${\mu }_{01}$ represent the total volume of each material over time, thus accounting for mass conservation. For monodisperse initial conditions, the collision frequency following the sum kernel ${\beta }_{ij}={\beta }_0(\sum {\mathbf{x}}_i+\sum {\mathbf{x}}_j)$ and $\mathit{\alpha }=\mathbf{1}$; the moment expressions for the solution derived by Scott [20] are summarized in [21,22] and given in Equations (14)–(16).

$$\begin{array}{c}{\mu }_{00}\left(t\right)=n_{0}exp(-{\beta }_0{n}_{0}t)\end{array}$$

(14)

$$\begin{array}{c}{\mu }_{10}\left(t\right)=v_{1,0}|{\mu }_{01}\left(t\right)=v_{2,0}\end{array}$$

(15)

$$\begin{array}{c}{\mu }_{11}\left(t\right)=\frac{v_{1,0}{v}_{2,0}(2-\varphi )\varphi }{n_0{(1-\varphi )}^2}\\ n_0={\int }_0^{\infty }n(\mathbf{x},0)\mathrm{d}\mathbf{x}|v_{m,0}={\int }_0^{\infty }n(\mathbf{x},0)x_m\mathrm{d}\mathbf{x}|\varphi =1-exp(-{\beta }_0{n}_{0}t)\end{array}$$

(16)

Figure 2a–c show the relative zeroth, first, and first cross moments over time for the analytical solution, the MC and dPBE model. All calculations used ${\beta }_0={10}^{-16}{\mathrm{m}}^3{\mathrm{s}}^{-1}$, $v_{1,0}=v_{2,0}={10}^{-5}$ and $d_{1,0}=d_{2,0}=1$µm and were run for $t_{\mathrm{A}}=800$ s. The MC simulations were run with $N_0=500$ and repeated $N_{\mathrm{MC}}=100$ times. Mean moment values are plotted and the error bars indicate a single standard deviation. The dPBE model follows a uniform discretization and its results are shown for two different grid parameters, $Z=12$ and $Z=120$. In general, all models closely follow the analytical solution up to $t\approx 200$ s, indicating no major errors in the implementation. For increasing agglomeration times, the results of the dPBE model differ strongly from the analytical solution, while increasing Z shifts the differences to larger values of t. Vividly, the agglomeration slows down and approaches a steady state. The MC model, on the other hand, closely follows the analytical solution, and the analytical moments are well within the confidence interval. Figure 2d visualizes the cumulative, volume-weighted agglomerate size distribution $Q_3$ at time $t_{\mathrm{A}}$ and explains the apparent problems of the dPBE model: In both cases, the grid parameter Z is insufficient to represent the full range of agglomerate sizes. When agglomerates accumulate in larger classes, they are no longer able to participate in agglomeration events and the entire process slows down. This is especially problematic for longer agglomeration times. It should be noted that other discretization schemes, such as the geometric grid, exist and have been implemented in previous studies [14], which make the calculation of broadly distributed particle systems manageable. However, they require more sophisticated numerical schemes like the cell average technique [21] and introduce discretization errors due to underlying assumptions.

Figure 2 highlights the general advantage of MC models, which are inherently grid-free: they neither suffer from numerical errors due to grid limitations nor require elaborate discretization schemes. Additionally, increasing the number of material systems M does not significantly change the MC algorithm, while it requires a different scheme and implementation for the dPBE. Coupled with the fact that only MC models provide the possibility to derive a 2D grade efficiency function $T(\widehat{\phi },\widehat{d})$, the case studies in Section 4 are only performed with the MC model.

3.2. Grid Study

Another interesting point of comparison between dPBE and MC is the required simulation time. Figure 3 shows the relative calculation times (normalized to the overall minimum value) for MC simulations with varying number of initial agglomerates $N_0=\{500,1000,3000\}$ and for dPBE simulations with varying grid parameters $Z=\{8,12,15\}$. The calculations use the same initial and boundary conditions as for the validation study in Section 3.1. Additionally, the one-component and three-component system were evaluated. All calculations were performed on identical hardware.

The calculation times of MC simulations are strongly dependent on the number of simulated agglomerates, i.e., on the sample size N. First, the MC algorithm itself requires more iterations due to a shorter inter-event time (see Equation (6)). However, the main reason for the large increase in computation time with N lies in the calculation of the collision frequency and collision efficiency arrays B and A. As mentioned in Section 2.2, they store all unique collision events and require $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.N(N-1)$ calculations each ($\mathcal{O}(N^2$)). In addition, each successful agglomeration event requires their recalculation or adjustment due to newly added and lost agglomerates. Xu et al. [23] summarizes several techniques to speed up these calculations, such as smart bookkeeping [24] or majorant kernels [25,26], which reduce the order to $\mathcal{O}\left(N\right)$. Calculations can further be parallelized to achieve a more efficient simulation. Since the present study aims at investigating the general capabilities of MC models for simulating multidimensional MSF experiments, and performance is not a top priority, neither efficient numerical schemes nor optimized code are implemented at this stage. The MC results in Figure 3 should, therefore, be interpreted as a worst-case estimate for the scaling behavior of MC models. Increasing the number of materials M has only a marginal effect on the overall simulation time. As the population matrix X increases in size, the memory requirements and internal calculations both for the collision case model and newly formed agglomerates slightly increase. However, neither the dimension of B nor A is affected and their calculation accounts for the majority of required time.

The calculation time for the dPBE model increases with increasing grid size Z. During the numerical solution of the ODE, each evaluation of the differential equation requires calculations for every possible combination of all discrete classes. Since the number of discrete classes scales according to $Z^M$, the solution of the dPBE is of order $\mathcal{O}\left(Z^{2M}\right)$. This exemplifies the underlying problem of the dPBE model in multicomponent systems, which is also clearly seen in Figure 3: Although the grid parameter is relevant, its role is strongly outweighed by the number of materials M. Calculations for more than three components quickly become infeasible for realistic values of Z.

4. Case Studies

4.1. Study Structure

Multidimensional separation with MSF is studied theoretically for two different cases of varying surface property $\phi$. Simple model functions for the collision efficiency matrix ${\mathit{\alpha }}_m$ are used which capture the underlying affinities discussed in the previous publication [9]. A “charge-based” separation is modeled for case 1 by

$${\alpha }_{ij}=\left\{\begin{array}{cc}0\hfill & {\zeta }_i{\zeta }_j\ge C_1\hfill \\ \frac{{\zeta }_i{\zeta }_j-C_1}{C_2}\hfill & C_1+C_2<{\zeta }_i{\zeta }_j<C_1\hfill \\ 1\hfill & {\zeta }_i{\zeta }_j\le C_1+C_2\hfill \end{array}\right.$$

(17)

which results in an attraction (${\alpha }_{ij}>0$) for oppositely charged particles that increases with the magnitude of the zeta potential $\zeta$. This reflects that the separation is based on a dissimilarity in the separation criterion. Particle pairs exceeding the critical zeta potential product (${\zeta }_i{\zeta }_j\le C_1+C_2$) are fully destabilized (${\alpha }_{ij}=1$). Equally charged particles can only agglomerate up to a certain charge (${\zeta }_i{\zeta }_j<C_1$) and are stable otherwise. The calculations are performed for 10 nonmagnetic materials with zeta potential values linearly spaced between $-30$ mV and $+30$ mV and one magnetic material with a zeta potential of $-30$ mV ($M=11$). Choosing the model parameters as $C_1=25{\mathrm{mV}}^2$ and $C_2=-900{\mathrm{mV}}^2$ results in the collision efficiency matrix shown in Figure 4a.

A “hydrophobicity-based” separation is modeled as case 2, according to

$${\alpha }_{ij}={\mathrm{RH}}_i{\mathrm{RH}}_j,$$

(18)

where $\mathrm{RH}\in \{0,1\}$ is the relative hydrophobicity of the surface, further discussed below. The underlying affinity arises from the fact that strong and long-ranged hydrophobic interactions have been shown to exist between two similar, hydrophobic particle surfaces [27,28], whereas dissimilar surfaces, i.e., one hydrophilic and one hydrophobic surface, have been shown to exhibit either no long-ranged or no attraction at all [28,29]. Separation is, therefore, based on a similarity in the separation criterion. A previously published study on MSF [10] supports this hypothesis, as a selective separation between hydrophilic cellulose and hydrophobic microplastic particles was performed by using a hydrophobic magnetic material. It should be noted that the discussion about the hydrophobic particle–particle interactions is still ongoing, and that no overall theory has been proposed to explain all of the observed effects. The bridging force of nanobubbles on the particle surface [30,31] is currently the most promising explanation. There are several measures of surface hydrophobicity, such as contact angle or surface energy, but their experimental determination is often challenging [32]. Additionally, no general, mechanistic model for the hydrophobic interactions based on these properties has been proposed as of yet. Therefore, Equation (18) is sufficient for the purpose of this study, as it captures the general dependencies. The relative hydrophobicity $\mathrm{RH}$ is a practical mathematical simplification that defines the minimum hydrophobicity of two interacting surfaces to achieve a complete destabilization ${\alpha }_{ij}=1$. Analogous to case 1, the calculations are performed for 10 nonmagnetic materials, five of which have an increasing degree of hydrophobicity ($\mathrm{RH}>0$). The relative hydrophobicity of the remaining five materials is set to $\mathrm{RH}=0$, although, conceptually, they can differ by being more or less hydrophilic without affecting their agglomeration behavior. For a selective separation, the magnetic material is strongly hydrophobic ($\mathrm{RH}=1$). The collision efficiency matrix is visualized in Figure 4b.

The collision frequency was calculated according to Equation (4) with a mean shear rate of $\overline{G}=60{\mathrm{s}}^{-1}$, and the agglomeration time was $t_{\mathrm{A}}=120\mathrm{s}$ in both cases. The parameter of the magnetic separation model in Equation (7) was set to $C_{\mathrm{mag}}=5\times {10}^{-13}{\mathrm{m}}^2$, indicating that an agglomerate with diameter $2$µm requires at least a magnetic volume fraction of ${\tilde{x}}_{\mathrm{mag}}=0.5$ to be separated. It is assumed that at $t=0$ no hetero-agglomeration has taken place; thus, the population matrix X is initialized with pure, i.e., single-material particles. A monodisperse magnetic material with $d=2$µm was used. The volume density distribution $q_3\left(d\right)$ for each of the 10 nonmagnetic materials followed a ${\log }_{10}$-normal distribution with mean value $\mu =2$µm and standard deviation $\sigma =0.15$. During initialization, $q_3\left(d\right)$ is transformed into a number density distribution $q_0\left(d\right)$, from which samples are drawn and added to the population matrix X. To increase numerical stability, only particles in the size range $1 \mathrm{\mu m}<d<4 \mathrm{\mu m}$ are allowed. All 11 materials had identical initial number concentrations of $n_{0,m}=2.39\times {10}^{13}{\mathrm{m}}^{-3}$ (corresponding to a volume concentration of $c_v={10}^{-4}$). Thus, the number of sampled particles was identical for each material at $\raisebox{1ex}{$N_0$}\!\left/ \!\raisebox{-1ex}{$M$}\right.$. A total number of sampled particles of $N_0=1000$ was used in the simulations, which is expected to give reliable results based on the validation in Section 3.1. To reduce statistical noise and increase the significance of the observed trends, the MC simulations were run $N_{\mathrm{MC}}=1000$ times and the results combined. Effectively, 1000 small, individual agglomeration volumes with 1000 particles each were simulated, which cannot interact with each other during agglomeration but are combined after $t_{\mathrm{A}}$, resulting in a total number of initial particles of ${10}^6$. This procedure massively reduces the computational time due to the properties of the MC algorithm discussed in Section 3.2. The resulting initial conditions for both cases are visualized in Figure 5.

4.2. Results Case 1: Charge-Based Separation

The results of the charge-based case are presented in Figure 6a–d. Figure 6a shows the volume-weighted histogram for the partial volumes of all 10 nonmagnetic materials after agglomeration at $t=t_{\mathrm{A}}$. Strongly charged materials tend to accumulate in larger agglomerates, while the sign of the charge has little effect. Materials closer to the isoelectric point ($\zeta =0$) agglomerate to a limited extent and remain predominantly in the initial size range. Performing the magnetic separation and tracing the separation back to the initially present particles yields the 2D grade efficiency visualized in Figure 6b. All investigated classes were separated to some degree. An increase in zeta potential (regardless of sign) or particle diameter results in improved separation. This trend is further illustrated in Figure 6c, where the overall separation efficiency $A_m$ for each material is plotted for four discrete agglomeration times. Naturally, the separation efficiency increases with agglomeration time for all materials. Again, the increase in separation efficiency with increasing zeta potential (regardless of sign) is clearly visible, but for any given magnitude of zeta potential $\left|\zeta \right|$, negatively charged particles are separated slightly less than their positively charged counterparts. These trends are also reflected in Figure 6d, where the material-specific selectivity $S_m$ is shown for the same time steps. Overall, relatively low selectivity values are obtained. Positively charged materials are separated most selectively, followed by negatively charged materials, while materials closer to the isoelectric point show the lowest selectivity values. The impact of process time is profound: For the shortest agglomeration time $t=12$ s, the separation is more selective towards positively charged materials with larger differences between materials. As process time increases, the selectivities level out and approach a value of $S_m=0.1$, which represents the complete loss of selectivity loss, where all materials are separated equally.

In essence, Figure 6b illustrates the feasibility of multidimensional separation, i.e., classification by size and surface charge, with MSF. However, the selectivity in a charge-based system is limited, as all materials undergo some degree of separation. This agrees well with the previous experimental study [9], where hetero-agglomeration phenomena between nonmagnetic target particles were shown to result in the formation of multimaterial agglomerates that are separated. Specifically, negatively charged nonmagnetic particles, generally not attracted to the magnetic particles (see Figure 4a), agglomerate to positively charged nonmagnetic particles. These newly formed agglomerates, in turn, form connections to the magnetic material and lead to separation. The fact that this is a two-step process, with only the second step leading to separation, explains the following two observations: First, positively charged particles are separated more effectively than their negatively charged counterparts because they are directly attracted to the magnetic particles. Second, at low agglomeration times, the initial and direct agglomeration processes between primary particles dominate, resulting in increased selectivity toward positively charged particles. With increasing agglomeration time, more multimaterial agglomerates are formed and selectivity is reduced. It follows that the separation of any material m depends on the concentration and charge of all other materials, i.e., on the initial state of the suspension $n(\mathbf{x},0)$.

To prove this, three additional simulations, C1.1, C1.2, and C1.3, were run with different initial state but otherwise identical parameters to case 1 (C1). C1.1 contained only magnetic and strongly negatively charged particles with ${\zeta }_1=-30\mathrm{mV}$ ($M=2$). C1.2 and C1.3 were simulations in the three-component system ($M=3$): C1.2 contained two nonmagnetic materials with ${\zeta }_1=-30\mathrm{mV}$ and ${\zeta }_2=+10\mathrm{mV}$, and C1.3 contained two nonmagnetic materials with ${\zeta }_1=-30\mathrm{mV}$ and ${\zeta }_2=+30\mathrm{mV}$. The resulting material-specific separation efficiencies $A_m$ are shown in Figure 7 and compared with the results of C1 from Figure 6c. As expected, no separation of nonmagnetic material was achieved in case C1.1, since the collision efficiency between both materials is ${\alpha }_{ij}=0$ (see Figure 4a). Comparing C1.2 and C1.3 with C1 shows that the positively charged nonmagnetic materials (${\zeta }_2=+10\mathrm{mV}$ and ${\zeta }_2=+30\mathrm{mV}$) are separated similarly to the full 11-component system, with only slightly increased separation efficiency. However, the negatively charged nonmagnetic material (${\zeta }_1=-30\mathrm{mV}$) is significantly affected: When paired only with the strongly positively charged material in C1.3, $A_m$ decreases from $A_m\approx 76\%$ to $A_m\approx 67\%$. When paired with the weakly positively charged material in C1.2, the decrease in separation efficiency is more drastic, reaching $A_m\approx 27\%$. However, the comparison with C1.1 shows that the presence of some form of positively charged material is required to achieve any separation at all. This reinforces the discussion above that negatively charged particles rely on agglomeration with other nonmagnetic materials to be separated. Furthermore, Figure 7 shows that the initial state of the suspension significantly affects the separation.

Regarding the size dependence, larger particles simply agglomerate more often due to their increased collision frequency (see Equation (4)). They are therefore generally more prone to hetero-agglomeration and their probability of coming into contact with magnetic particles is increased. This was also shown experimentally in the previous, experiment study [9], although it should be noted that the size dependence of the agglomeration process is multifaceted and by no means trivial: Large and/or highly unequal particle sizes lead to reduced collision efficiency $\alpha$ and affect magnetic separation. Both effects are not included in the current model, as they need to be calibrated to experimental results. However, with all these partially opposing effects, it is most likely that there is always some form of size dependence in any given parameter range that can be used to achieve multidimensional separation.

4.3. Results Case 2: Hydrophobicity-Based Separation

The results for case 2 are shown in Figure 8a–d with corresponding graph layouts to case 1. The volume-weighted histogram (Figure 8a) after agglomeration $t=t_{\mathrm{A}}$ shows that only hydrophobic particles ($\mathrm{RH}>0$) are able to agglomerate and leave the initial particle size range, while increasing relative hydrophobicity leads to stronger agglomeration. Accordingly, the 2D grade efficiency shown in Figure 8b indicates that hydrophilic particles ($\mathrm{RH}=0$) are not separated at all. The grade efficiency increases for more hydrophobic and larger particles, indicating a selective, multidimensional separation. The trend of increasing separation with increasing hydrophobicity is also visible in the material-specific separation efficiency in Figure 8c. Naturally, increasing the agglomeration time leads to higher $A_m$. The selectivity trends in Figure 8d are remarkable: For short agglomeration times $t=12$ s, the separation is more selective between strongly and weakly hydrophobic particles. With increasing agglomeration time, the differences decrease and the selectivities approach $S_m=0.2$, which represents the case where all hydrophobic materials are separated equally. Again, hydrophilic particles ($\mathrm{RH}=0$) are not separated at all, resulting in $A_m=0$ and $S_m=0$, and, thus, a perfect selectivity between hydrophobic and hydrophilic materials.

Figure 8 highlights the advantages of an agglomeration mechanism that is based on dissimilarity in the separation criterion: In the investigated model case, MSF is able to achieve a perfect separation between hydrophilic and hydrophobic materials, since hydrophilic particles are attracted neither to magnetic nor to other hydrophobic nonmagnetic particles. This agrees well with a previously published experimental study, where hydrophobic microplastic particles were selectively separated from hydrophilic cellulose [10]. Similar to case 1, the material-specific selectivity of the MSF process is most pronounced for low agglomeration times. Again, this is due to hetero-agglomeration between hydrophobic nonmagnetic particles and a resulting increase in separation of weakly hydrophobic materials. The size dependence is analogous to case 1 and has already been discussed in detail there.

5. Discussion and Conclusions

The results demonstrate that, in theory, a multidimensional separation with respect to both surface properties and particle size is feasible with MSF. The presented theoretical study is in perfect agreement with experimental observations from previous studies [9,10]. Specifically, a separation based on dissimilarity in the separation criterion, represented by the charge-based separation of case 1, is strongly limited in material-specific selectivity. Naturally, all nonmagnetic target particles must differ in zeta potential for a selective separation, which directly results in an unavoidable affinity towards each other. This results in the formation and separation of multicomponent hetero-agglomerates, causing the separation of materials not intrinsically attracted to the magnetic component. Conversely, a separation based on similarity in the separation criterion, as exemplified by the hydrophobicity-based separation of case 2, does not encounter this problem, since hydrophilic particles are neither attracted to hydrophobic magnetic nor hydrophobic nonmagnetic particles and remain in suspension.

This work highlights the importance of considering the affinity of a material not only towards the magnetic particles, but also towards all other nonmagnetic materials. In both cases 1 and 2, nonmagnetic materials with no intrinsic attraction to the magnetic particles ($\alpha =0$) are investigated that would not be separated on their own. However, these materials are still separated to a large extent in case 1 due to an undesirable hetero-agglomeration with other nonmagnetic materials. Thus, the 2D grade efficiency function T depends not only on the separation criteria $\phi$ and d or general process parameters like pH and ionic strength, but also, and even more importantly, on the initial state of the suspension. Additional simulations illustrated this effect in Figure 7. Obviously, this effect is also present in other solid–liquid separation processes. During sedimentation, the settling velocity is often hindered by high solid concentrations [33,34,35], or during cake filtration, low solid concentrations result in slow cake formation and low initial performance [36]. However, this phenomenon is much more pronounced in MSF and presumably all agglomeration-based separation techniques. The initial state can essentially invert the selectivity of the process, posing a previously unrecognized challenge in separation technology. This raises critical questions about the applicability of the current definition of $T(\phi ,d)$ for agglomeration-based separation techniques, as $T(\phi ,d)$ shifts and changes depending on the initial state $n(\mathbf{x},0)$ it is applied to.

However, if a general 2D grade efficiency function $T(\phi ,d)$ cannot be defined, the following question arises: How can agglomeration-based processes be designed and controlled? This work presents MC simulations as an invaluable tool capable of predicting process outcomes based on static material functions like ${\mathit{\alpha }}_m$, process parameters, and—most importantly—initial conditions. Not only are MC simulations easy to set up, but they also outperform dPBE simulations in multicomponent systems. However, since dPBE models are significantly faster in lower dimensions (1D, 2D), a combined approach is proposed. Here, binary kernels ${\alpha }_{ij}$ are estimated from simple experiments in the one- and two-component system with the dPBE model, collected in ${\mathit{\alpha }}_m$, and subsequently transferred to the multicomponent system with the MC model. Future studies need to evaluate the predictive power and validity of the MC method with respect to experimental kernel data. Finally, it should be stressed that although discussed for MSF, all methods and results are applicable to any agglomeration-based process.