Real-Time Modeling of Volume and Form Dependent Nanoparticle Fractionation in Tubular Centrifuges

Abstract

A dynamic process model for the simulation of nanoparticle fractionation in tubular centrifuges is presented. Established state-of-the-art methods are further developed to incorporate multi-dimensional particle properties (traits). The separation outcome is quantified based on a discrete distribution of particle volume, elongation and flatness. The simulation algorithm solves a mass balance between interconnected compartments which represent the separation zone. Grade efficiencies are calculated by a short-cut model involving material functions and higher dimensional particle trait distributions. For the one dimensional classification of fumed silica nanoparticles, the numerical solution is validated experimentally. A creation and characterization of a virtual particle system provides an additional three dimensional input dataset. Following a three dimensional fractionation case study, the tubular centrifuge model underlines the fact that a precise fractionation according to particle form is extremely difficult. In light of this, the paper discusses particle elongation and flatness as impacting traits during fractionation in tubular centrifuges. Furthermore, communications on separation performance and outcome are possible and facilitated by the three dimensional visualization of grade efficiency data. Future research in nanoparticle characterization will further enhance the models use in real-time separation process simulation.

1. Introduction

Processing and purification of nanoscale products is becoming increasingly important in modern industry and science. One of the main reasons lies in the fact that nowadays, the connection between specific nanoparticle (NP) traits and general improvements of key product properties are well known [1,2,3,4,5,6,7,8]. Unit-operations regarding bottom-up synthesis [9,10,11] as well as top down formulation procedures [12,13,14] are therefore tailored towards high yields of narrowly distributed and beneficial particle sizes, shapes or material compositions. In many cases the use of additional separation steps to sort nanoparticles and remove impurities are mandatory [15,16]. In light of this, the demand for a multi-dimensional separation- and grade efficiency quantification in particle technology is increasing [17,18]. This means that several geometric and material parameters of a particle collective have to be evaluated simultaneously [19,20,21]. Following this goal, initial affords expand traditional schemes such as one dimensional (1D) particle size distributions (PSD) to e.g., two dimensional (2D) property distributions [22,23]. To clarify the notation used in this paper, the concept of particle trait distribution (PTD) is introduced. Traditional PSDs are identical to 1D-PTD because only one feature is highlighted. The consideration of multiple particle properties lead to higher dimensional PTDs, e.g., a three dimensional (3D) PTD.

In the field of solid–liquid separation, tubular bowl centrifuges are capable of processing nanosuspensions on a bench-scale with high throughput rates. Due to high centrifugal forces the apparatus can address challenging fractionation tasks of ultrafine particles and biological products [24,25,26,27,28,29]. The evaluation of the apparatus performance is based on the physical properties of a fine (centrate) and coarse fraction (sediment). Monitored process information usually include the solids mass fraction and a grade efficiency according to one major separation criterion (e.g., particle size). In an effort to include an additional criterion, Winkler et al. [30] studied the simultaneous fractionation according to particle size and solids density of mixed nanosuspensions experimentally. Moreover, the authors established a real-time monitoring methodology to evaluate the fine fraction composition [31]. The latter is necessary because tubular bowl centrifuges operate semi-continuously. The reason for this is that the accumulated sediment affects the separation result at longer process times. The experimental studies listed in this paragraph show that tubular bowl centrifuges can be successfully used for the fractionation of submicron- and nanosuspensions. However, they also underline that a large number of experiments, online monitoring and laboratory analyses are necessary to reliably characterize influencing variables.

To help understand the physical behavior within the centrifuge and gather vital process knowledge efficiently, mathematical modeling is important. Build upon the idea of flowsheet simulation [32,33], recent efforts established numerical models for time-dependent separation efficiency and sediment build-up computations in solid bowl centrifuges. These so called dynamic short-cut models describe complex physical conditions during the separation process on the basis of empirical equations and model assumptions. In contrast to a macroscopic process-chain consideration, these meso-scale models incorporate apparatus geometries and specific material functions to characterize the nanosuspension behavior. The well-documented methods of Gleiss et al. and Menesklou et al. [34,35,36,37,38,39] serve as the basis for a novel tubular centrifuge model (TCM) developed in this work. Regardless of the type and implementation of the process model, they share the same advantages. They can be used in combination with experimental data to optimize individual unit operations by identifying advantageous process parameters. Moreover, they reduce experimental efforts and assist in the design of experiments. Lastly, compared to simulations based on Computational Fluid Dynamics (CFD) and the Discrete Element Method (DEM), the dynamic process models demand a low investment in computational power. Consequently, this leads to a possible application in model predictive control strategies that demand real-time process simulation speeds [40].

One commonality of the highlighted studies is the exclusive focus on particle size as the major classification characteristic. In other words, solids were assumed to be perfect spheres, resulting in 1D grade efficiency curves. Revisiting the novel concept of multi-dimensional separation- and grade efficiency quantification stresses the point that existing short-cut models for tubular bowl centrifuges should be expanded to address e.g., form dependent fractionation. The presented study aims to merge the concept of form dependent sedimentation and multi-dimensional trait distributions with the short-cut modeling of NP fractionation in tubular bowl centrifuges. The adapted simulation tool is able to process both 1D and higher dimensional PTDs as a model input. Incorporating material functions and operating parameters, the spatial and time-dependent separation of NPs during centrifugation is computed.

This article is organized as follows: Section 2 outlines the theoretical background on form dependent NP sedimentation in tubular bowl centrifuges. Additionally, mathematical formulations of 3D joint distribution density functions and their application as PTDs are described. In Section 3, vital information on the TCM implementation, its working assumptions and parameters are presented. In a first simulation case, a 1D classification setup is used to validate the TCM functionality and output (Section 4.2) with experimental data. Here, commercially available fumed silicon dioxide (SiO${}_2$) NPs serve as a real-world particle system (RPS). The novel approach of 3D separation efficiency modeling is tested with the help of virtually generated SiO${}_2$ particles, further abbreviated as virtual particle system (VPS). Their creation and statistical evaluation according to particle form and volume is described in Section 3.4. The VPS and its 3D-PTD serves as new input for a second TCM simulation case. Material parameters and functions from the 1D case are reused. Based on the model output for 3D grade efficiencies, the new TCM approach is discussed thoroughly in Section 4.

2. Theory

2.1. Characterization of Particle Shape Based on Inertia Ellipsoids

The most frequently used isotropic quantity to describe a solid particle of arbitrary shape is the so called equivalent diameter,

$$d=\sqrt[3]{\frac{6}{\pi }{V}_{\mathrm{P}}},$$

(1)

of a sphere with identical volume $V_{\mathrm{P}}$. While it provides no morphological information, it can easily convey the particle mass if an affiliated solid density ${\rho }_{\mathrm{p}}$ is given. The particles morphology, however, can differ greatly in its degree of complexity. Both smooth, convex shapes such as spheres, plates, cuboids and ellipsoids as well as more complex aggregated structures and rough surfaces are possible. An extensive overview of different form parameters found in the literature is listed in [41]. Barrett et al. [42] mentions the three independent parameters form, roundness and surface texture one may consider for describing particle shape. Here, particle form is the easiest to express and measure since almost any given form parameter is estimated based on the longest $x_{\mathrm{L}}$, intermediate $x_{\mathrm{I}}$ and shortest $x_{\mathrm{S}}$ axis length. As shown in Zingg’s diagram [43], the two independent axis ratios of elongation,

$$e=x_{\mathrm{I}}/x_{\mathrm{L}}$$

(2)

and flatness

$$f=x_{\mathrm{S}}/x_{\mathrm{I}},$$

(3)

can express and visualize shape on a macroscopic level but exclude information on roundness and surface texture. Recent studies expanded this classification system and described a more detailed strategy to index non-spherical particles based on their elongation and flatness [44] differentiating between compact, bladed, elongated and flat objects.

In the scope of this work, the model strategy involves a simplification of any complex particle morphology based on a surrounding convex hull, hereby named Legendre${}^{\prime }$s ellipsoid. Rooted in the field of mechanics, it is also referred to as the inertia ellipsoid of a given compact body sharing the same moments of inertia and center of mass [45,46]. Figure 1 shows a 2D view of three complex 3D particle aggregates. Following the introduced classification criteria by Barrett et al. [42], their surface may be considered smooth, however, their overall structure is irregular and concave. The schematic also illustrates all associated inertia ellipsoids whose calculation is described in Section 3. From left to right, the first particle is elongated, the second represented by an oblate ellipsoid, and the one on the right is very compact. Throughout this paper, these scalene hulls and their axis lengths $x_{\mathrm{L}}$, $x_{\mathrm{I}}$ and $x_{\mathrm{S}}$ represent the global form of any given particle.

2.2. Form Dependent Particle Settling in a Centrifugal Field

Along the radial axis of rotation, a smooth particle in creeping motion within a liquid medium inside a centrifuge experiences mass and frictional forces. Assuming uncharged particles in an infinitely diluted suspension with no particle-particle interactions and negligible particle diffusion, the state of force equilibrium between drag, buoyancy, and centrifugal force yields

$$\begin{array}{c}\hfill \left(\frac{\mathrm{d}r}{\mathrm{d}t}\right)=u_0=\sqrt{\frac{V_{\mathrm{P}}r{\omega }^2\Delta \rho }{c_{\mathrm{D}}A\frac{1}{2}{\rho }_{\mathrm{l}}}}.\end{array}$$

(4)

Equation (4) provides an expression for the terminal settling velocity $u_0$ of a single particle with volume $V_{\mathrm{p}}$. Here, r defines the distance from the axis of rotation, $\omega$ is the constant angular velocity inside a stationary liquid pond and t stands for the elapsed time. Equation (4) describes the particle velocity as a function of the density difference between fluid and solid $\Delta \rho ={\rho }_{\mathrm{p}}-{\rho }_{\mathrm{l}}$ and the area A projected by the particles surface in the direction of radial flow. Lastly, the drag coefficient $c_{\mathrm{D}}$ is a dimensionless parameter which quantifies the solids hydrodynamic drag. Regarding the assumptions made to introduce Equation (4), the influencing factors of this parameter include the particle Reynolds Number,

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\mathrm{p}}=\frac{{\rho }_{\mathrm{l}}d{u}_0}{{\eta }_{\mathrm{l}}},\end{array}$$

(5)

geometric shape and orientation. Here, d is introduced as the equivalent diameter of a volume-equivalent sphere and ${\eta }_{\mathrm{l}}$ the dynamic viscosity of the fluid. Due to the small length-scales of colloids, a creeping motion with low Reynolds numbers $\left({\mathrm{Re}}_{\mathrm{p}}\ll 1\right)$ can be expected. In this regime, Stokes [47] found an analytical solution for spherical bodies with the definition $c_{\mathrm{D},\mathrm{S}}=24/{\mathrm{Re}}_{\mathrm{p}}$ for the drag coefficient, using Equation (1) and the projection area $A=\frac{\pi }{4}{d}^2$. Substituting these expressions in Equation (4) leads to the well-known expression of the Stokes settling velocity

$$\begin{array}{c}\hfill u_{0,\mathrm{S}}=\frac{d^2\Delta \rho r{\omega }^2}{18{\eta }_{\mathrm{l}}},\end{array}$$

(6)

for a spherical particle. In terms of particle shape, Equation (6) describes sedimentation in one dimension. The consideration of multiple geometric particle properties by drag correlations in several ranges of Reynolds numbers is extensively studied and discussed in the literature [41,48,49,50,51,52,53]. A recently published article by Trunk et al. [41] compiles multiple drag correlations for non-spherical particles and provides additional information on the drag correlation approximation for Reynolds numbers outside the Stokes regime. In the context of this work the mathematical formulation of the corrected drag coefficient,

$$\begin{array}{c}\hfill c_{\mathrm{D}}^{*}=f\left(R{e}_{\mathrm{p}},k_{\mathrm{S}}\right)=c_{\mathrm{D},\mathrm{S}}{k}_{\mathrm{S}},\end{array}$$

(7)

is defined as the product of $c_{\mathrm{D},\mathrm{S}}$ and the so called Stokes drag correction factor $k_{\mathrm{S}}$ in accordance with the definition used by Ganser et al. [49]. This simplification is only valid in the Stokes regime but applicable for this work, since the presented model is restricted to NP sedimentation. With the identical transformation strategy that was used for Equation (6), Equations (4) and (7) yield a approximation for the terminal settling velocity of non-spherical particles:

$$\begin{array}{c}\hfill u_0\left(k_{\mathrm{S}}\right)=\frac{d^2\Delta \rho r{\omega }^2}{k_{\mathrm{S}}18{\eta }_{\mathrm{l}}}.\end{array}$$

(8)

Besides Stokes’ settling theory, analytical solutions for $k_{\mathrm{S}}$ were obtained for spheroids [54] and other axisymmetric particle shapes [55,56]. The equations listed use synonyms of the form descriptors $x_{\mathrm{L}}$, $x_{\mathrm{I}}$ and $x_{\mathrm{S}}$ introduced in Section 2.1. Building upon this well documented knowledge, recent studies of Bagheri et al. [53] dealt with the prediction of the translational friction coefficient in gaseous and liquid media for a wide range of sub-critical Reynolds Numbers ($R{e}_{\mathrm{p}}<3\times {10}^5$). The authors claim, that shape descriptors defining $k_{\mathrm{S}}$ should be easy to measure and not dependent on the complex determination of e.g., sphericity in 3D scenarios. Consequently, they are defined by form- rather than roundness- or surface texture parameters which suits the purpose of this study. The presented definition of the drag correction factor

$$k_{\mathrm{S},\mathrm{B}}=\frac{1}{2}\left({\left(f{e}^{1.3}\right)}^{\alpha }+{\left(f{e}^{1.3}\right)}^{-\alpha }\right)$$

(9)

is merely a function of the particles elongation and flatness. The exponent $\alpha$ is set to $1/3$ following the assumption that particles tend to settle with neither a preferred orientation but rather with a random rotation in the Stokes regime [53]. Migration due to Brownian motion is common for NPs [57,58] which explains this assumption. Hence, Equation (9) is applicable when approximating the NP settling velocity in tubular bowl centrifuges. Following the model theory, the drag correction factor for strongly elongated and flat particles is high. Substituting $k_{\mathrm{S},\mathrm{B}}$ in Equation (8) with $\alpha =1/3$ yields

$$u_{0,\mathrm{B}}=\frac{d^2\Delta \rho }{6{\eta }_{\mathrm{l}}\left({\left(f{e}^{1.3}\right)}^{1/3}+{\left(f{e}^{1.3}\right)}^{-1/3}\right)}r{\omega }^2=s_{0,\mathrm{B}}r{\omega }^2,$$

(10)

for the sedimentation rate of arbitrary shaped, randomly oriented particles in creeping motion. Its diameter d and density define the mass of an equivalent spherical body while its shape induced deceleration is approximated by a scalene inertia ellipsoid with elongation e and flatness f. Additionally, Equation (10) introduces the sedimentation coefficient $s_{0,\mathrm{B}}$ as a particle trait dependent scalar with unit $s^{-1}$. Consequently, this parameter is dependent on the physical properties of both the solid and the surrounding fluid as well as the particle shape.

The presented expressions for the sedimentation coefficient are valid for highly diluted suspensions which have low technical relevance. Both density and viscosity of the suspension increase with increasing solids volume fraction. In addition, a solvent backflow is generated by the moving particles [59]. At high dilution, these hydrodynamic interactions induce a first order dependence of particle friction on the solids volume fraction $\varphi$ [60,61]. As a consequence an expression for the apparent sedimentation coefficient,

$$\begin{array}{c}\hfill s_{\mathrm{B}}=s_{0,\mathrm{B}}{\left(1+w\varphi \right)}^{-1},\end{array}$$

(11)

may be defined, where w is the hindered settling coefficient. For hard spheres, empirical observations are reported in the literature describing the influence of concentration on the settling behavior of mono- and polydisperse suspensions [62,63,64]. In the case of non-isometric particles shapes, only few experimental studies are reported in the literature. There are studies that have investigated the settling behavior of variously shaped macro-molecules or colloidal platelets. Again a linear relationship between the concentration and the normalized sedimentation velocity was found for low solid volume fractions [65]. Particularly noteworthy is the increase in particle resistance due to the decreasing sphericity and the resulting increased backflow of the solvent. Such results emphasize that the influence of concentration may not be neglected neither in experimental studies, nor in any model approach of colloidal particle sedimentation.

2.3. Joint Particle Trait Distributions

In particle technology, probability density distributions or PTD in general are used both for the characterization of a particulate product as well as for the quantification of a separation experiment outcome. In most cases, a PSD is derived either from the indirect analysis of a physical property such as the attenuation or scattering of light [66,67], by sieving, or by direct image processing [68,69]. Here, the equivalent particle diameter d is commonly used as a random variable $Z=\{d_1,d_2,...,d_M\}$, representing the arithmetic mean of each discrete particle size class in a predefined sample space. The distribution weighting is defined by the subscript $\mathrm{r}$. The probability P of a particle occurrence in interval $[d_1,d_2]$ is defined by:

$$\begin{array}{c}\hfill P\left(d_1\le Z\le d_2\right)=\underset{d_1}{\overset{d_2}{\int }}{q}_{\mathrm{r}}\left(d\right)\mathrm{d}d,\end{array}$$

(12)

for Z being a continuous random variable. In probability theory, $q_{\mathrm{r}}\left(d\right):\mathbb{R}\to \mathbb{R}$ is defined as a probability density function of Z which has to satisfy $q_{\mathrm{r}}\left(d\right)\ge 0$ and $\underset{-\infty }{\overset{\infty }{\int }}{q}_{\mathrm{r}}\left(d\right)\mathrm{d}d=1$. However, a dynamic process model for separation based on arbitrary shaped particles and their respective inertia ellipsoids calls for at least two random variables. Let Y be a second discrete random variable that takes on values $\{e_1,e_2,...,e_L\}$ and X be a third discrete random variable that takes on values $\{f_1,f_2,...,f_K\}$ representing the form parameters elongation and flatness. Then, the joint probability density function of the continuous random vector $(Z,Y,X)$ is $q_{\mathrm{r}}(f,e,d):{\mathbb{R}}^3\to \mathbb{R}$. The situation is analogous to the 1D case where the probability to meet a particle in a given interval,

$$\begin{array}{c}\hfill P\left(f_1\le X\le f_2,e_1\le Y\le e_2,d_1\le Z\le d_2\right)=\underset{f_1}{\overset{f_2}{\int }}\underset{e_1}{\overset{e_2}{\int }}\underset{d_1}{\overset{d_2}{\int }}{q}_{\mathrm{r}}(f,e,d)\mathrm{d}f\mathrm{d}e\mathrm{d}d,\end{array}$$

(13)

is equal to the volume of an infinitesimal cuboid of length $\mathrm{d}f$, width $\mathrm{d}e$ and height $\mathrm{d}d$ times the probability density at $(f,e,d)$. Again, $q_{\mathrm{r}}(f,e,d)\ge 0$ and

$$\begin{array}{c}\hfill \underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}{q}_{\mathrm{r}}(f,e,d)\mathrm{d}f\mathrm{d}e\mathrm{d}d=1\end{array}$$

(14)

must hold [70,71]. This means that for a particle trait distribution to be normalized neither the area nor the volume under a curve or plane must equal unity, but rather the sum of a four dimensional hyperspace. In Section 4, this is of interest for visualizing a particle trait distribution with three defining geometric variables.

Equations (12) and (13) can be formulated analogously as discrete distributions with sums of individual property classes, which are better suited for applications in separation technology. This is because numeric calculations involving PTDs are usually formulated in a discrete manner. In addition, experimentally determined distributions are given in discrete values, which result in an adequate histogram. In the scope of this work, volume weighted distributions are used for both particle trait distribution visualization and calculations in the presented dynamic process model for tubular bowl centrifuges.

One major difficulty associated with 3D distributions is the fact that their representations are often convoluted and lack clarity. To overcome these challenges in data visualization, marginalization is an adequate tool to reduce the dimensional complexity and potentially enhance the readability of higher dimensional data. In the following, the mathematical computation needed for the marginalization of 3D data to 2D and 1D datasets is shown. Referring to strategies used throughout this study we define:

$$\begin{array}{c}\hfill q_3(f,e)=\underset{d_{\min }}{\overset{d_{\max }}{\int }}{q}_3(f,e,d)\mathrm{d}d,\end{array}$$

(15)

as the marginalized volume weighted distribution density by integration over the equivalent particle diameter d. Analogously, performing a two step integration over elongation e and flatness f, a marginalized distribution density,

$$\begin{array}{c}\hfill q_3\left(d\right)=\underset{f_{\min }}{\overset{f_{\max }}{\int }}\underset{e_{\min }}{\overset{e_{\max }}{\int }}{q}_3(f,e,d)\mathrm{d}f\mathrm{d}e,\end{array}$$

(16)

in one dimension can be derived. The method of marginalization is used in Section 4 for discussion and data visualization purposes.

3. Materials and Methods

This section introduces the tubular bowl centrifuge used in both experimental and numerical investigations. Furthermore, the dynamic simulation tool is thoroughly described and all necessary balance equations and model assumptions regarding the time dependent NP separation are given. It builds upon the work of Gleiss [72] and published work of his research group at the IMVM in Karlsruhe. Lastly, physical information on two SiO${}_2$ NP suspensions are given. These particle systems are used to validate the 1D model output and help to understand the challenges in 3D and multi-dimensional separation setups.

3.1. Tubular Centrifuge Design Equations

This section deals with the mathematical description of the tubular centrifuge design equations used to model the separation process regarding single particles in motion. Figure 2 depicts a schematic cross section of a tubular centrifuge rotor. Typically, in most commercially available devices, the rotor is mounted in a suspended position and accelerated to a desired angular velocity $\omega$. The second important operating parameter is the volume flow that forces the suspension through the vessels bottom inlet. During operation, the centrifugal force forms a liquid pond, the height of which is set by the rotor weir at the outlet. Additionally, a gaseous core centered around the rotor axis is present. The exemplified trajectory in Figure 2 of particles with different sizes visualizes the apparatus general classification use. Coarser fractions accumulate at the rotor wall whereas smaller particles are transported beyond the overflow weir and enter the centrifuge downstream. Depending on the initial feed concentration, this build-up sediment can quickly reduce the available liquid cross section and diminish the separation efficiency over time.

Abstracted in two dimensions, the particles settling path in the fluid reservoir is set by its radial and axial velocity. The ratio of these two velocities is decisive for the separation efficiency approximation. By assuming plug flow inside any compartment with length $L_j$ in the liquid pond, the particles residence time,

$$\begin{array}{c}\hfill t_{\mathrm{v},j}=\frac{V}{\dot{V_{\mathrm{f}}}}=\frac{\pi \left(r_{\mathrm{b},j}^2-r_{\mathrm{w}}^2\right)L_j}{\dot{V_{\mathrm{f}}}},\end{array}$$

(17)

is altered by the liquid throughput ${\dot{V}}_{\mathrm{f}}$ and the surface area $A=\pi \left(r_{\mathrm{b},j}^2-r_{\mathrm{w}}^2\right)$ of the liquid layer. The weir $r_{\mathrm{w}}$ and boundary $r_{\mathrm{b},j}$ radius as well as the total length of the settling zone L are set by the accumulated sediment height and the rotor geometry.

Equating the particles residence and settling time by combining Equations (10), (11) and (17) yields an expression for a critical settling position,

$$\begin{array}{c}\hfill r_{\mathrm{c},j}=r_{\mathrm{b},j}·\left(exp\left(-s_{\mathrm{B}}·{\omega }^2·\frac{\pi \left(r_{\mathrm{b},j}^2-r_{\mathrm{w}}^2\right)L_j}{{\dot{V}}_{\mathrm{f}}}\right)\right),\end{array}$$

(18)

for each individual particle fraction. All material and fluid characteristics influencing the particles sedimentation are comprised in the introduced model for $s_{\mathrm{B}}$. This derivation is similar to that in publications of Gleiss et al. [35,36], who defined a grade efficiency:

$$\begin{array}{c}\hfill T_j\left(s_{\mathrm{B}}\right)=\frac{A_{\mathrm{sep}}}{A}=\frac{r_{\mathrm{b},j}^2-r_{\mathrm{c},j}^2}{r_{\mathrm{b},j}^2-r_{\mathrm{w}}^2},\end{array}$$

(19)

assuming a homogeneous distribution of each particle fraction along the radial axis. Here, $A_{\mathrm{sep}}$ is the area of the annular gap in which particles with an associated s-value will reach the rotor wall at the end of their residence time and are thus separated. As in Gleiss et al. [34] and Menesklou et al. [38], a compartment grade efficiency,

$$\begin{array}{c}\hfill T_j\left(s_{\mathrm{B}}\right)=T_j(d,e,f)=\frac{r_{\mathrm{b},j}^2}{r_{\mathrm{b},j}^2-r_{\mathrm{w}}^2}·\left[1-{\left(exp\left(-s_{\mathrm{B}}·{\omega }^2·\frac{\pi \left(r_{\mathrm{b},j}^2-r_{\mathrm{w}}^2\right)L_j}{{\dot{V}}_{\mathrm{f}}}\right)\right)}^2\right],\end{array}$$

(20)

is defined by inserting Equation (18) in (19). This expression for the grade efficiency gives an approximation for the separation probability as a function of the sedimentation coefficient, operating parameters and geometric boundaries. Additionally, this leads to a multivariate grade efficiency formulation centered at the three particle traits equivalent diameter d elongation e and flatness f that define $s_{\mathrm{B}}$. The consideration of two additional dimensions requires modifications to the discretization strategy used in the transient model and balance equations described in the next section.

3.2. Mathematical Model and Discretization

Building on the physical processes described in Section 3.1, Figure 3 highlights the spatial discretization of the separation zone in J compartments of equal volume. These sections are coupled by mass balances involving both suspended solids in the suspension zone (I) and deposited material in the sediment zone (II). Consequently, both particle sedimentation and sediment build-up are calculated locally with a desired spatial resolution. By using a finite amount of compartments, the temporal change in sediment height and particle residence time are approximated with low numerical effort [35].

One resulting advantage is that this strategy allows real-time modeling of the transient separation process. Referring to the illustration of compartment j in Figure 3, the mass balance equation of the suspension zone,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{m}_{\mathrm{I},j}}{\mathrm{d}t}={\dot{m}}_{\mathrm{I},j-1}-{\dot{m}}_{\mathrm{I},j}-{\dot{m}}_{j,\mathrm{sep}},\end{array}$$

(21)

describes the temporal mass change. Assuming a constant particle density, their mass can be substituted by the product of sedimentation zone volume $V_{\mathrm{I},j}$ and solids volume fraction ${\varphi }_{\mathrm{I},j}$ resulting in the expression:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left({\varphi }_{\mathrm{I},j}{V}_{\mathrm{I},j}\right)}{\mathrm{d}t}=\dot{V_{\mathrm{f}}}{\varphi }_{\mathrm{I},j-1}(1-E_{j,\mathrm{sep}}\left(t\right))-\dot{V_{\mathrm{f}}}{\varphi }_{\mathrm{I},j}.\end{array}$$

(22)

Here, the ratio of separated to incoming particle mass is introduced as the local separation efficiency $E_{j,\mathrm{sep}}\left(t\right)$. Analogously, the change in sediment volume $V_{\mathrm{II},j}$ over time,

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left({\varphi }_{\mathrm{II},j}{V}_{\mathrm{II},j}\right)}{\mathrm{d}t}=E_{j,\mathrm{sep}}\left(t\right)·\dot{V_{\mathrm{f}}}·{\varphi }_{\mathrm{I},j-1},\end{array}$$

(23)

is equal to the mass flow rate of the solids which do not migrate to the next compartment and are thus separated. For the calculation of $E_{j,\mathrm{sep}}\left(t\right)$ is calculated from $\dot{V_{\mathrm{f}}}$, $\omega$ and $s_{\mathrm{B}}$ as shown in Equation (20).

The general, discrete formulation of the separation efficiency,

$$\begin{array}{cc}\hfill E_{j,\mathrm{sep}}& \hfill ={\dot{m}}_{j,\mathrm{sep}}{\left({\dot{m}}_{\mathrm{I},j-1}\right)}^{-1}=\sum _{k=0}^K\sum _{l=0}^L\sum _{m=0}^M\left(T_j\left(f_k,e_l,d_m\right)q_{3,j}\left(f_k,e_l,d_m\right)\right),\hfill \end{array}$$

(24)

will now be further explained based on the introduced distribution theory in Section 2.3. Given any particle collective, a discrete, mass weighted formulation of Equation (12) yields:

$$m_j\left(d_m\right)=q_{3,j}\left(d_m\right)\Delta d\sum _{m=0}^M{m}_j\left(d_m\right),$$

(25)

for the mass fraction $m_j\left(d_m\right)$ of a particle with diameter $d_m$ in a class with uniform interval width $\Delta d$. M is the number of discrete diameter classes. Consequently, both distribution density $q_{3,j}\left(d_m\right)$ and compartment grade efficiency $T_j\left(d_m\right)$ (Equation (20)) are vectors with size $\left(1\times M\right)$. When incorporating two additional dimensions, the mass fraction,

$$m_j\left(f_k,e_l,d_m\right)=q_{3,j}\left(f_k,e_l,d_m\right)\Delta f\Delta e\Delta d\sum _{k=0}^K\sum _{l=0}^L\sum _{m=0}^M{m}_j\left(f_k,e_l,d_m\right),$$

(26)

includes $q_{3,j}\left(f_k,e_l,d_m\right)$ and $T_j\left(f_k,e_l,d_m\right)$ as 3D arrays of shape $\left(K\times L\times M\right)$. Both in the traditional 1D case and in the extended 3D consideration, $E_{j,\mathrm{sep}}$ is calculated from a numerical integration of the product of the separation efficiency and the distribution density as shown in Equation (24). The departing PTD is linked to the one in the prior compartment $\left(j-1\right)$ or in case of $j=1$ to the initialized feed distribution $q_{3,\mathrm{in}}$ as shown in Equation (27).

$$\begin{array}{c}\hfill q_{3,j}\left(f_k,e_l,d_m\right)=q_{3,j-1}\left(f_k,e_l,d_m\right)\frac{\left(1-T_j\left(f_k,e_l,d_m\right)\right)}{\left(1-E_{j,\mathrm{sep}}\right)}\end{array}$$

(27)

Consequently, Equations (22) and (23) form a linear system of equations solved by the algorithm depicted in Figure 3b. For its implementation, the python3 programming language is used. The list of mandatory repositories include the numpy [73] package and function fsolve [74,75] from the scipy package. At the start of a simulation run, the compartment count J, the initial PDT, solid and liquid parameters as well as the rotor geometry is set. The introduced system of balance equations is solved implicitly for a user defined time step size $\Delta t$ and total simulation time $t_N$. A solution is found for each $t_n$ if the relative residual error between two consecutive iterates is at most ${10}^{-6}$. Due to the progressive reduction of the volume $V_{\mathrm{I},j}$, the boundary radii $r_{\mathrm{b},j}$ increase and the local grade efficiencies must be recalculated for each time step. The first salient model output is the centrates solids volume fraction expressed as ${\varphi }_{\mathrm{out}}\left(t_n\right)$. With the known feed concentration ${\varphi }_{\mathrm{in}}$ an expression for the global separation efficiency,

$$E_g\left(t_n\right)=1-\frac{{\varphi }_{\mathrm{out}}\left(t_n\right)}{{\varphi }_{\mathrm{in}}},$$

(28)

is set. The second model output calculation centers around the global grade efficiency,

$$\begin{array}{ccc}\hfill T_g\left(d_m,t_n\right)& \hfill =1-\frac{E_g\left(t_n\right)q_{3,\mathrm{out}}\left(d_m,t_n\right)}{q_{3,\mathrm{in}}\left(d_m\right)}\mathrm{and}\hfill & \hfill \left(1\mathrm{D}\right)\end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\tilde{T}}_g\left(f_k,e_l,d_m,t_n\right)& \hfill =1-\frac{{\tilde{E}}_g\left(t_n\right){\tilde{q}}_{3,\mathrm{out}}\left(f_k,e_l,d_m,t_n\right)}{{\tilde{q}}_{3,\mathrm{in}}\left(f_k,e_l,d_m\right)}\hfill & \hfill \left(3\mathrm{D}\right),\end{array}$$

(30)

for one up to three dimensions following the fundamental grade efficiency descriptions found in the mechanical process engineering literature [76]. With Equations (28)–(30) each time step provides information on the actual separated mass per particle trait class based on the outgoing, mass weighted density distribution $q_{3,\mathrm{out}}$.

Regarding transparency and the discussion of model caveats in later sections, the following list of assumptions and model simplifications is noted:

(1)

The flow pattern inside the centrifuge is plug flow. solid–liquid or solid-solid interactions are modeled solely in one direction by the hindered settling coefficient (Equation (11));

(2)

Particles are small compared to the rotor length and weir perimeter. Fluid transport is therefore not influenced by the dispersed NPs;

(3)

A pronounced radial dispersion within the liquid layer is assumed. In each compartment the incoming particle fractions are back-mixed ideally;

(4)

Any formed sediment is incompressible and maintains a constant porosity

$\left(\frac{\mathrm{d}{\varphi }_{\mathrm{II},j}}{\mathrm{d}t}=0\right)$ simplifying the implicit solution of Equations (22) and (23) and reducing the computational effort;

(5)

The gas-liquid interface has no influence on the separation process;

(6)

The suspension is stabilized and particles retain their geometric shape and material density. Hence, no breakage or agglomeration takes place resulting in an immutable PTD both in the 1D and 3D case.

One note concerning the code structure and its execution is that both classification and form dependent fractionation can be addressed. Inside the computational domain, multiple TCM cases with different initial boundary conditions, inputs and discretization parameters can be set up. If the input PTD is a 1D-array holding discretization variable d, a classification is calculated and the assumption $e=f=1$ for each particle class holds. The second considered setup is that the input PTD represents a tensor. Here, d, e and f define the volume and form properties of the respective particle classes. All relevant cases presented for this paper and their parameters, inputs and outputs are described in Section 4.

3.3. Experimental Setup and Model Particle System Characterization

The experimental setup includes a Z11-type tubular centrifuge designed by Carl Padberg Zentrifugenbau GmbH (CEPA, Lahr, Germany). All geometric boundaries used in the simulation tool are listed in

Table 1. Overview of the centrifuge geometry, fluid and particle system properties.

Parameter	Symbol	Value	SI unit
weir radius	$r_{\mathrm{w}}$	0.0073	$\mathrm{m}$
wall radius	$r_{\mathrm{b}}\left(t_{n=0}\right)$	0.0215	$\mathrm{m}$
rotor length	L	0.175	$\mathrm{m}$
angular velocity	$\omega$	40,000	${\mathrm{s}}^{-1}$
flow rate	${\dot{V}}_{\mathrm{f}}$	$5\times {10}^{-6}$	${\mathrm{m}}^3{\mathrm{s}}^{-1}$
solids density	${\rho }_{\mathrm{p}}$	2200.0	$\mathrm{kg}{\mathrm{m}}^3$
mean particle size	$d_{50}$	64.39	$\mathrm{nm}$
solids volume fraction	${\varphi }_{\mathrm{in}}$	$3.0\times {10}^{-4}$	−
liquid density	${\rho }_{\mathrm{l}}$	998.207	$\mathrm{kg}{\mathrm{m}}^3$
liquid viscosity	${\eta }_{\mathrm{l}}$	$1.0027\times {10}^{-3}$	$\mathrm{kg}{\mathrm{m}}^{-1}{\mathrm{s}}^{-1}$

and predefined by the apparatus rotor itself. An eccentric screw conveyor (NETZSCH Pumpen und Systeme GmbH, Waldkraiburg, Germany) is able to transport a suspension with low pulsation and constant flow rate to the rotor inlet. For further information regarding the process layout for NP separation, the reader is referred to a preliminary study [31] which incorporated a similar setup.

To test and validate the adapted process model for 1D-PTDs, a non-spherical particle system is used in this study. The trade name Aerodisp^® W7512 S specifies a hydrophilic, fumed silicon dioxide (SiO${}_2$) suspension manufactured by Evonik Industries AG (Hanau, Germany). The product is synthesized in a flame hydrolysis process forming primary particles which then collide and partially attach to each other. This aggregation step is followed by an agglomeration of the non-spherical aggregates due to weak inter-particular interactions. Some examples showing the branched structure of the individual aggregates can be found here [29,77,78]. The SiO${}_2$ particles are well dispersed in demineralized water resulting in $30\mathrm{nm}\le d\le 130\mathrm{nm}$ for the expected equivalent aggregate diameter. To ensure this dispersion state, the suspension is sonified (Ultrasound processor UP400St, Hielscher Ultrasonics, Teltow, Germany) prior to any classification experiment. The volume weighted, 1D-PTD was measured in a CPS 24,000 disk centrifuge (CPS Instruments Inc., Prairieville, LA, USA). The measuring principle of this analytical disk centrifuge (ADC) is based on an extinction screening at the outer edge of a rotating disk. A laser with wavelength 470 nm traverses through the transparent disk and sample at a small area before hitting the detector. Large particles reach this position earlier, successively followed by smaller particles. The advantage of the ADC measurement is that a density gradient in the disc effectively separates multiple size fractions, thereby reducing inter-particular forces by dilution and stabilizing the sedimentation during analysis. Using predefined procedures for each material, the raw extinction signal is translated into a 1D-PTD by combining Equation (6) and Mie scattering theory. Hence, only the equivalent diameter d and its distribution $q_{3,\mathrm{in}}\left(d\right)$ is measurable. This inevitably reduces the aggregates characterization to 1D. Nevertheless, referring to Figure 3 and Equation (25), the measured 1D-PTD of the RPS serves as an input for the 1D-TCM.

In order to assess the influence of particle concentration on the sedimentation behavior, the multisample analytical centrifuge (AC) LUMiSizer (LUM GmbH, Berlin, Germany) was utilized. The SiO${}_2$ suspension with a defined solids volume fraction was filled into a small cuvette and inserted onto a rotating platform. Their filling volume was 0.4 ml, spinning with a rotational speed of 4000 min${}^{-1}$ at a constant temperature of $\delta =$ 10 °C. Lowering $\delta$ in the measurement chamber helps to reduce particle diffusion. During sedimentation analysis a light source emits parallel rays of light [79] with the same wavelength as the ADC ($\lambda =470$ nm). While the rotation forms a meniscus in the perimeter of the cuvette, the AC measures an extinction profile along its entire height. Logging a new curve every 30 s provides a progression of sedimentation over time. The raw AC data can be used to calculate an apparent sedimentation coefficient distribution $g^{*}\left(s_{\mathrm{ext}}\right)$ either by a spatial or temporal differentiation of these extinction profiles [80,81]. In this study, the method published by Schuck et al. was used to derive $g^{*}\left(s_{\mathrm{ext}}\right)$ by means of least-squares boundary modeling (LSBM) [82]. The described procedure was first translated to a python code which automated extinction data input, pre-processing and $g^{*}\left(s_{\mathrm{ext}}\right)$ calculation. Note that the procedure neglects diffusion effects and involves an extinction weighting of the sedimentation coefficient. Unlike $s_{\mathrm{B}}$ introduced in Section 2.2, the extinction weighted distribution $g^{*}\left(s_{\mathrm{ext}}\right)$ can not be used directly for the TCM model but is sufficient for an empirical approximation of the sedimentation hindrance coefficient w in Equation (11).

3.4. Virtual Particle System Creation and Characterization

This section describes the workflow used to create a virtual SiO${}_2$ particle system (VPS) and compute an affiliated 3D-PTD. The entire method is structured in three sub-modules depicted in Figure 4.

The first module (M1) tackled the computer simulation of particle aggregation. From a general point of view, the formation of SiO${}_2$ aggregates during flame synthesis can be described by different kinetic growth models [83,84]. For the presented study, a custom diffusion-limited-aggregation (DLA) algorithm was used.

The computation starts with a seed particle in the center of a cubic grid. A second particle is released in the domain and moves step-wise to one of 26 possible lattice coordinates. If the randomized path intersects the initial seed perimeter, a collision is recorded. For a randomized total primary particle count, this process is repeated until a final 3D structure is build. Particles that wander outside the boundary are deleted and respawn [85,86]. A unique feature of the presented algorithm is that three individual spawn locations are defined. Usually, a particle trajectory starts at the surface of a sphere that surrounds the seed particle (option 1). Two additional options restrict the particles entry positions to the spheres equatorial circumference (option 2) and both polar regions (option 3). In each new run, the algorithm randomly chooses one option. This additional degree of randomization ensures that each particle has a unique structure with a lesser occurrence of compact aggregates. A more detailed and visual description of the custom DLA algorithm is provided as Supplementary Material of this manuscript. The setup was used to construct a collection of $P=3771$ aggregates for the VPS. Visual examples for some of the generated particles in the collection are shown in Figure 5.

The second module (M2) conducts a volumetric analysis of multiple 3D objects at once. Its functionality is inspired by quantitative computer tomography (CT) measurements of particles on the micro- and nanoscale [87]. In particle technology, nano- and micro-CT measurements are reconstructed into 3D images of a representative sample in a defined environment. Subsequently, particles in 2D image slices are segmented and counted. Each observed particle can be labeled with, for example, a vector of geometric properties. Detailed information on the state-of-the-art execution of CT-measurements and image processing workflows can be found here [88].

To recreate this methodology for the VPS, the open source content creation suite Blender [89] and its build-in Python application programming interface (API) was used to place roughly 100 particles out of the whole collection in a rectangular volume with a base area of 3.3 µm${}^3$. Thereby, care was taken that the individual particles did not overlap. At the end of this preparation procedure, 37 samples were exported in the Standard Triangle Language (STL) format and imported separately into the image processing software Fiji [90]. Afterwards, a plugin called SciView [91] created image stacks of each individual sample. An image resolution of 600×600 pixels results in a voxel size of 5.5 nm. At this point it should be mentioned that the technically feasible spatial resolutions in nano-CT are generally lower than the one applied in this module. To give an example, Ditscherlein et al. [88] reported applicable voxel sizes of 16 nm and 64 nm depending on the nano-CT measurement type. The applied spatial resolution in this study ensures an accurate recording and quantitative analysis of the nanocluster form. Processing the resulting binary image stacks in the 3D Objects Counter-plugin [92] created an indexed voxel representation of every virtual sample. The obtained object maps were then analyzed by the well-known mathematical morphology analysis plugin MorphoLibJ [93]. Of particular interest is the automated determination of the 3D inertia ellipsoid and voxel volume. This volume in combination with a predefined solids density provides an approximation of the mass of each individual particle in the collection. As stated in Section 2, both particle volume and the form of a convex hull defined by the inertia ellipsoid are needed for sedimentation coefficient calculations. Therefore, the transverse and polar radii were used to calculate the particles elongation and flatness. Finally, the discrete particle information for $V_{\mathrm{P}},f,e$ and d were compiled in a list with $P=3771$ entries.

The last module (M3) performs a statistical particle trait analysis based on a given input list of discrete information and outputs a user defined PTD. In accordance to Section 2.3 a 3D-PTD in discrete notation is a tensor of size ($K\times L\times M$). Each of the three variables e, f and d have user defined limits and discretize one dimension. During computation, an algorithm sorts each particle n with respect to their trait combination ${\{f,e,d\}}_n$ into a matching bin. Afterwards, the relative particle mass $m\left(f_k,e_l,d_m\right)$ in each class is calculated by adding up all particle volumes and multiply the sum class-wise with a predefined solids density ${\rho }_{\mathrm{p}}$. In the final step, Equation (26) is used to compute the volume weighted 3D-PTD. This routine was applied to the VPS trait database output of M2 resulting in a multi-dimensional density distribution ${\tilde{q}}_{3,\mathrm{in}}\left(f_k,e_l,d_m\right)$ which served as one input for the 3D-TCM model.

The assignment of discrete property vectors based on high fidelity image analysis is a common practice in material science. For example, methods such as transmission (TEM) [94,95] or scanning electron microscopy (SEM) [96,97] are able to statistically characterize complex particle systems on a broad size range. In case of 3D image processing, Ditscherlein et al. [98] introduced an open access archive for particle-discrete CT. One major advantage of module M3 in the presented workflow is that it can process such database entries independently of the measurement methodology used. This means that, in future research, real-world particle data can be easily translated into a 3D-PTD and implemented into a 3D-TCM.

Modeling the fractionation of NPs is not the only application for particles emerging from the presented workflow. As an additional advantage, the universal data structures of the discrete particles leads to possible applications in other research areas. One example is the permeability investigation of porous media, which is largely determined by the spatial structure of the capillaries [99]. Similarly, numerical simulations on the performance of lithium ion batteries by Chauhan et al. can be highlighted in this regard. In their work, the authors generate idealized microstructures which mimic the spatial distribution of key materials such as carbon black and active materials in the cathode [100]. Here, the presented workflow of this study could offer an alternative to generate discrete 3D particles with variable microstructures. Consequently, the volumetric and statistical analysis (M2, M3) would add an additional factor to the performance simulation of lithium ion batteries: distributed particle features in more than one dimension.

4. Results

The presented results are divided into two parts. First, the experimental determination of the hindered settling coefficient is presented. Afterwards, a 1D- and 3D-TCM case output is shown. In the 1D setup, a 1D-PTD with 500 discrete bins for d serves as a density distribution input measured by an ADC. In the 3D case, a PTD with a total of 27,000 bins is utilized. For its construction, the VPS trait database provides all required information. Note that for reference and clarity all model outputs of the 3D case are marked with an accented tilde (e.g., $\tilde{T}$).

4.1. Hindered Settling Evaluation of Aerodisp^® W7512 S

Following the steps described in Section 3.3, the RPS suspension was diluted with demineralized water and analyzed in the AC. Raw extinction data import and processing via LSBM leads to an apparent, extinction weighted sedimentation coefficient distribution. The advantage of this method is that the sedimentation does not need to be significantly advanced in order to derive a solution for $g^{*}\left(s_{\mathrm{ext}}\right)$, resulting in less time required per measurement. The evaluation set-point is chosen so that all particles are in motion, jet no complete separation of larger fractions takes place. By numerical integration of $g^{*}\left(s_{\mathrm{ext}}\right)$, the cumulative sum $G^{*}\left(s_{\mathrm{ext}}\right)$ was determined and plotted with its standard deviation in Figure 6a. Each curve is based on four analyzed samples and the mean value is plotted.

As described in Section 2.2, $s_{\mathrm{B}}$ is a function of the sedimentation coefficient $s_{0,\mathrm{B}}$ in infinite dilution, the solids volume fraction and the hindrance coefficient. In a subsequent step, the characteristic of the abscissa at $G^{*}\left(s_{25,\mathrm{ext}}\right)=0.25$, $G^{*}\left(s_{50,\mathrm{ext}}\right)=0.5$ and $G^{*}\left(s_{75,\mathrm{ext}}\right)=0.75$ were determined and plotted reciprocally against the solids volume fraction in Figure 6b. In other words, the intersections between each of the horizontal lines and the $G^{*}\left(s_{\mathrm{ext}}\right)$ cumulative sum distribution in Figure 6a lead to the 18 plotted datapoints shown in Figure 6b. Assuming the linear dependence,

$$\frac{1}{s_{\mathrm{ext}}}=\frac{w}{s_{0,\mathrm{ext}}}\varphi +\frac{1}{s_{0,\mathrm{ext}}},$$

(31)

after rearranging Equation (11), a least squares fit yields w and $s_{0,\mathrm{ext}}$ as functional parameters [101,102]. A parallel evaluation for all three distribution reference points leads to a mean hindered settling coefficient of $w=32.44\pm 9.40$ for the SiO${}_2$ suspension. The uncertainties in w can be explained by error propagation. The averaging of four individual measurements from the LSBM leads to noticeable slope deviations marked as gray areas beneath the linear fit curves in Figure 6b. Nevertheless, the AC analysis provides a correction for the sedimentation conditions based on inter-particle interactions, which are enhanced with an increased solids volume fraction. Therefore it is used in both the 1D- and 3D-TCM for sedimentation hindrance modeling.

4.2. TCM Results for 1D Classification

This part illustrates a 1D classification experiment of the RPS used to validate the 1D-TCM. Utilizing the experimental setup and parameters highlighted in Section 3.3 and Table 1, the analyzed 1D-PTDs of both feed suspension and centrate are shown in Figure 7a. Additionally, ADC measurements yield an experimental value for the global separation efficiency $E_{\mathrm{g}}$ used to calculate a 1D graph of $T_{\mathrm{g}}$ (Equations (28)–(30)). In general, experimentally determined values are shown as points whereas the 1D-TCM simulation output is marked as solid lines. The separation experiment was performed twice and both ADC and gravimetric analysis were performed in triplicate to derive the mean and standard deviation for ${\varphi }_{\mathrm{in}}$, ${\varphi }_{\mathrm{out}}$ and $T_{\mathrm{g}}$. Results for both $q_{3,\mathrm{out}}$ and grade efficiency are in good agreement with comparable classification scenarios found in the literature [27].

Continuing with the validation procedure, identical boundary conditions and the acquired 1D-PTD input $q_{3,\mathrm{in}}\left(d\right)$ supply a working 1D-TCM case as proposed schematically in Figure 3. The total simulation time $t_N$ is limited to 360 s to match the suspension sampling time in the experimental study. Following a preliminary evaluation procedure described in Appendix A.1, the spatial discretization is set to $J=20$ compartments and the simulation timestep is set to $\Delta t=0.5s$.

Figure 7b highlights the history of $E_{\mathrm{g}}$ (solid line) as well as the experimental reference acquired by differential weighing (DW) and ADC measurements. The initial model output indicates a separation efficiency of one since the rotor is pre-filled with demineralized water. After 30 seconds, the suspension and sediment zone gradually populate with solid material due to particle entrance and deposition. After roughly one minute, the process runs practically stationary with almost no change in $E_{\mathrm{g}}$. This is explained by the low feed concentration of the SiO${}_2$ suspension resulting in no significant sediment buildup and effect on sedimentation behavior. Likewise, the influence of sedimentation hindrance is barely noticeable since measured uncertainties in w result in a 0.48% fluctuation of the simulated mean cut size $d_{T,50}$. This conclusion surfaced after a sensitivity analysis in which w was adjusted multiple times within the scope of its standard deviation. This is in accordance to experimental results conducted by Spelter et al. [29] describing a negligible effect of sedimentation hindrance for ${\varphi }_{\mathrm{in}}\le 0.5v\%$ for the identical particle system. In summary, numerical simulation results and experimental data are in a very good agreement. Most discernible deviations can be explained by a very low concentration of NPs at the beginning ($d\ge$ 150 nm) and at the end ($d\le$ 35 nm) of a 1D-PTD analysis with the ADC. Overall, these findings underline the validity of the current results in case of the 1D-TCM model. Chances made to the initial dynamic process model introduced by Gleiss [72] did not affect the congruence of experiment and simulation.

The presented TCM model version is tailored to flexible, multi-dimensional PTD inputs and short simulation times. This means that if long-term experiments are to be simulated, model accuracy could be worse. The importance of sediment build-up consideration has already been proven experimentally by several authors [29,103]. Therefore, Gleiss et al. [37] and Menesklou et al. [39] implemented an adequate method building upon empirical measurements for the consolidation mechanism in decanter centrifuges. For now, the need for long-term model accuracy is secondary to the successful implementation of 3D-PTD into the TCM model. Future research, however, needs to incorporate polydisperse, multi-component sediment build up and consolidation with adequate material functions. Nonetheless, in the presented work, a high model accuracy was achieved for the numerical simulation of 1D classification in tubular centrifuges.

4.3. Visualization of 3D-PTDs

The pipeline output presented in Section 3.4 is best described as a 3D array that holds the PTD information. Unlike 1D-PTDs usually found in particle technology, their graphical visualization is more challenging. Figure 8a shows the computed PTD of the VPS. Here, discretization settings with $K=L=M=30$ nets 27,000 individual trait classes. Data presentation resembles a cloud where each individual bin is color in accordance to a certain probability density ${\tilde{q}}_3(f,e,d)$. High values (yellow color) indicate that particles in this trait class are more likely to be found during a quantitative sample analysis of the whole collection. Empty bins are rendered as transparent volumes. The color-map on the left indicates the volume weighed density distribution. An advantage of this visualization method is a fast identification of elongated, flat or bladed NP populations. In the given example, both compact particles with e.g., $f=e\ge 0.9$ and elongated aggregates with $e\le 0.3$ and $f\ge 0.8$ are present in the sample and can be detected with the naked eye. A severe disadvantage poses the noticeable cluttering of density distribution data. Unless the 3D plot is viewed from multiple angles, parts of the distribution are obscured. At the same time, the interior data are not visible at all. As described in Section 2.3, one possible solution is distribution marginalization. As an example, Figure 8c shows the 2D-PTD obtained by Equation (15). In this graphical representation of the VPS, one can identify the presence of both elongated and more compact aggregates in the range of $30nm\le d\le 120nm$. Hence, in combination with the other two marginalization options ${\tilde{q}}_{3,\mathrm{in}}(f,d)$ and ${\tilde{q}}_{3,\mathrm{in}}(e,d)$, the discussion of discrete, higher dimension probability data is simplified.

Another point to address is the distribution resolution or in parameter expression, the total discrete class count. One negative aspect of a high PTD resolution is linked to an elevated computational effort in the 3D-TCM (see Figure A2). However, a low resolution with a total class size below 100 inevitably leads to higher discretization errors, since detailed particle trait information is lost. A preliminary study justifying the discretization settings used for the 3D-TCM case is given in Appendix A.2.

4.4. TCM Results for 3D Fractionation

This section covers the 3D-TCM case output for the VPS. The simulation domain is identical to the presented 1D process simulations. This means that all boundary settings mentioned in Section 4.2 stay fixed. The simulation time is set to 360 s, the total compartment count is $J=20$ and the time step size equals 0.5 s. Additionally, the empirical hindered settling coefficient (see Section 4.1) is reapplied. The 3D-TCM output results are outlined in two figures: The presentation in Figure 9 is identical to the one in Figure 8, since both ${\tilde{T}}_{\mathrm{g}}\left(f,e,d,t=360\right)$ and ${\tilde{q}}_{3,\mathrm{out}}\left(f,e,d\right)$ share the same data structure and array shape. The 3D grade efficiency illustration, acquired by Equation (29), implies the same negative aspects regarding readability on 2D prints. Nevertheless, it can be seen that particle classes with larger values for d tend to be separated more effectively in the centrifuge. This is indicated by a more saturated coloration of the individual classes. Since the equivalent diameter reflects the particles volume (Equation (1)) which dominantly impacts particle sedimentation rates according to Equation (4), this observation is expected. The same trend can be noticed in Figure 9b. Here, the grade efficiency,

$$\begin{array}{cc}\hfill {\tilde{T}}_g\left(e_l,d_m,360s\right)& \hfill =1-\frac{{\tilde{E}}_g\left(360s\right){\tilde{q}}_{3,\mathrm{out}}\left(e_l,d_m,360s\right)}{{\tilde{q}}_{3,\mathrm{in}}\left(e_l,d_m\right)},\hfill \end{array}$$

(32)

is calculated from the marginalized PTD of the feed and centrate suspension, isolating the traits elongation and particle volume. A marginalization to 2D draws attention to the fact that SiO${}_2$ grade efficiencies are slightly larger for more elongated particles ($e\to 0$). For a guide to the eye, classes that hold a value of $0.475\le {\tilde{T}}_g\left(e_l,d_m,360s\right)\le 0.525$ are marked with a single black dot.

Revisiting the presented form analysis (Section 3.4) of all 3D aggregates which populate the VPS can explain these observations. Figure 10 shows $k_{\mathrm{S},\mathrm{B}}$-values in dependence of elongation and flatness calculated with Equation (9). Additionally, for each individual particle, $P=3771$ black dots are plotted. Their position is set according to their e and f values which were recorded by the volumetric analysis of M2. Looking at the spatial distribution of all dots, only a small minority of the particles generated by the custom DLA algorithm can be associated with $k_{\mathrm{S},\mathrm{B}}$-values greater than 1.2. Consequently, the sedimentation coefficient (Equation (10)) for a majority of the particles is similar to that of a sphere with equal volume. As a result, the overall impact of particle form on the model output is almost negligible. Despite this fact, it should be positively emphasized that the separation of several different form fractions can be identified with this method. To give an example, one could imagine a particle collective of small, compact and larger rod-like particles. Due to their differences in elongation, their fractions are easily identified in a 3D-PTD. Referring to the presented output of any 3D-TCM model, their separation efficiency can be calculated and visualized intuitively with the presented methodology. In summary, it can be said that the integration of Equation (10) into the TCM mass balance calculation is functional. Both volume and particle form, represented by a surrounding inertia ellipsoid with elongation e and flatness f, is considered in the calculation of the sedimentation coefficient. The simulation-tool is able to visualize the outcome of a tubular centrifuge separation experiment in more than just one dimension even for small deviations in particle form.

Closing the gap between the 3D-TCM and a traditional separation process outcome, the marginalized 1D-PTD and the associated grade efficiency curve is shown in Figure 11. A bar-plot is used to view the 30 classes which discretize the equivalent diameter. Blue bars and their height reference the inlet PTD whereas red bars show the product of ${\tilde{E}}_{\mathrm{g}}$ and the marginalized outlet PTD. One important question to address is whether the global grade efficiency in the 3D case is comparable to the preceding RPS classification experiment. This should be the case since both the RPS and VPS share the same solid density and a similar equivalent diameter range. To compare both conducted TCM cases, the 1D global grade efficiencies $T_{\mathrm{g}}$ and ${\tilde{T}}_{\mathrm{g}}$ are added as line plots to Figure 11. As expected, the data are congruent with a slight deviation for higher particle volumes. This outcome can be explained by the presence of elongated particles with $d\ge 80nm$ in the collection. For a visual reference, attention is drawn to the upper left quadrant of ${\tilde{q}}_{3,\mathrm{in}}\left(d\right)$ in Figure 8b and again to Figure 10. According to the fundamentals of form dependent NP settling, these discrete classes have trait combinations which lead a lower value for $s_{\mathrm{B}}$.

Hence, during iterative short-cut computation of $T_j\left(s_{\mathrm{B}}\right)$, these mass fractions are separated less effectively. In the 1D setup, however, discrete sedimentation rates are higher due to the strict limitation of $k_{\mathrm{S},\mathrm{B}}=1$. As a result, the mean cut size is shifted slightly towards smaller particles. Since the VPS is considered a unique particle system, there are no comparable experimental or model results present in the literature. Hence, only the plausibility of the results are validated in this study. A detailed validation of 3D grade efficiencies will demand precise and adapted methodologies regarding the multi-dimensional analysis of NP properties used in future research. In light of this, the presented work establishes an adaptable routine for the creation and detailed trait analysis of NP collectives. Consequently, the proposed workflow and working adjustments made to the TCM facilitate these future efforts to conduct an experimental validation in multi-dimensional fractionation.

The results imply that the novel simulation routine is capable of addressing the NP form fractionation in tubular bowl centrifuges at a fixed set of boundary condition. However, the added benefit of the adapted multi-dimensional TCM does not end here. Model parameters such as fluid viscosity, angular velocity and volumetric flow rate can be varied with little effort resulting in an updated separation result. A completely automated adjustment of these boundary conditions is also conceivable and can be implemented by only a few changes in the program code. This is an important prerequisite for the efficient performance of parameter studies that can quantify the influence of individual process parameters. Since the novel model also calculates and displays 3D separation efficiency data (see Figure 9), optimization with respect to one or more particle features is possible. Lastly, the feed material can be easily modified, either based on real-world or computer-generated particles derived from simulated (see Section 3.4 and Figure 4) or real-world discrete particle databases [98].

5. Conclusions

Traditionally, the separation efficiency in tubular bowl centrifuges is evaluated as a function of particle size, represented by the radius of a volume equivalent sphere. However, current research draws attention to advantageous product properties linked to, for example, the form distribution of a particle collective. Therefore, this work develops a simulation tool for fractionation modeling of multi-dimensional NP systems in tubular bowl centrifuges. To facilitate the numerical calculation of both one- and multi-dimensional grade efficiencies, an existing model by Gleiss [72] is modified. In its core, the dynamic model solves mass balance equations for a finite number of compartments that spatially discretize the sedimentation zone and define the solids residence time.

Investigations first addressed the 1D classification of a real-world SiO${}_2$ nanosuspension. In a direct comparison, numerical simulations predicted the separation outcome with very high accuracy. In both the 1D-TCM and reference measurements for the achievable grade efficiency only the equivalent diameter was considered a particle trait.

For a second simulation case, 3771 virtual SiO${}_2$ aggregates were created by a custom DLA algorithm and characterized individually by direct image analysis. Their form and volume is discretized by three variables: elongation, flatness and equivalent diameter. By means of an automated statistical trait analysis, a volume weighted PTD of the particle system is calculated which served as the new TCM model input in a 3D fractionation setup. Here, form-dependent settling was implemented by empirical correlations found in the literature. Numerical simulations confirm that particle form has a much smaller influence on their sedimentation rate compared to particle volume. The presented discussion implies that the separation result per particle trait class is only influenced by a more pronounced elongation or flatness $\left(e<<1,f<<1\right)$. With regard to the precise adjustment of the product properties through, for example, modifications to the form distribution, a major challenge in multi-dimensional fractionation emerges. The majority of these separation tasks are extremely difficult in a singular centrifugal field of a tubular bowl centrifuge. This is because similarities in the sedimentation rate occur frequently even for a broad range of non-spherical particles. Future approaches could include the addition of a second force field, whose parameters show a stronger correlation regarding the NPs elongation and flatness. Adjustments of the tube centrifuge geometry are also conceivable, which can be studied and optimized with the presented dynamic simulation tool focusing on the apparatus performance in multi-dimensional fractionation.

Nevertheless, added value by the presented TCM emerges from a facilitated communication on multi-dimensional grade efficiencies for tubular bowl centrifugation processes. Moreover, the proposed algorithm that computes 3D-PTDs is applicable to any discrete dataset of particle traits which represent physical properties of nanosuspensions. This includes statistical data acquired from, e.g., micro-CT and TEM image analysis. Consequently, the workflow enables a fast and flexible numerical assessment of separation processes for a variety of different nanosuspensions. Furthermore, advancements in multi-dimensional characterization of particulate products in separation technology will assist with future model validation attempts. Lastly, real-time computation speeds of the 3D-TCM highlight its potential use in model predictive control scenarios. The novelty here comes with the simultaneous consideration of two or more separation criteria in process optimization.