Modeling the Batch Sedimentation of Calcium Carbonate Particles in Laboratory Experiments—A Systematic Approach

The design of continuous thickeners and clarifiers is commonly based on the solid flux theory. Batch sedimentation experiments conducted with solid concentrations still provide useful information for their application. The construction of models for the velocity of settling allows the estimation of the flux of solids throughout time, which can, in turn, be used to find the area of the units required to achieve a given solid concentration in the clarified stream. This paper addresses the numerical treatment of data obtained from batch sedimentation experiments of calcium carbonate particles. We propose a systematic framework to fit a model that is capable of representing the process features that involve (i) the numerical differentiation of data to generate initial estimates for the instantaneous velocity of settling; (ii) the integration of a differential equation to fit the model for the velocity of settling; and (iii) the assessment of the quality of the fit using common statistical indicators. The model used for demonstration has a theoretical basis combined with an empirical component to account for the effect of the particle concentrations and their state of aggregation. The values of the numerical parameters obtained are related to the characteristic dimensions of the aggregates and their mass-length fractal dimensions.


Motivation
The analysis of batch settling velocities provides useful information for the design of continuous gravity settlers, for the separation of particulate solids from a slurry. This practice has its roots in the classic development of sedimentation theory and continues to be used because batch settling data are easy to obtain from laboratory experiments and they provide lower bounds for sizing [1] (Chapter 5). The measurement of the height of the interface between the diluted and concentrated phases formed in controlled experiments in a glass cylinder allows for estimating the velocity of settling, which, in turn, can be used to determine the solid flux and the area necessary to safely operate the unit [2]. Typically, very fine particles settle slowly, while if aggregation occurs, larger velocities are expected. In practice, many factors that depend on the (i) particle properties; and (ii) flow characteristics influence the flow dynamics and the velocity of settling. Apart from those factors, the assumption that a particle settles freely without interference from other particles is inadequate for high-concentration suspensions. Often, the term free settling is applied to systems where the particles are able to settle individually, and hindered settling (or thickening) is used to describe the behavior at high particle concentrations. For many systems of practical interest, the sedimentation rates depend on the concentration and the state of aggregation of the particles, rather than the particle size alone [3]. The incorporation of aggregation in steady-state models describing the velocity of settling in a thickener was considered by Usher et al. [4] and later improved by Zhang et al. [5]. Our study focuses on the region of low solid concentrations of micrometer-sized calcium carbonate particles with irregular shapes, where the occurrence of aggregation is required to explain the behavior of the settling velocities observed experimentally.
After collecting their experimental data, several authors use least squares or weighted least squares regression to fit the height of the interface to a given model; see Torfs et al. [6] for an example. The fitted models range from purely empirical forms [7,8] to phenomenologicalbased forms [9,10]. Depending on their structure, the models can be complex to fit, since the height of the interface is often represented by one or more nonlinear ordinary differential equations (ODEs) or partial differential equations. A strategy to circumvent this difficulty consists in using numerical differentiation techniques to determine the instantaneous velocities of settling from consecutive experimental heights of the interface. This simplifies adjusting the velocities since the optimization problem is then described by nonlinear algebraic equations, thus avoiding the numerical integration of the ODEs. However, this numerical differentiation usually also introduces errors in the estimates of the velocity since the measurement errors are also differentiated; original data uncertainty can be drastically amplified by calculating the time derivatives of the data in the absence of adequate smoothing/regularization methods. This point is worthy of attention, and a numerical approach combining the advantages of both strategies might reduce the numerical issues associated with the model fitting while improving its accuracy. Practically, a first estimate of the velocity of settling obtained through the use of finite differences is easy to achieve. This might provide useful initial estimates for the model parameters, which can be later estimated more reliably through the direct integration of the model's ODEs. Furthermore, in the case of constrained least squares, the efficient estimation of the local derivatives of the differential model can be included in the optimization problem as well, making this procedure more efficient [11].
Another numerical aspect requiring further analysis is related to the assessment of the quality of the model [12]. Quite often, the parameters obtained are highly correlated, and this may lead to uncertain estimates of the velocity of settling when the model is used for prediction [13]. To characterize the confidence bounds of the values predicted by these models, the correct use of the confidence limits determined for the parameters is also crucial. However, much of the literature involving model validation of batch sedimentation data disregards these numerical issues. This paper intends to help fulfill this gap through the use of a systematic approach for treating the regressed data obtained from batch sedimentation experiments. The approach and the tools proposed herein are general concerning their application and can be applied to experiments with larger data volumes as well as to other systems, such as those involving sludge sedimentation. This paper contains four elements of novelty: (i) A model including the particle aggregation is used to describe the batch sedimentation of calcium carbonate particles in glass cylinders. The model stands on phenomenological principles but is aligned with the empirical representation of Richardson and Zaki [14] for the sedimentation of concentrated suspensions. (ii) The aggregation of particles is explicitly included in the model using fractal metrics. (iii) A numerical approach based on the combination of differential and integral methods commonly used to fit kinetic data is proposed. The velocity of settling is firstly fitted to find initial parameter estimates, and subsequently used to fit the height of the interface. (iv) Confidence intervals are constructed for the parameters and the quality of the fitting is assessed using common performance metrics, such as the parametric cross-correlation and the visual analysis of 95 % confidence ellipsoids for pairs of parameters.
The paper is organized as follows. Section 2 introduces the model used to describe batch sedimentation in glass cylinders. Section 3 describes the numerical approach used for model fitting. First, the architecture is presented; next, each step is detailed. In Section 4, we characterize the materials involved in the experiments and the experimental procedure. Section 5 presents the results obtained for a battery of tests carried out at different solid concentrations. Section 6 provides a final overview of the work and a summary of the results.

Particles Settling Model
In our notation, boldface lowercase letters represent vectors, boldface capital letters stand for continuous domains, blackboard bold capital letters are used to denote the discrete domains, and capital letters represent matrices. Finite sets containing ι elements are compactly represented by ι ≡ {1, · · · , ι}. The transpose operation of a matrix or vector is represented by " ".
First, we introduce the theoretical model for the sedimentation of uniformly sized particles in a column of liquid. Next, the model is extended to non-uniformly sized particles of different shapes. Let ϑ ϑ ϑ be the velocity vector field of the particles, assumed to have a characteristic diameter d p , falling through a Newtonian fluid under gravity. Their density is designated as ρ p , the density of the liquid is ρ l , and its viscosity µ l . The momentum balance representing the forces interacting with the particle is given by Navier-Stokes (see Equation (1)) [15] ρ l ∂ϑ ϑ ϑ ∂t where t stands for time, P for pressure, g for the vector field representing the gravity force, ∇ for the gradient operator, and ∇ 2 for the Laplacian. The first term on the left-hand side of (1) represents the moment accumulation per unit volume, the second is for the momentum loss by convection. The first term on the right-hand side represents the loss of momentum due to the pressure, the second the viscous force, and the third the gravitational force. At a steady state, the accumulation term is null. At low flow rates, corresponding to small Re( 1) numbers, the term ρ l ϑ ϑ ϑ · ∇ϑ ϑ ϑ µ l ∇ 2 ϑ ϑ ϑ and Equation (1) reduces to: When the inertia and viscous forces are in equilibrium, the net acceleration of the particles is zero and their velocity (in the y-direction) becomes constant. This implies that the viscous force also becomes constant and the bodies continue to move in this direction with constant or terminal velocity, commonly designated as the Stokes velocity, v St (see Equation (3)): In opposition, when the accumulation term in Equation (1) is considered, the dynamics of the velocity of the sphere, v(t), is obtained through the solution of the equation [16] dv dt As an illustration, when the particles are of micrometer size and the fluid is not significantly viscous, starting from a state of rest, the particles are able to achieve 99 % of v St in the interval of about 1 ms. This implies that the characteristics associated with the terminal velocity will be more important for the overall description of the system. Moreover, different behaviors are observed for multiple spheres falling in a liquid-filled tube due to their possible interactions with each other. These aspects, together with the interactions with the walls of the tube used, are not considered in the above model.
Happel and Brenner [17] analyzed the sedimentation of multiple spheres in a fluid, assuming that each falls in a spherical unit cell of liquid. The steady-state solution of Equation (2) for this conceptual framework has the form where v t is the terminal velocity of a single sphere in the population, also known as the settling velocity, is the volumetric fraction of fluid (or media porosity), and f (•) expresses a functional relation. This study also evidences that the settling velocities decrease with the concentration of suspensions, emphasizing the dependence between v t and the concentration of solids. A theoretical relationship between the velocity of settling and the concentration of suspensions was initially proposed by Maude and Whitmore [16]. For an overview of the different functions f (•), the reader is referred to Davies [18]; the most commonly used was developed by Richardson and Zaki [14], who proposed f ( ) = n , with n being a constant depending on the particles' Reynolds number. For systems where Re < 0.2, n was found to be equal to 4.65 [19]. The batch sedimentation experiments considered were performed in graduated cylinders ( Figure 1a). Initially (at time t 0 ), a perfectly mixed suspension with a concentration of solids c 0 = c(0) (here, expressed in g/L and corresponding to a volumetric concentration of c 0 /ρ p ) was added to a cylinder provided with a convenient height scale. The initial height of the suspension, h(t 0 ) = h 0 , was recorded. Then, a chronometer was switched on, and the reading of the height of the interface between the settling particles and the supernatant liquid, designated as h i = h(t i ), at a given set of time instants, t i , was carried out; see the scheme in Figure 1b. The test was extended until no relevant change was seen in the height of the interface. We designate this height by h ∞ and the corresponding time instant is t ∞ .  Above the interface, a clarified zone is formed where the concentration of solids is zero. The concentrated zone (shortened to Concentr. zone in Figure 1b) is also assumed to remain perfectly mixed, and a single characteristic diameter d p was previously considered, meaning that all particles settle at the same rate. In this case, a steady-state balance of the solids enables the characterization of their respective concentrations throughout the duration of the experiment [3]:

Laboratório de Processo de Separação
The porosity of the uniformly mixed suspension in the concentrated zone is, thus, also dependent on time, and is given by: Therefore, monitoring the position of the interface allows for estimating the instantaneous hindered settling velocity where the empirical relation of Richardson and Zaki [14] is used to account for the effect of the particle concentration. The parameter n can be estimated by fitting v(t) from h(t).
Clearly, even after assuming a single characteristic diameter d p , the possible effects of the particle shapes when irregular particles are used and the occurrence of particle aggregation are not included in this model. To account for the agglomeration of particles, we consider that the characteristic dimensions of the aggregates are their fractal dimensions, as in Meakin [20]. A common cause for aggregation is the presence of electric charges in the particles, resulting in the inclusion of ions, water, and other smaller particles in the aggregates that form through coagulation. Let d agg denote the diameter of the aggregates and D f denote its mass-length fractal dimension [21,22]. The effect of the formation of aggregates on the effective particle diameter in Equation (7) is explicitly introduced in the model using a new parameter, k. An additional parameter η is also needed to denote the fact that the final composition in the concentrated zone also includes a significant fraction of water; this means that the final value of (t ∞ ) is different from zero when the settling velocity reaches zero. In this analysis, the value of n in Equation (7) was kept fixed at 4.65, so that the remaining parameter values obtained in this study can be more directly compared with reference values determined using the relation postulated by Richardson and Zaki [19]. The resulting model is which coincidentally has a form similar to the model proposed by Michaels and Bolger [23]. It is worth mentioning that k and η are, respectively, equal to non-dimensional numbers (d agg /d p ) D f −1 and (d agg /d p ) 3−D f , considered in Nieto et al. [22] (Equation (2)). For compactness, let θ θ θ = (k, η) be the vector containing the parameters to be fitted. Each of them is constrained to a compact domain, i.e., where θ L i is the lower bound for parameter i, and θ U i is the upper bound. n θ is the number of parameters to fit in Equation (8), which in this context is 2.

Numerical Strategy for Model Fitting
In this section, we present the strategy used to fit Equation (7) to data resulting from laboratory batch experiments. Similar to the analysis of chemical kinetics rate laws, kinetic data can be analyzed employing the integral method [24,25] combined with robust numerical optimization procedures [26]. This method fits the kinetic rates expressed as a set of differential-algebraic equations to cumulative data (i.e., the reaction extent) calculated by sampling the reaction system during the experiment at a given set of time instants. This approach can be numerically challenging, as it requires a numerical ODE solver combined with an optimizer that handles the nonlinear least squares problem; an adequate initialization of the parameters is crucial to achieving the convergence of the optimizer and increasing the computational efficiency. Apart from this issue, this technique is rather appealing as it allows fitting the data as it was recorded in a bias-free adjustment (i.e., without any previous treatment)-see Schittkowski [27], Edsberg and Wedin [28] for applications.
Alternatively, one may use the differential method, which fits the kinetic rates to increments in each reaction advance [29]. Here, the model to be fitted is represented by algebraic rate equations, which can be linear or nonlinear, and naturally become easier to handle within the least squares procedure than the corresponding model in the integral method. However, the numerical differentiations of the sampled data are required prior to the fitting task. This operation is numerically unstable and amplifies the error of the derivatives with respect to the measurement errors of the acquired signals [30]. This approach was considered by Cremers and Hübler [31], among others. Recently, Vertis [32] proposed a scheme, combining both methods in the context of reaction network identification. The strategy inspired us to systematize the architecture of the model fitting tool used in this study, which appears in Figure 2. Now, we briefly analyze each step of the algorithm in Figure 2 and their connections; the next sections present an in-depth analysis of the tools involved. Let us consider that the data obtained from the experiments consist of records including the time instants at which the system was observed and the corresponding heights of the interface, i.e., D n o = {t i , h(t i ), i ∈ n o }, n o being the number of data points. These data are used to estimate the instantaneous hindered velocities using Lagrange interpolation (Step 2). The extended data set containing the original data and the estimates of the velocity of settling, i.e., D n o ,exp = {t i , h(t i ), v t (t i ), i ∈ n o }, are then used to create initial estimates of the model parameters with the differential method (see Step 3). Here, we use the least squares method to fit the algebraic form of the model and obtain starting parameter estimates, which are fed to Step 4, where the integral method is applied. Least squares is also used in Step 4 to regress the differential equation representing the model, and produce final estimates. The information produced here (i.e., the parameter estimates, the parametric covariance matrix, called S(θ θ θ est ), and the local prediction errors given by e = h − h est ) is fed to Step 5, where additional model accuracy statistics and metrics for evaluating the quality of fit are determined.

Numerical differentiation via 6 th order Lagrange polynomials
Model fitting: Differential method Model fitting: Integral method Step 4 Performance metrics: conf. intervals & correlation Step 5 Section 3.1 describes the technique used for constructing the estimates of the velocity of settling using Lagrange polynomial approximations. Section 3.2 states the mathematical formulations of the optimization problems solved to fit the data. Finally, Section 3.3 introduces the metrics used to analyze the quality of fit and the model's predictive accuracy.

Numerical Differentiation-Step 2
This section describes the approach used for estimating the velocity of settling at sampling time instants t i . For this purpose, we use Lagrange polynomials as they can easily cope with data measured at unevenly spaced points. Given an (ordered) set of k + 1 time instants represented by τ τ τ = {t 1 , · · · , t k+1 }, all distinct, at which the system was sampled, the k th degree Lagrange polynomial approximation for the height of the interface at t ∈ [t 1 , t k+1 ] is as follows: Then, we derive (9) with respect to t and construct local estimates for the velocity of the interface. The approximation produced is [33]: In our context, all the instants t of interest (those at which dh(t)/dt is estimated) belong to τ τ τ; consequently, t ∈ {t 1 , · · · , t k+1 }. To minimize the approximation error, we assume that t coincides with the middle point of the series of k + 1 points, corresponding to t (k+1)/2 +1 ; • is the floor operator. For the first and the last (k + 1)/2 points of D n o , t coincides with the point closer to the center of [t 1 , t k+1 ]. Equation (10) can be used with unequal time intervals t i − t i−1 .
In our application, we use Lagrange polynomials of the sixth degree. This means that in Equations (9) and (10), we have k = 6, the derivative in each point includes information of k + 1 = 7 contiguous points, and the middle point (the point at which the derivative is calculated) is the fourth of each set {t 1 , · · · , t 7 }.

Model Fitting-Steps 3 and 4
Here, we fit the models for the velocity of settling (Step 3 of the algorithm) and the height of the interface (in Step 4). Model fitting is an LS problem, where the sum of squared errors of the predictions is minimized. Nonlinear programming problems describing the model parametrization in Steps 3 and 4 are, respectively, and min θ θ θ Equations (11a) and (12a) are the objective functions in the LS procedure, (11b) and (12b) are the equalities representing the model, (11c) and (12d) are the boundedness conditions of the parameters and (12c) is the initial condition for the ODE representing the height of the interface. The problems are solved with the Levenberg-Marquardt algorithm; the details are in Nielsen and Madsen [34]. The problem (12) requires an ODE solver, and we use an algorithm of variable order and variable steps based on the numerical backward differentiation formulae. The numerical tolerance for the optimization was set to 1 × 10 −5 . Practically, to handle problem (12), we use a sequential approach, which sequentially iterates between the model solution and the likelihood minimization until convergence [35].

Quality of Fit Metrics-Step 5
We now present the performance indicators used to evaluate the quality of the fitting (see Step 5 of the algorithm).
To measure the predictive accuracy of the model, i.e., the ability to reproduce the data, we use the adjusted R-squared correlation coefficient, designated here as R 2 adj , as proposed by Seber and Wild [36]. To analyze the collinearity of the parameters, we use the covariance matrix S(θ θ θ est ) (given by the linearized model from the LS algorithm at convergence) to construct an approximation to the parametric correlation matrix. Let B(θ θ θ est ) be a diagonal (square) matrix of size n θ (=2) containing the square roots of the diagonal elements of S(θ θ θ est ), i.e., B i,i = S i,i , i ∈ n θ and B i,j = 0, i, j ∈ n θ , i = j. The correlation matrix of the parameter estimates is then given by Ideally, the off-diagonal elements of R(θ θ θ est ) should be as low as possible but several authors mention that 0.95 is a cut-off value above which the parameters cannot be estimated uniquely. In our context, the element R 1,2 of R(θ θ θ est ) holds special significance as it measures the interdependence between parameters k and η, and it is necessary for R 1,2 ∈ [−0.95, 0.95] to ensure the independence of the parameters from each other.
The 100(1 − α) % confidence level intervals for the parameters are as follows [36]: where F(α, n θ , n o + 1 − n θ ) is the value of the F distribution with (n θ , n o + 1 − n θ ) degrees of freedom at a 1 − α significance level and s is an approximation of the standard error of the prediction. The visual analysis of the 95 % confidence ellipsoids for pairs of parameters allows the identification of problems related to the capability of estimating the parameters uniquely. In practice, its orientation, which is a function of the elements of the parametric covariance matrix, is commonly used to diagnose parameter collinearity and local model identifiability issues [37].

Materials and Experimental Method
In this section, we characterize the material used in the sedimentation experiments and the experimental methodology.
All the experiments were carried out with calcium carbonate (CaCO 3 ) particles. In the experiments, we used commercial calcium carbonate of analytical grade (full specification: Calcium carbonate agr, reag.ph.eur., Labbox Export, Barcelona, Spain). The product has MW = 100.09 Da, bulk density equal to 0.760 kg/L, and meets the requirements for identifications A and B. The microscopy-based image of the particles is shown in Figure 3a. We note that the particles have irregular shapes. This effect is accounted for in the sedimentation model by parameter k. Figure 3b presents the particle size distribution (PSD) obtained by laser diffraction, using a Mastersizer 3000 (Malvern Instruments Co., Ltd., Malvern, UK) granulometer. The average particle diameter is 4.51 µm. The density of the calcium carbonate was measured by helium pycnometry in a MicroMeritics AccuPyc 1330 Pycnometer (Micromeritics, Norcross, GA, USA). The value obtained was ρ p = 2.5327 g/cm 3 . The water used in the experiments was previously distilled and was carried out at nearly 25°C. At this temperature, the density and viscosity of the water were ρ l = 0.99705 g/cm 3 and µ l = 8.891 × 10 −4 Pa s, respectively. For these conditions, v St = 1.91 × 10 −5 m/s. Previous to the experiments, a homogenized suspension with a given concentration, c 0 , in CaCO 3 , was produced. Then, a graduated cylinder was filled up to a given height with the suspension. There, the suspension was homogenized and the initial height, h 0 , was recorded. Then, the height of the interface was monitored throughout the experiment and its value was measured and recorded at each sampling time. During the initial stages, the frequency of the measurement is 1 min −1 , and in the later stages, when its variation is small, the frequency decreases to a lower value. The resolution of the scale used to measure the height is 1 mm and the internal diameter of the cylinder is 4.5 cm. For convenience, in the remaining sections, we adopt derived units of the I.S. to represent the height of the interface and the velocity of the interface. The former is quantified in cm (of the cylinder) and the second in cm/min.

Numerical Results
In this section, we present the results obtained for a sequence of experiments, where the concentration of CaCO 3 in the suspension varies. The values tested are 15 g/L, 20 g/L, 25 g/L, 30 g/L, and 35 g/L. Figure 4 compares the data obtained for c(0) = 25 g/L and the respective model predictions when the velocity of settling and the height of the interface are fitted. Figure 4a compares the estimated settling velocity obtained using the numerical differentiation technique introduced in Section 3.1 with the model predictions obtained by solving the problem (11). Figure 4b compares the experimental data obtained for the height of the interface and the model predictions obtained by solving the problem (12). In both cases, we present the data, the model predictions, and the respective 95 % confidence limits for the prediction. We note that the estimates based on v t (t) are less accurate than those based on h(t), as a result of the numerical differentiation introduced in Step 2 of the Algorithm. Consequently, the values of the settling velocity are somewhat dispersed and the 95 % confidence limits appear relatively large. Contrarily, the 95 % confidence limits for the height of the interface are remarkably tight, which is a consequence of the small error in the estimates. Figure 4b demonstrates that the model can adequately capture the main features of the system behavior.  Now, we compare the data with the model estimates for all solid concentrations tested. Figure 5a compares the predictions for the velocity of settling with data obtained via numerical differentiation and Figure 5b is for the height of the interface. The designations est and exp in the plots are for "estimated" and "experimental data", respectively. The model prediction in the differential method is worse than the prediction by the integral method, due to the uncertainty introduced in the numerical differentiation of data. This behavior indicates that the direct use of the differential form of the model (in the integral method) produces more accurate parameters, as expected. These parameters are obtained in the sense of maximum likelihood (free of bias). However, the successful application of the integral method depends on the availability of good initial parameter estimates to avoid local solutions and other convergence difficulties. Hence, the initial application of the differential method to provide reasonable estimates of the numerical derivatives seems essential for the reliable use of the integral method.  Table 1 presents the values of the parameters and their respective confidence intervals, obtained for all concentrations. Column 1 contains the solid concentrations of the suspensions tested. Columns 2 and 3 present the estimates for the parameters k and η, respectively, and their respective 95 % confidence ranges (after the ± sign). Column 4 presents the quality of the fit described by R 2 adj , column 5 is for the parametric cross-correlation between k and η, column 6 contains the values of d agg and column 7 contains the values of the mass-length fractal dimension D f . These two quantities are obtained by relating k and η with the non-dimensional numbers (d agg /d p ) D f −1 and (d agg /d p ) 3−D f (see Equation (13) below) introduced by Heath et al. [21] and Nieto et al. [22]: Table 1 shows that parameter k has a notorious decrease as the concentration of solids increases; this means that higher concentration suspensions produce more aggregates with average smaller sizes. Contrarily, the η values are almost independent of the initial concentration of solids. The agreement between the estimates and the data is remarkably high (with R 2 above 0.9980 for all the concentrations), which indicates the accuracy of the prediction model. The parametric cross-correlation is below 0.95 for all the concentrations, evidencing that the parameters are independent of each other, making the model locally identifiable. The average diameter of the aggregates is between 10 and 17 times larger than the equivalent dimension of the individual particles and tends to decrease as the concentration increases. Finally, D f is almost independent of the initial concentration of the suspension, as expected from a mass-length characteristic dimension when the same material is used. The values of D f are in good agreement with that obtained by Heath et al. [38] for calcite (the value indicated is 2.4).  Figure 6 compares the 95 % confidence regions for η versus k. We note that although the ellipse orientation indicates a certain degree of collinearity between the parameters, they can still be uniquely estimated.

Conclusions
We considered the problem of modeling the batch sedimentation of calcium carbonate particles from laboratory tests in graduated glass cylinders. Various experiments at different concentrations were carried out to analyze the effect of the initial concentration on the hindered settling velocity. The height of the interface between the supernatant zone and the clarified zone was monitored throughout the experiments. We proposed a systematic approach, introduced in Section 3, to fit a model to batch sedimentation data. This tool combines the differential and integral methodologies also used in the analysis of chemical kinetic data. Its structure involves: (i) The numerical differentiation of the height of the interface with a Lagrange-based polynomial approximation, to find local estimates of the velocity of settling, see Section 3.1; (ii) The fitting of the velocity of settling using the least squares algorithm with an integral approach, in Section 3.2; and (iii) The construction of numerical and graphical indicators to evaluate the model accuracy and the quality of fit, see Section 3.3.
The tool is demonstrated for the case of small-sized calcium carbonate particles, where the model used to describe the dynamics of the velocity of settling combines theoretical fundamentals with empirical knowledge. The phenomenological component follows from the balance of forces interacting in the particle; the empirical knowledge is used to account for the interactions between particles where two distinct effects are considered (i) the concentration of particles; and (ii) the aggregation of particles. This resulted in a simple but flexible two-parameter model, which can be related to other hindered settling velocity models in the literature. The tool proposed herein is general enough to find application in the treatment of batch sedimentation data for different systems, such as sludge, colloids, and metal oxide sedimentation, eventually with different or more complex kinetic models.
The error associated with the model structure is minimized by using a physically consistent representation that relies on (i) first principles-based knowledge; and (ii) empirical knowledge developed for similar systems. The observational error is taken into account to set the architecture behind the proposed tool. Specifically, our algorithm avoids fitting velocity data obtained by numerical differentiation, which may increase the influence of the observational noise; contrarily, we use the velocity data to provide a good parameter estimate that is then used as the initial solution to fit the original data to a differential model. The quality of fit metrics show that (i) the observational error allows identifying both model parameters independently (see the correlation matrix of the parameters), and (ii) the metrics capture the main features underlying the data (see the low magnitude of the standard error of the prediction). Finally, the numerical error is minimized by using accurate robust algorithms combined with good parameter initial estimates to fit the differential model.