The main purpose of the recently introduced Bürger–Diehl simulation model for secondary settling tanks was to resolve spatial discretization problems when both hindered settling and the phenomena of compression and dispersion are included. Straightforward time integration unfortunately means long computational times. The next step in the development is to introduce and investigate time-integration methods for more efficient simulations, but where other aspects such as implementation complexity and robustness are equally considered. This is done for batch settling simulations. The key findings are partly a new time-discretization method and partly its comparison with other specially tailored and standard methods. Several advantages and disadvantages for each method are given. One conclusion is that the new linearly implicit method is easier to implement than another one (semi-implicit method), but less efficient based on two types of batch sedimentation tests.

benchmark simulation, Bürger–Diehl simulation model, clarifier, compressive settling, linearly implicit time integration, numerical efficiency

Benchmark simulations of entire wastewater treatment plants (WWTPs) are today performed with one-dimensional simulation models of the secondary settling tank (Gernaey et al. 2014; Li & Stenstrom 2014). In the model by Bürger et al. (2011, 2013), sometimes referred to as ‘Bürger–Diehl model’, the physical phenomena of hindered settling, volumetric bulk flows, compression of the sludge at high concentrations and dispersion of the suspension near the feed inlet can be included in a flexible way. Each phenomenon is associated with a separate constitutive function with its model parameters, and can be activated or deactivated at the user's discretion. The possibility to include sludge compression is particularly important, since this improves the predictive power considerably (De Clercq et al. 2008; Ramin et al. 2014; Guyonvarch et al. 2015; Torfs et al. 2015; Van Loosdrecht et al. 2015). However, the inclusion of physical phenomena that result in second-order derivative terms in the model partial differential equation (PDE), e.g., compression, means that straightforward time discretization by easy-to-implement, explicit methods such as the explicit Euler method leads to long computational times.

Bürger et al. (2013) presented implementation details of a numerical algorithm, which gives reliable simulations with respect to the underlying physical principles and can be obtained with a user-defined accuracy. For long-time simulations of entire WWTPs, it is important to keep discretization errors small. This calls for a fine resolution in both space (many layers in the settler) and time (short time steps), which implies long computational times. Conversely, fast computations obtained with a low resolution in space and time come at the cost of poor accuracy.

In scientific computing, the numerical error is a measure of how close a numerical solution is to the exact solution or reference solution (obtained by a very fine discretization) of the model, i.e., the governing differential equation. The efficiency of a numerical method is assessed by relating the numerical errors to the computational (central processing unit; CPU) times necessary to obtain the numerical solutions for different discretizations. It is then said that one numerical method is more efficient than another if it allows a numerical solution with a determined numerical error to be obtained in less CPU time, or equivalently, a given budget of CPU time allows one to obtain a more accurate numerical solution by the first method than by the second.

The simulation model by Bürger et al. (2013) is based on a method-of-lines formulation of the underlying nonlinear PDE. This means a system of time-dependent ordinary differential equations (ODEs), one for each layer of the settler. Simulations of PDEs are stable and reliable if a so-called CFL (Courant–Friedrichs–Lewy) condition is satisfied. This gives a maximal time step for each given layer thickness ⁠. If only hindered settling and bulk flows are included, then the CFL condition means that can be chosen proportional to ⁠, i.e., ⁠, which results in fast simulations. When compression or dispersion is included and standard ODE solvers used, the CFL condition states that ⁠, which means very small when the error should be reduced (small is chosen).

The purpose of this contribution is to investigate different time-integration methods, of which one is new, with respect to efficiency and other aspects, such as implementation complexity. There is a qualitative difference in the numerical treatment depending on whether only hindered settling and bulk flows are considered or, in addition, compression or dispersion are included. For clarity of presentation, we limit ourselves here to batch sedimentation in a vessel with a constant cross-sectional area, and for which the depth z is measured from the suspension surface downwards to the bottom at ⁠. The model PDE is
formula
1
where the flux function contains the hindered settling velocity function and the compression function is ⁠, where K is a constant containing the solid and fluid mass densities and the acceleration of gravity, and is the effective solid stress function, which is zero below a critical concentration and increasing above.

The investigated methods are used for the simulation of two different batch settling tests in a vessel with m. While the present study is limited to batch settling in a column, the method proposed herein should be appropriate to handle secondary settling tanks as well (under suitable modifications that would mainly affect the convective flux, but not the compression term for whose discretization the present method has been tailored). It is therefore essential to demonstrate that the method can robustly handle several scenarios, in particular a scenario that exhibits some spatial variation of the solution. For this reason we consider a conventional Kynch test (KT; Kynch 1952), that is, the settling of an initially homogeneous suspension, and a Diehl test (DT; Diehl 2007), where a body of concentrated suspension is initially located above clear liquid, separated by, for example, a membrane. The DT shows a stronger variation of the solution than the KT. Both tests provide complementary information that can be used to identify large portions of the flux function f (Betancourt et al. 2014).

As hindered settling velocity function we choose for simplicity with m/h and m3/kg, and the effective solids stress chosen is
formula
where m2/s2 and kg/m3. To obtain a working numerical method an important preparation is to compute the primitive
formula
of the compression function ⁠. In the present case, with the functions and chosen herein, we can obtain in closed algebraic form; if this is not possible (for other choices of these functions), the function can be obtained numerically (Bürger et al. 2013). Now (1) can be written as
formula
For batch sedimentation, this PDE should be complemented with initial data and zero-flux boundary conditions
formula
2
Suitable numerical schemes for the approximate solution of (1) are based on subdividing the depth interval into a number N of layers of equal thickness ⁠. A discretization in space leads to the following method-of-lines formulation (Bürger et al. 2013), which is a system of N ODEs:
formula
3
formula
4
formula
5
where is the concentration of layer j, denotes the numerical convective flux due to hindered settling between layer j and layer ⁠, which can be chosen, for instance, as the Godunov flux, and is the compressive numerical flux chosen as
formula
Note that (2) implies that the boundary fluxes ⁠. Implementation details on how to compute these fluxes are provided by Bürger et al. (2013). The last term in (4) now becomes the usual second-order difference approximation of ⁠, i.e.,
formula
6
The method-of-lines formulation (3)–(5) is converted into a fully discrete scheme by time discretization. Although a variety of methods could be used for (3)–(5), there is nothing to gain by applying a standard ODE solver of accuracy higher than first order, e.g., a Runge–Kutta method. On the contrary, this will only cause longer computational times producing the same error as explicit Euler time stepping (Diehl et al. 2015).

For the discussion of the various time-stepping methods, which will take the numerical solution from to ⁠, it will be important to carefully distinguish between numerical fluxes and concentrations that are evaluated at the old time step and those evaluated at ⁠, which we mark by the respective upper index n and ⁠. The time-integration methods compared herein are the following, where are constants that depend on the choice of the functions f and d.

The explicit Euler method

The explicit Euler method was used for all simulations in Bürger et al. (2011, 2013). The CFL condition of the fully discrete scheme can be captured by when is small. The method is easy to implement since all terms in the right-hand side of (3)–(5) are evaluated at ⁠, and therefore each concentration is an explicit function of the known ones at the previous time ⁠. Thus, the fully discrete version of (4) is
formula
with analogous formulas replacing the boundary updates (3) and (5).

The semi-implicit method

The semi-implicit (SI) method (Bürger et al. 2006) is described in detail by Diehl et al. (2015). The CFL condition is ⁠. To advance the solution from to we must solve the following nonlinear system of algebraic equations, where ⁠, ⁠, are the unknowns:
formula
7
supplemented with analogous formulas replacing the boundary updates (3) and (5). These equations are solved iteratively, for example, by the Newton–Raphson method.

The linearly implicit method

The idea of the linearly implicit (LI) method goes back to Berger et al. (1979) and is based on first considering the contribution from the compressive flux term (6) from time to ⁠. The purpose of the LI method is to avoid the numerical solution of the nonlinear system of Equation (7). The nonlinearities in (7) are found in the evaluations of the function D. The idea is to replace in (7) by (where is a parameter connected to the convergence of the method) so that
formula
8
and to find a simple update formula for ⁠, ⁠, to be executed first in each time step. A stable implicit Euler time step implies the formula
formula
This means that a linear system of equations should be solved for at the next time step. Then (8) can be updated explicitly to obtain ⁠.
To completely describe the LI method, we define
formula
where is a (nominal) maximum concentration and is a parameter. This parameter can be chosen at the user's discretion. To provide some guidance we mention that based on theoretical considerations (not detailed here) and numerical evidence (presented below), the closer is chosen to the value 1, the more accurate is the scheme. However, the admissible time step behaves as when ⁠. Consequently, the choice of is subject to the competing goals of accuracy (small errors) and speed (short CPU times) of the simulation. A detailed discussion of the optimal choice of under these restrictions is provided below.

To advance the solution from to ⁠, one proceeds as follows:

  1. For ⁠, set ⁠.

  2. Solve the following linear system for ⁠:
    formula
    9
  3. Calculate from
    formula
    formula
    formula

Note that the linear system (9) can be written as follows, where ⁠:
formula
This tridiagonal linear system of equations can easily be solved, for example by the Thomas algorithm (cf., e.g., Diehl et al. 2015). A preliminary analysis on the stability of the scheme (not presented here) implies that the CFL condition of the LI scheme is
formula
10

Other time-stepping methods

Any ODE solver could be used for the method-of-lines system (3)–(5). Possibly competitive methods are adaptive step-size methods and we choose here the ODE solver ode15s in Matlab (2014), which is an implicit multi-step method of variable order with step-size control. Of several standard ODE solvers investigated by Diehl et al. (2015), this was the second most efficient one (after the SI method) in the investigations with both stand-alone settler simulations and benchmark simulations for the entire activated sludge process.

For both the KT and the DT we measure the performance of the numerical methods in terms of numerical error and CPU time. The error of a total simulated solution with N layers up to a time T is calculated by comparing with a reference solution obtained by Euler's method and layers ⁠. The relative error is calculated as
formula
For the KT we choose an initial concentration of and ⁠, which means that almost steady state has been reached. This test is also used to study the dependence of the numerical solutions produced by the LI method on the choice of the parameter ⁠. For each run, the time step is chosen according to the CFL condition (10) (with equality). Note that the right-hand side of (10) is a decreasing function of and it tends to infinity as approaches 1 from above. Figure 1 shows numerical solutions obtained for three different values of ⁠. This figure and Figure 2 indicate that for a given number of layers ⁠, the values of closest to 1 produce the solutions with the smallest errors, but the highest CPU times. We note that, for instance, for a given value of N the time step for can be chosen twice larger than that for ⁠, so the total CPU time for should be about half of that for ⁠. To assess which value of is optimal, an efficiency plot (of relative error versus CPU time) is helpful, see Figure 2. The curves are roughly L-shaped, indicating that should be chosen close to the point corresponding to the point of ‘bend’, since smaller or larger values would lead to larger numerical errors (at only slightly smaller CPU times) or to smaller numerical errors at much increased CPU times. Based on these considerations, the choice seems a good compromise for further comparison with other methods.
Figure 1
Figure 1. Kynch test simulated by the LI method with layers and the indicated values of .
View largeDownload slide

Kynch test simulated by the LI method with layers and the indicated values of ⁠.

Figure 1
Figure 1. Kynch test simulated by the LI method with layers and the indicated values of .
View largeDownload slide

Kynch test simulated by the LI method with layers and the indicated values of ⁠.

Close modal
Figure 2
Figure 2. LI method applied to KT: efficiency plots (relative error versus CPU time) for the indicated values of N, determined for (from top to bottom, the curve for is labeled for illustration).
View largeDownload slide

LI method applied to KT: efficiency plots (relative error versus CPU time) for the indicated values of N, determined for (from top to bottom, the curve for is labeled for illustration).

Figure 2
Figure 2. LI method applied to KT: efficiency plots (relative error versus CPU time) for the indicated values of N, determined for (from top to bottom, the curve for is labeled for illustration).
View largeDownload slide

LI method applied to KT: efficiency plots (relative error versus CPU time) for the indicated values of N, determined for (from top to bottom, the curve for is labeled for illustration).

Close modal
For the DT we choose the initial data
formula
Figure 3 shows the numerical solution by the LI method with and a contour plot of the reference solution. The right plot illustrates that in those regions where ⁠, therefore ⁠, and (1) reduces to a first-order hyperbolic PDE, iso-concentration curves are straight lines, in complete agreement with the corresponding theory.
Figure 3
Figure 3. Numerical solution of the DT. Left: solution by the LI method with and . Right: contour plot of the reference solution, showing areas of (light grey; blue in online version), (medium grey; light brown) and (dark grey; dark brown) and iso-concentration lines corresponding to the annotated values of C [kg/m3]. The full colour version of this figure is available online.
View largeDownload slide

Numerical solution of the DT. Left: solution by the LI method with and ⁠. Right: contour plot of the reference solution, showing areas of (light grey; blue in online version), (medium grey; light brown) and (dark grey; dark brown) and iso-concentration lines corresponding to the annotated values of C [kg/m3]. The full colour version of this figure is available online.

Figure 3
Figure 3. Numerical solution of the DT. Left: solution by the LI method with and . Right: contour plot of the reference solution, showing areas of (light grey; blue in online version), (medium grey; light brown) and (dark grey; dark brown) and iso-concentration lines corresponding to the annotated values of C [kg/m3]. The full colour version of this figure is available online.
View largeDownload slide

Numerical solution of the DT. Left: solution by the LI method with and ⁠. Right: contour plot of the reference solution, showing areas of (light grey; blue in online version), (medium grey; light brown) and (dark grey; dark brown) and iso-concentration lines corresponding to the annotated values of C [kg/m3]. The full colour version of this figure is available online.

Close modal
Figure 4 shows the efficiency curves for all methods investigated plus those produced by employing Matlab's ode15s solver. It turns out that the latter method is the least efficient, and that for both tests the implicit methods are most efficient, as expected from the corresponding CFL conditions. Moreover, for a given number of layers N, the LI method is the fastest, although not necessarily the one with the smallest error. Note that the efficiency curve for the SI method for the DT (the right plot of Figure 3) is composed of four symbols only. In fact, no information is available for the run with since the Newton–Raphson iterations did not converge for that case. While this situation could be easily overcome by the ad–hoc remedy of further reducing (say, to 80% of its maximal value determined by the CFL condition), it alerts us to a more fundamental problem observed with the SI method, namely that the convergence of (iterative) solvers for the nonlinear equations is not ensured a priori.
Figure 4
Figure 4. Efficiency plots for the KT (left) and the DT (right) and the Euler, SI and LI methods and Matlab's ode15s solver. In the right plot the cross corresponding to for the SI method is not shown since the Newton–Raphson iterations did not converge, i.e., the method failed.
View largeDownload slide

Efficiency plots for the KT (left) and the DT (right) and the Euler, SI and LI methods and Matlab's ode15s solver. In the right plot the cross corresponding to for the SI method is not shown since the Newton–Raphson iterations did not converge, i.e., the method failed.

Figure 4
Figure 4. Efficiency plots for the KT (left) and the DT (right) and the Euler, SI and LI methods and Matlab's ode15s solver. In the right plot the cross corresponding to for the SI method is not shown since the Newton–Raphson iterations did not converge, i.e., the method failed.
View largeDownload slide

Efficiency plots for the KT (left) and the DT (right) and the Euler, SI and LI methods and Matlab's ode15s solver. In the right plot the cross corresponding to for the SI method is not shown since the Newton–Raphson iterations did not converge, i.e., the method failed.

Close modal

For accurate simulation results, should be chosen small, whereby the CFL condition implies that the explicit Euler method requires much smaller time steps than the other methods. On the other hand, the SI method needs more computations at every time step and requires the evaluation of the Jacobian matrix of the nonlinear algebraic system of equations. However, even when this system is solved by the Newton–Raphson method, the SI method has turned out to be far more efficient than explicit Euler (Diehl et al. 2015). To reduce the implementation complexity and to remove the necessity to solve numerically a system of nonlinear algebraic equations at every iteration in the SI method, whose convergence depends on the Newton–Raphson iterations, the LI method can be used instead. The advantage of solving a linear system only in each iteration, which also means an easier implementation, is paid fort by the price of a larger numerical error.

With the aim of finding an efficient time-integration method for the simulation of sedimentation with compression, we have in this work confined ourselves to two batch sedimentation tests for the comparison of some methods. The ideas of the new LI method and its implementation have been provided. When evaluated on the batch settling tests, the following advantages and disadvantages can be concluded.

  • Explicit Euler method:

    • + Implementation easy.

    • + Convergence proof exists also for continuous sedimentation (numerical approximate solutions converge to the PDE solution as ⁠) by Bürger et al. (2005). Robust method.

    • − Not efficient.

  • SI method:

    • + Most efficient of the investigated methods. Under the assumption that the Newton–Raphson iterations find a solution each time step, convergence of the numerical method was proved for batch sedimentation by Bürger et al. (2006).

    • − Implementation more complex than Euler and LI. At each time step, a nonlinear system of algebraic equations is solved, e.g., by Newton–Raphson iterations, which require tolerance parameter to be set. There is no guarantee that these iterations converge.

  • LI method:

    • + Implementation easy. The ingredient in addition to the Euler method is basically that a linear system of equations is solved at each time step. Robust method.

    • ± Second most efficient for for batch sedimentation. The efficiency can be adjusted to some extent by a parameter. Fastest method for a given ⁠, but least accurate. Convergence proof in preparation.

  • Matlab's ode15s:

    • + Ready-to-use standard time solver.

    • − Implementation most complex of the investigated methods. Well-established robust ODE solver for stiff problems, but not developed for solutions containing discontinuities.

    • − Least efficient of the investigated methods.

RB is supported by Fondecyt project 1130154; Conicyt project Anillo ACT1118 (ANANUM); Red Doctoral REDOC.CTA, MINEDUC project UCO1202; BASAL project CMM, Universidad de Chile and CI2MA, Universidad de Concepción; and Fondef ID15I10291. RB and CM are also supported by CRHIAM, Proyecto Conicyt Fondap 15130015. SD acknowledges support from The ÅForsk Foundation, Sweden.

Berger
A. E.
Brézis
H.
Rogers
J. C. W.
1979
A numerical method for solving the problem ut − Δf(u) = 0
.
RAIRO Anal. Numér.
13
,
297
–
312
.
Google Scholar
Crossref
Search ADS
 
Betancourt
F.
Bürger
R.
Diehl
S.
Mejas
C.
2014
Advanced methods of flux identification for clarifier-thickener simulation models
.
Minerals Eng.
63
,
2
–
15
.
Google Scholar
Crossref
Search ADS
 
Bürger
R.
Karlsen
K. H.
Towers
J. D.
2005
A model of continuous sedimentation of flocculated suspensions in clarifier-thickener units
.
SIAM J. Appl. Math.
65
,
882
–
940
.
Google Scholar
Crossref
Search ADS
 
Bürger
R.
Diehl
S.
Farås
S.
Nopens
I.
Torfs
E.
2013
A consistent modelling methodology for secondary settling tanks: a reliable numerical method
.
Water Sci. Technol.
68
,
192
–
208
.
Google Scholar
Crossref
Search ADS
PubMed
 
Gernaey
K. V.
Jeppsson
U.
Vanrolleghem
P. A.
Copp
J. B.
2014
Benchmarking of Control Strategies for Wastewater Treatment Plants
.
IWA Scientific and Technical Report No. 23
.
IWA Publishing
,
London
,
UK
.
Kynch
G. J.
1952
A theory of sedimentation
.
Trans. Farad. Soc.
48
,
166
–
176
.
Google Scholar
Crossref
Search ADS
 
MATLAB
2014
Release 2014b
.
The MathWorks, Inc.
,
Natick, MA
,
USA
.
Ramin
E.
Wágner
D. S.
Yde
L.
Binning
P. J.
Rasmussen
M. R.
Mikkelsen
P. S.
Plósz
B. G.
2014
A new settling velocity model to describe secondary sedimentation
.
Water Res.
66
,
447
–
458
.
Google Scholar
Crossref
Search ADS
PubMed
 
Van Loosdrecht
M. C. M.
Lopez-Vazquez
C. M.
Meijer
S. C. F.
Hooijmans
C. M.
Brdjanovic
D.
2015
Twenty-five years of ASM1: past, present and future of wastewater treatment modelling
.
J. Hydroinformatics
17
,
697
–
718
.
Google Scholar
Crossref
Search ADS
 
© IWA Publishing 2016