- Email: [email protected]
Classical models of secondary settlers revisited
16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process Systems Engineering W. Marquardt, C. Pantelides (Editors) © 2006 Published by Elsevier B.V.
677
Classical M o d e l s of S e c o n d a r y Settlers Revisited R. David a, A. Vande Wouwer a, p. Saucez a, and J.-L. Vasel b aService d'Automatique - Service de Mathrmatiques, Facult6 Polytechnique de Mons, 31 Boulevard Dolez, 7000 Mons, Belgium bService Assainissement et Environnement, Universit6 de Libge, 185 Avenue de Longwy, 6700 Arlon, Belgium Secondary settlers, which ensure the separation of the activated sludge from the treated effluent, are described by distributed parameter models taking temporal and spatial variations into account. Over the years, a number of mathematical models have been proposed, including Kynch, Tak~ics and Hamilton models. The objective of this study is threefold: (a) to highlight the influence of the model formulation on the numerical solution procedure, (b) to propose a modified expression of the settling velocity, which avoids sharp spatial variations reported in the literature, and (c) to develop a numerical solution procedure following a method of lines strategy rather than the classical "tanks-in-series" approach.
Keywords: Secondary settler, Distributed parameter model, Method of lines 1. I N T R O D U C T I O N The performance of the activated sludge process (in a wastewater treatment plant) strongly depends on the performance of the secondary settler, which separates the sludge from water by gravity sedimentation. For process optimization, it is therefore of great interest to understand the dynamic behavior of the sedimentation process. In the past several decades, several one-dimensional models of sedimentation have been discussed. One of the pioneering works in this area is due to Kynch [2], who derived a basic mass balance partial differential equation (PDE) from Stoke's law. However, this simple model appears to be quite challenging in terms of numerical solution due to the existence of steep moving fronts and sharp spatial transitions. In order to alleviate these problems, Tak~ics and coworkers [5] proposed a discretization of the settler into several layers (following a standard "tank-in-series" formulation), associated with constraints on the material fluxes between layers, and a mathematical formulation of the settling velocity in terms of the concentration of solid particles. Tak~ics model is certainly one of the most popular nowadays. Some variations of the PDE have been proposed, such as the introduction of spatial dispersion in order to smooth off sharp gradients [1]. Recent simulation studies [3], [7] highlight the difficulties associated with the numerical simulation and parameter estimation of these mathematical models. The purpose of this paper is to analyze the origin of the reported difficulties, and to develop a numerical simulation tool based on a slightly modified formulation of the settling velocity and a method of lines solution procedure of the mass balance PDE, as implemented in the MATLAB toolbox MATMOL [6]. This paper is organized as follows. The next section introduces the basic one-dimensional model of sedimentation, the mathematical formulation of the settling velocity due to Tak4cs, and the classical spatial discretization into n layers ("tanks-in-series" formulation). Section 3 presents numerical simulation results obtained with Kynch, Tak4cs and Hamilton models and discusses potential improvements.
678
R. David
Section 4 proposes a modification of the expression of the settling velocity, implying that the settling velocity vanishes at a finite value of the solid concentration (in contrast to Takfics law), and discusses the influence of this modification on the classical models. Section 5 develops a numerical solution procedure following a method of lines strategy. Finally, Section 6 draws some conclusions. 2. SECONDARY S E T T L E R M O D E L L I N G E ~ clear wat~r Qf A
l
f
l
• Qw_ q _ S - - ~v
~
Q':, ~- _ q .
.
The secondary settler is a tank which separates the treated effluent (Qw) from the activated sludge (Qs) thanks to the effect of gravity (see Fig. 1). For simplicity, a simple geometry is considered in this study (tank of constant section A and depth ZL, fed at the rate Qf at the level zf ).
Figure 1: Settler sheme
2.1. Mass balance equation The considered mathematical models are based on a mass balance PDE for the solid particles (in concentration C) settling down with a sedimentation (or settling) velocity 7"
OC Ot
O(Fs+Cq) =Cfqf6(z-zf), Oz
(1)
where F~ - C Vs is the sedimentation flux and q is the hydraulic velocity: q~ in the lower part of the settler (z >_ zf) and - q w in its upper part (z < zf); 6 is the Dirac impulse function.
2.2. Settling velocity A fundamental parameter in the model formulation is the settling velocity ~. Takdcs and coworkers [5] proposed the following double-exponential law: vs = max (0, min (v~, Vo (e -rh(c-Cmin) - e-rp(C-Cmi")))). The velocity first increases due to the gravity acceleration exerted on the solid particles (influence of rp), reaches a maximum (v~) and then decreases as the particles are hindered by the other ones (formation of the sludge, influence of rh). V0 is the theoretical maximum value of the velocity and Cminthe minimum concentration needed for settling.
2.3.
97
Tanks-in-series 97 formulation
The traditional "tanks-in-series" formulation amounts to a discretization of the mass balance PDE with a first-order finite volume method. The settler is spatially discretized into n layers, on which a mass balance is calculated. The numerical concentration values represented in the next figures correspond to the middle of each layer. The flux Fs,i between two adjacent layers i and i + 1 is equal to Vs,iCi, where vs,i = vs(Ci). The mass balance equations can be written as follows (see also [7]): First layer: Az (~1 = qw C2 - qw C1 - Fs,1 i-th layer above the f e e d level (2 < i < m - 1): Az Ci = qwG+l - qw G + Fs,i-1 - Fs,i Feed level (layer m) • Az Cm = q f C f - qwCm - qsCm + Fs,m-1 - Fs,m i-th layer below the f e e d level (m + 1 < i < n - 1): Az (?i = q s G - 1 - q s G + Fs,i-1- Fs,i Last layer (tank's bottom, layer n) • AZ Cn = qsCn-l - qsCn + Fs,n-1
e t al.
Classical Models of Secondary Settlers Revisited
679
3. N U M E R I C A L SIMULATION STUDY The system of mass balance ordinary differential equations (ODEs) can be solved using the solver ODE15s of MATLAB [4]. The numerical values of the parameters used in this study are taken from [7]. The initial conditions correspond to a settling tank filled with clear water, i.e. G"= 0 for i = 1, ...,n. 3.1. Model of Kynch Equation (1) is the PDE model which was originally proposed by G. Kynch in 1952 [2]. Numerical simulation results are presented in Fig. 2. A steep concentration front moves through the settler and, after some time, the predicted spatial concentration profile becomes uniform in the lower part of the settler, with the exception of a sharp spatial transition at the outlet boundary, as observed in [3]. This sharp spatial transition at the outlet boundary is not physical, but results from the discretization scheme. The shape of the spatial concentration profiles can be explained by considering the fourth "tanks-in-series" equation at steady-state (Ci- 0): Ci = Ci-i -~ Fs,i-l-rs,i. The concentration profile is uniform since the
qs
term Fsi-~-Fs,i vanishes in steady state (as the sedimentation flux is constant). At the outlet boundary,
qs
the last "tanks-in-series" equation at steady-state Cn - 0 ) gives: ( C n - C n - 1 ) -- Fs,n-i The observed qs discontinuity results from the second term (sedimentation flux from the previous layer). 3.2. Model of Takfics Tak~ics [5] introduced constraints on the gravitational flux Fs, as an ad-hoc procedure applicable to the "mixed-tanks-in-series" formulation, based on the comparison between Fs values for successive layers: - in the upper zone of the settler (above the feed level): Fs,i = min (Vs,iCi, Vs,i+1Ci+ 1) if Ci+ 1 > Ct, Fs,i = Vs,iCi if Ci+l _< Ct, where Ct is a threshold concentration (G = 3000 g/m3 in [7]); - in the lower zone of the settler (below the feed level): Fs,i = min(Vs,iCi, Vs,i+l Ci+l). The simulation results are more realistic (accumulation of sludge at the bottom of the settler), but there is an interplay between the model formulation (conditions on F~) and spatial discretization (number of layers), as shown in [7]. Consequently, the predicted concentration profiles change when the number of layers n is increased, and there is no convergence to the desired physical profiles, but to the concentration profiles predicted by Kynch model! (see Fig. 3). This makes the model of Tak~ics very delicate to use in practice as an ad-hoc choice of a number of layers and of the numerical values of the model parameters has to be made a priori.
3.3. Model of Hamilton Hamilton and coworkers [ 1] proposed to introduce a term accounting for material dispersion in Kynch model: -07 •c + O(Fs+Cq) ~z D ~o32C = C f q f t~(z -- Zf ). The "tanks-in-series" formulation is modified accordingly. The dispersion term introduces back-mixing effects, which smooth off the traveling waves and give realistic concentration profiles (see Fig. 4). To understand how this is translated numerically, consider {/qw+~)(--~)fs,i-l-fs,i the second"tanks-in-series" equation at steady-state" G-G+I kqw+2~ -~-Ci-1 qw+2-~ + qw+2-~" The value of Ci is influenced by C/+1 (as in Kynch model), but also by C/_ 1. In Kynch model, the coefficient of C/+I is equal to 1 while in Hamilton model the influence of neighboring layers is weighted by coefficients, whose sum is equal to 1.This explains the smooth aspect of the concentration profiles along the settler. Numerical results converge for increasing numbers of layers, in contrast to Tak~ics results.
3.4. Summary and discussion The mathematical simplicity of Kynch model contrasts with the difficulty associated with its numerical solution. A discretization of the spatial domain into n layers - following a standard "tanks-in-series" approach - yields unrealistic steady-state profiles, with a sharp spatial transition at the outlet bound-
R. D a v i d et al.
680
ary. The numerical procedure suggested by Takfics to alleviate this problem has the drawback that the constraints on the sedimentation flux are intimately related to spatial discretization. The concentration profiles change dramatically when the number of layers is increased, and ultimately turn back to the profiles predicted by Kynch model. The number of layers recommended by Takfics (n - 10) is insufficient to ensure convergence and leads to numerical simulation programs in which the issues of numerical accuracy (solution convergence) and model identification (estimation of the physical model parameters from experimental observations) cannot be separated (which should be the case). The model proposed by Hamilton introduces a dispersion term, which smooths off the spatial gradients and facilitates the numerical solution procedure (constraints on the sedimentation fluxes are no longer required), but this model does not include a clear definition of the outlet boundary condition, which are required for a second-order (in space) equation. These observations suggest the following: - the best way to clearly separate the model formulation and the numerical solution procedure would be to define an appropriate mass balance PDE and its associated initial and boundary conditions, and to solve the resulting set of equations using an advanced numerical procedure such as the method of lines; - the definition of an outlet boundary condition implies the consideration of a vanishing settling velocity, which would, following Takfics law, lead to infinite concentrations (indeed the second exponential in Takfics law vanishes for infinite concentrations). This latter point is addressed in the next section. 4. M O D I F I E D EXPRESSION OF THE SETTLING V E L OC IT Y 4.1. A n e w formulation
The original law implies that for Vs = 0, C is infinite. Physically, the sedimentation velocity should vanish at a finite concentration Csat, representing the concentration at which the sludge forms a solid aggregate. The second exponential in the original law is thus replaced by a parabola. Here, as an example value, Csat -- 11685g/m 3, which corresponds to the maximum value obtained with Kynch equation (see Fig. 2, t = 1000h) and is of the same order of magnitude as the maximum concentration reported in [5], and validated with experimental data. The first part of the sedimentation law is reduced to a simple exponential law with one parameter re (set to 1.494510 -3 m3//gin the study). The new sedimentation law is given by:
Vs = m a x ( 0 , m i n ( v r ,
Vo (e rp(C-Cmin) -
(C--Csat)4 (Cv~ -Csat) 4
VS : 4
1)))
for C < CvO for Cv(, -<- C __ < Csat
(2)
4v~
Vs = 0
for C > Csat
Cv~ is a concentration value corresponding to ~ and is the boundary value between the exponential law and the parabolic law. 4.2. I m p a c t on the previous m o d e l s
The results obtained with the models of Kynch and Hamilton are slightly modified, but remain basically the same. The biggest change is obtained with Takfics model, which now gives results comparable to those obtained with Hamilton's model. Furthermore, the numerical results now converge as expected with increasing n, see Fig. 5. A sharp, unrealistic transition is observed at the top of the settler though, which is the consequence of the discretization method and the condition on the solid flux. This anomaly is in fact also present in the numerical results obtained with Hamilton model, even though the introduction of the diffusion term softens the resulting curves. However, the saturated settler should not present this drastic reduction in concentration at the clear water outlet. Numerically, this reduction can be explained by the first Hamilton "tanks-in-series" equation at steady-state C1 --0),
Classical Models of Secondary Settlers Revisited
which shows that: C1 - - C 2 -
68 1
Fs~, , i.e., the first layer is only influenced by the next layer, and C] will
qw+~
always be smaller than C2.
5. METHOD OF LINES SOLUTION
The method of lines is a standard numerical solution procedure for PDE models, which proceeds in two separate steps: (ct) the spatial derivatives are first approximated using finite difference, finite element or finite volume methods; (15) the resulting semi-discrete (discrete in space - continuous in time) equations are then integrated in time using an ODE or DAE solver. The model of Hamilton, with the modified settling velocity, is solved using the Method of Lines. To this end, the settler is decomposed into two zones, I and H, which are delimited by the feed point. Indeed, the feed flow rate is a pointwise input, which introduces a discontinuity in the concentration profile. To avoid this spatial discontinuity, the source term in (1) is replaced, through the system decomposition into two zones, by a boundary condition. The mass balance equations for these two zones can be written as: OCt
O(-qwCl+Fs 1)
OCtl
O(qsCll+Fs,ll)
Ot + Ot +
Ozl
02CI : 0
' - D-~/-z21
Oz11
- D
02CII = 0
Boundary and continuity conditions are defined at the top, feed and bottom levels of the tank: oct -for z = 0:-~-z/=0 -forz=zf:
qfCf + (v~,i-qw)Ci-
(Vs,ii+qs)Cii -- 0
- for z - ZL: V~,II = 0 Initial conditions are given by: Q(t0,z) --0, and Cii(to,z) --O. Spatial grids are defined in each zone, with nl and n2 uniformly distributed grid points, respectively. The spatial derivatives are approximated using higher-order finite difference methods, implemented as differentiation matrices in a MATLAB library called MATMOL, [6]. In particular, the first-order spatial derivatives (convective terms) are approximated with five-point biased upwind FD schemes method, and the second-order spatial derivatives (dispersive terms) are computed with five-point centered schemes. The resulting system of differential-algebraic equations (DAEs) is integrated in time using the MATLAB solver ODE 15s. The resulting concentration profiles are physically coherent and no anomaly is observed. Convergence is ensured for small numbers of grid points (nl -- 76 and n2 -- 76 in Fig. 6), and the computational load is very reasonable (a few seconds on a standard PC).
6. CONCLUSIONS Most of the difficulties experienced in the numerical simulation of secondary settlers result from the interplay between model formulation (e.g. constraints on the sedimentation flux) and the discretization of the spatial domain into finite volumes (classical "tanks-in-series" approach). In the present study, this problem is circumvented by focusing attention on: (a) the definition of a system of well-posed mass balance PDEs and their associated initial and boundary conditions, (b) their solution using an advanced numerical procedure such as the method of lines. In addition, the influence of the formulation of the sedimentation velocity is highlighted. Particularly, at high (but finite) solid concentrations, the settling velocity should vanish, whereas classical laws only vanish asymptotically (i.e. at infinite concentration). The proposed model and simulation tool provide the dynamic evolution of the concentration profiles with accuracy and efficiency (a few seconds CPU on a standard PC).
682
R. D a v i d et al.
REFERENCES 1.
2. 3. 4. 5. 6.
7.
E (.D
J. Hamilton and R. Jain and P. Antoniou and S.A. Svoronos and B. Koopman and G. Lyberatos. Modeling and pilot-scale experimental verification for predenitrification process, J. Environ. Eng. - ASCE 118(1), 1992, pp 38-55. G.J. Kynch. A theory of sedimentation, Trans. Faraday Soc. 48, 1952, pp 166-176. I. Queinnec and D. Dochain. Modelling and simulation of the steady-state of secondary settlers in wastewater treatment plants, Wat. Sci. Tech. 43(7), 2001, pp 39-46. L.E Shampine, I. Gladwell and S. Thompson. Solving ODES with MATLAB, Cambridge University Press, 2003. I. Tak~ics, G.G. Patry and D. Nolasco. A dynamic model of the clarification-thickening process, Wat. Res. 25(10), 1991, pp 1263-1271. A. Vande Wouwer, P. Saucez and W.E. Schiesser. Simulation of Distributed Parameter Systems using a MATLAB-based Method of Lines Toolbox - Chemical Engineering Applications, Industrial Engineering and Chemistry Research, 43, 2004, pp 3469-3477. L.B. Verdickt and J. E Van Impe. Simulation analysis of a one-dimensional sedimentation model, Preprints of the 15 th triennial IFAC World Congress, Barcelona, Spain, 2002.
t
i:ii,,oii: !'::!j,,o;i
E
. . . . . . . . . . . . . . . nil ~l ~ +~~ l ,~. . . . . . . . . . . . . . . ¢"9
...............
) ~............
. . / 4 , ~ - ~ ,,.rs . . . . . . . . . . . . . . i-i ~,"~,o ....... =
/
/,
n
.... .... / . l
3o = ,z~, 0
-e,
z
[m]
z
[m]
z
[m]
z
[ m ]
Figure 2. Dynamic evolution of the concentration profile Figure 3. Changes in the steady-state concentration propredicted by Kynch model with n - 30. file with increasing values of n, as predicted by Tak~ics model.
, .8
"
0
1.8
,,
0
,"
" ¢,51-,
__~oooo (..)
(~
0
,.8
.
.
.
.
.
.
0
,8
,,
z [m] z [m] z [m] z [m] Figure 4. Dynamic evolution of the concentration profile Figure 5. Convergence of the steady-state concentration predicted by Hamilton model with n - 30. profile obtained with Tak~ics model for increasing n (from 10 to 100).
' °hi ¸¸¸¸¸¸¸¸¸
('~
o
18
,
o
,8
4
o
18
,
o
18
4
0
,8
4
o
18
,
,
o
,
o
o
,~
z
[m]
18
z
[rn]
,8
z
4
[m]
Figure 6. Dynamic evolution of the concentration profile obtained with the method of lines solution procedure (with nz -- 76 and n2 = 76).
677
Classical M o d e l s of S e c o n d a r y Settlers Revisited R. David a, A. Vande Wouwer a, p. Saucez a, and J.-L. Vasel b aService d'Automatique - Service de Mathrmatiques, Facult6 Polytechnique de Mons, 31 Boulevard Dolez, 7000 Mons, Belgium bService Assainissement et Environnement, Universit6 de Libge, 185 Avenue de Longwy, 6700 Arlon, Belgium Secondary settlers, which ensure the separation of the activated sludge from the treated effluent, are described by distributed parameter models taking temporal and spatial variations into account. Over the years, a number of mathematical models have been proposed, including Kynch, Tak~ics and Hamilton models. The objective of this study is threefold: (a) to highlight the influence of the model formulation on the numerical solution procedure, (b) to propose a modified expression of the settling velocity, which avoids sharp spatial variations reported in the literature, and (c) to develop a numerical solution procedure following a method of lines strategy rather than the classical "tanks-in-series" approach.
Keywords: Secondary settler, Distributed parameter model, Method of lines 1. I N T R O D U C T I O N The performance of the activated sludge process (in a wastewater treatment plant) strongly depends on the performance of the secondary settler, which separates the sludge from water by gravity sedimentation. For process optimization, it is therefore of great interest to understand the dynamic behavior of the sedimentation process. In the past several decades, several one-dimensional models of sedimentation have been discussed. One of the pioneering works in this area is due to Kynch [2], who derived a basic mass balance partial differential equation (PDE) from Stoke's law. However, this simple model appears to be quite challenging in terms of numerical solution due to the existence of steep moving fronts and sharp spatial transitions. In order to alleviate these problems, Tak~ics and coworkers [5] proposed a discretization of the settler into several layers (following a standard "tank-in-series" formulation), associated with constraints on the material fluxes between layers, and a mathematical formulation of the settling velocity in terms of the concentration of solid particles. Tak~ics model is certainly one of the most popular nowadays. Some variations of the PDE have been proposed, such as the introduction of spatial dispersion in order to smooth off sharp gradients [1]. Recent simulation studies [3], [7] highlight the difficulties associated with the numerical simulation and parameter estimation of these mathematical models. The purpose of this paper is to analyze the origin of the reported difficulties, and to develop a numerical simulation tool based on a slightly modified formulation of the settling velocity and a method of lines solution procedure of the mass balance PDE, as implemented in the MATLAB toolbox MATMOL [6]. This paper is organized as follows. The next section introduces the basic one-dimensional model of sedimentation, the mathematical formulation of the settling velocity due to Tak4cs, and the classical spatial discretization into n layers ("tanks-in-series" formulation). Section 3 presents numerical simulation results obtained with Kynch, Tak4cs and Hamilton models and discusses potential improvements.
678
R. David
Section 4 proposes a modification of the expression of the settling velocity, implying that the settling velocity vanishes at a finite value of the solid concentration (in contrast to Takfics law), and discusses the influence of this modification on the classical models. Section 5 develops a numerical solution procedure following a method of lines strategy. Finally, Section 6 draws some conclusions. 2. SECONDARY S E T T L E R M O D E L L I N G E ~ clear wat~r Qf A
l
f
l
• Qw_ q _ S - - ~v
~
Q':, ~- _ q .
.
The secondary settler is a tank which separates the treated effluent (Qw) from the activated sludge (Qs) thanks to the effect of gravity (see Fig. 1). For simplicity, a simple geometry is considered in this study (tank of constant section A and depth ZL, fed at the rate Qf at the level zf ).
Figure 1: Settler sheme
2.1. Mass balance equation The considered mathematical models are based on a mass balance PDE for the solid particles (in concentration C) settling down with a sedimentation (or settling) velocity 7"
OC Ot
O(Fs+Cq) =Cfqf6(z-zf), Oz
(1)
where F~ - C Vs is the sedimentation flux and q is the hydraulic velocity: q~ in the lower part of the settler (z >_ zf) and - q w in its upper part (z < zf); 6 is the Dirac impulse function.
2.2. Settling velocity A fundamental parameter in the model formulation is the settling velocity ~. Takdcs and coworkers [5] proposed the following double-exponential law: vs = max (0, min (v~, Vo (e -rh(c-Cmin) - e-rp(C-Cmi")))). The velocity first increases due to the gravity acceleration exerted on the solid particles (influence of rp), reaches a maximum (v~) and then decreases as the particles are hindered by the other ones (formation of the sludge, influence of rh). V0 is the theoretical maximum value of the velocity and Cminthe minimum concentration needed for settling.
2.3.
97
Tanks-in-series 97 formulation
The traditional "tanks-in-series" formulation amounts to a discretization of the mass balance PDE with a first-order finite volume method. The settler is spatially discretized into n layers, on which a mass balance is calculated. The numerical concentration values represented in the next figures correspond to the middle of each layer. The flux Fs,i between two adjacent layers i and i + 1 is equal to Vs,iCi, where vs,i = vs(Ci). The mass balance equations can be written as follows (see also [7]): First layer: Az (~1 = qw C2 - qw C1 - Fs,1 i-th layer above the f e e d level (2 < i < m - 1): Az Ci = qwG+l - qw G + Fs,i-1 - Fs,i Feed level (layer m) • Az Cm = q f C f - qwCm - qsCm + Fs,m-1 - Fs,m i-th layer below the f e e d level (m + 1 < i < n - 1): Az (?i = q s G - 1 - q s G + Fs,i-1- Fs,i Last layer (tank's bottom, layer n) • AZ Cn = qsCn-l - qsCn + Fs,n-1
e t al.
Classical Models of Secondary Settlers Revisited
679
3. N U M E R I C A L SIMULATION STUDY The system of mass balance ordinary differential equations (ODEs) can be solved using the solver ODE15s of MATLAB [4]. The numerical values of the parameters used in this study are taken from [7]. The initial conditions correspond to a settling tank filled with clear water, i.e. G"= 0 for i = 1, ...,n. 3.1. Model of Kynch Equation (1) is the PDE model which was originally proposed by G. Kynch in 1952 [2]. Numerical simulation results are presented in Fig. 2. A steep concentration front moves through the settler and, after some time, the predicted spatial concentration profile becomes uniform in the lower part of the settler, with the exception of a sharp spatial transition at the outlet boundary, as observed in [3]. This sharp spatial transition at the outlet boundary is not physical, but results from the discretization scheme. The shape of the spatial concentration profiles can be explained by considering the fourth "tanks-in-series" equation at steady-state (Ci- 0): Ci = Ci-i -~ Fs,i-l-rs,i. The concentration profile is uniform since the
qs
term Fsi-~-Fs,i vanishes in steady state (as the sedimentation flux is constant). At the outlet boundary,
qs
the last "tanks-in-series" equation at steady-state Cn - 0 ) gives: ( C n - C n - 1 ) -- Fs,n-i The observed qs discontinuity results from the second term (sedimentation flux from the previous layer). 3.2. Model of Takfics Tak~ics [5] introduced constraints on the gravitational flux Fs, as an ad-hoc procedure applicable to the "mixed-tanks-in-series" formulation, based on the comparison between Fs values for successive layers: - in the upper zone of the settler (above the feed level): Fs,i = min (Vs,iCi, Vs,i+1Ci+ 1) if Ci+ 1 > Ct, Fs,i = Vs,iCi if Ci+l _< Ct, where Ct is a threshold concentration (G = 3000 g/m3 in [7]); - in the lower zone of the settler (below the feed level): Fs,i = min(Vs,iCi, Vs,i+l Ci+l). The simulation results are more realistic (accumulation of sludge at the bottom of the settler), but there is an interplay between the model formulation (conditions on F~) and spatial discretization (number of layers), as shown in [7]. Consequently, the predicted concentration profiles change when the number of layers n is increased, and there is no convergence to the desired physical profiles, but to the concentration profiles predicted by Kynch model! (see Fig. 3). This makes the model of Tak~ics very delicate to use in practice as an ad-hoc choice of a number of layers and of the numerical values of the model parameters has to be made a priori.
3.3. Model of Hamilton Hamilton and coworkers [ 1] proposed to introduce a term accounting for material dispersion in Kynch model: -07 •c + O(Fs+Cq) ~z D ~o32C = C f q f t~(z -- Zf ). The "tanks-in-series" formulation is modified accordingly. The dispersion term introduces back-mixing effects, which smooth off the traveling waves and give realistic concentration profiles (see Fig. 4). To understand how this is translated numerically, consider {/qw+~)(--~)fs,i-l-fs,i the second"tanks-in-series" equation at steady-state" G-G+I kqw+2~ -~-Ci-1 qw+2-~ + qw+2-~" The value of Ci is influenced by C/+1 (as in Kynch model), but also by C/_ 1. In Kynch model, the coefficient of C/+I is equal to 1 while in Hamilton model the influence of neighboring layers is weighted by coefficients, whose sum is equal to 1.This explains the smooth aspect of the concentration profiles along the settler. Numerical results converge for increasing numbers of layers, in contrast to Tak~ics results.
3.4. Summary and discussion The mathematical simplicity of Kynch model contrasts with the difficulty associated with its numerical solution. A discretization of the spatial domain into n layers - following a standard "tanks-in-series" approach - yields unrealistic steady-state profiles, with a sharp spatial transition at the outlet bound-
R. D a v i d et al.
680
ary. The numerical procedure suggested by Takfics to alleviate this problem has the drawback that the constraints on the sedimentation flux are intimately related to spatial discretization. The concentration profiles change dramatically when the number of layers is increased, and ultimately turn back to the profiles predicted by Kynch model. The number of layers recommended by Takfics (n - 10) is insufficient to ensure convergence and leads to numerical simulation programs in which the issues of numerical accuracy (solution convergence) and model identification (estimation of the physical model parameters from experimental observations) cannot be separated (which should be the case). The model proposed by Hamilton introduces a dispersion term, which smooths off the spatial gradients and facilitates the numerical solution procedure (constraints on the sedimentation fluxes are no longer required), but this model does not include a clear definition of the outlet boundary condition, which are required for a second-order (in space) equation. These observations suggest the following: - the best way to clearly separate the model formulation and the numerical solution procedure would be to define an appropriate mass balance PDE and its associated initial and boundary conditions, and to solve the resulting set of equations using an advanced numerical procedure such as the method of lines; - the definition of an outlet boundary condition implies the consideration of a vanishing settling velocity, which would, following Takfics law, lead to infinite concentrations (indeed the second exponential in Takfics law vanishes for infinite concentrations). This latter point is addressed in the next section. 4. M O D I F I E D EXPRESSION OF THE SETTLING V E L OC IT Y 4.1. A n e w formulation
The original law implies that for Vs = 0, C is infinite. Physically, the sedimentation velocity should vanish at a finite concentration Csat, representing the concentration at which the sludge forms a solid aggregate. The second exponential in the original law is thus replaced by a parabola. Here, as an example value, Csat -- 11685g/m 3, which corresponds to the maximum value obtained with Kynch equation (see Fig. 2, t = 1000h) and is of the same order of magnitude as the maximum concentration reported in [5], and validated with experimental data. The first part of the sedimentation law is reduced to a simple exponential law with one parameter re (set to 1.494510 -3 m3//gin the study). The new sedimentation law is given by:
Vs = m a x ( 0 , m i n ( v r ,
Vo (e rp(C-Cmin) -
(C--Csat)4 (Cv~ -Csat) 4
VS : 4
1)))
for C < CvO for Cv(, -<- C __ < Csat
(2)
4v~
Vs = 0
for C > Csat
Cv~ is a concentration value corresponding to ~ and is the boundary value between the exponential law and the parabolic law. 4.2. I m p a c t on the previous m o d e l s
The results obtained with the models of Kynch and Hamilton are slightly modified, but remain basically the same. The biggest change is obtained with Takfics model, which now gives results comparable to those obtained with Hamilton's model. Furthermore, the numerical results now converge as expected with increasing n, see Fig. 5. A sharp, unrealistic transition is observed at the top of the settler though, which is the consequence of the discretization method and the condition on the solid flux. This anomaly is in fact also present in the numerical results obtained with Hamilton model, even though the introduction of the diffusion term softens the resulting curves. However, the saturated settler should not present this drastic reduction in concentration at the clear water outlet. Numerically, this reduction can be explained by the first Hamilton "tanks-in-series" equation at steady-state C1 --0),
Classical Models of Secondary Settlers Revisited
which shows that: C1 - - C 2 -
68 1
Fs~, , i.e., the first layer is only influenced by the next layer, and C] will
qw+~
always be smaller than C2.
5. METHOD OF LINES SOLUTION
The method of lines is a standard numerical solution procedure for PDE models, which proceeds in two separate steps: (ct) the spatial derivatives are first approximated using finite difference, finite element or finite volume methods; (15) the resulting semi-discrete (discrete in space - continuous in time) equations are then integrated in time using an ODE or DAE solver. The model of Hamilton, with the modified settling velocity, is solved using the Method of Lines. To this end, the settler is decomposed into two zones, I and H, which are delimited by the feed point. Indeed, the feed flow rate is a pointwise input, which introduces a discontinuity in the concentration profile. To avoid this spatial discontinuity, the source term in (1) is replaced, through the system decomposition into two zones, by a boundary condition. The mass balance equations for these two zones can be written as: OCt
O(-qwCl+Fs 1)
OCtl
O(qsCll+Fs,ll)
Ot + Ot +
Ozl
02CI : 0
' - D-~/-z21
Oz11
- D
02CII = 0
Boundary and continuity conditions are defined at the top, feed and bottom levels of the tank: oct -for z = 0:-~-z/=0 -forz=zf:
qfCf + (v~,i-qw)Ci-
(Vs,ii+qs)Cii -- 0
- for z - ZL: V~,II = 0 Initial conditions are given by: Q(t0,z) --0, and Cii(to,z) --O. Spatial grids are defined in each zone, with nl and n2 uniformly distributed grid points, respectively. The spatial derivatives are approximated using higher-order finite difference methods, implemented as differentiation matrices in a MATLAB library called MATMOL, [6]. In particular, the first-order spatial derivatives (convective terms) are approximated with five-point biased upwind FD schemes method, and the second-order spatial derivatives (dispersive terms) are computed with five-point centered schemes. The resulting system of differential-algebraic equations (DAEs) is integrated in time using the MATLAB solver ODE 15s. The resulting concentration profiles are physically coherent and no anomaly is observed. Convergence is ensured for small numbers of grid points (nl -- 76 and n2 -- 76 in Fig. 6), and the computational load is very reasonable (a few seconds on a standard PC).
6. CONCLUSIONS Most of the difficulties experienced in the numerical simulation of secondary settlers result from the interplay between model formulation (e.g. constraints on the sedimentation flux) and the discretization of the spatial domain into finite volumes (classical "tanks-in-series" approach). In the present study, this problem is circumvented by focusing attention on: (a) the definition of a system of well-posed mass balance PDEs and their associated initial and boundary conditions, (b) their solution using an advanced numerical procedure such as the method of lines. In addition, the influence of the formulation of the sedimentation velocity is highlighted. Particularly, at high (but finite) solid concentrations, the settling velocity should vanish, whereas classical laws only vanish asymptotically (i.e. at infinite concentration). The proposed model and simulation tool provide the dynamic evolution of the concentration profiles with accuracy and efficiency (a few seconds CPU on a standard PC).
682
R. D a v i d et al.
REFERENCES 1.
2. 3. 4. 5. 6.
7.
E (.D
J. Hamilton and R. Jain and P. Antoniou and S.A. Svoronos and B. Koopman and G. Lyberatos. Modeling and pilot-scale experimental verification for predenitrification process, J. Environ. Eng. - ASCE 118(1), 1992, pp 38-55. G.J. Kynch. A theory of sedimentation, Trans. Faraday Soc. 48, 1952, pp 166-176. I. Queinnec and D. Dochain. Modelling and simulation of the steady-state of secondary settlers in wastewater treatment plants, Wat. Sci. Tech. 43(7), 2001, pp 39-46. L.E Shampine, I. Gladwell and S. Thompson. Solving ODES with MATLAB, Cambridge University Press, 2003. I. Tak~ics, G.G. Patry and D. Nolasco. A dynamic model of the clarification-thickening process, Wat. Res. 25(10), 1991, pp 1263-1271. A. Vande Wouwer, P. Saucez and W.E. Schiesser. Simulation of Distributed Parameter Systems using a MATLAB-based Method of Lines Toolbox - Chemical Engineering Applications, Industrial Engineering and Chemistry Research, 43, 2004, pp 3469-3477. L.B. Verdickt and J. E Van Impe. Simulation analysis of a one-dimensional sedimentation model, Preprints of the 15 th triennial IFAC World Congress, Barcelona, Spain, 2002.
t
i:ii,,oii: !'::!j,,o;i
E
. . . . . . . . . . . . . . . nil ~l ~ +~~ l ,~. . . . . . . . . . . . . . . ¢"9
...............
) ~............
. . / 4 , ~ - ~ ,,.rs . . . . . . . . . . . . . . i-i ~,"~,o ....... =
/
/,
n
.... .... / . l
3o = ,z~, 0
-e,
z
[m]
z
[m]
z
[m]
z
[ m ]
Figure 2. Dynamic evolution of the concentration profile Figure 3. Changes in the steady-state concentration propredicted by Kynch model with n - 30. file with increasing values of n, as predicted by Tak~ics model.
, .8
"
0
1.8
,,
0
,"
" ¢,51-,
__~oooo (..)
(~
0
,.8
.
.
.
.
.
.
0
,8
,,
z [m] z [m] z [m] z [m] Figure 4. Dynamic evolution of the concentration profile Figure 5. Convergence of the steady-state concentration predicted by Hamilton model with n - 30. profile obtained with Tak~ics model for increasing n (from 10 to 100).
' °hi ¸¸¸¸¸¸¸¸¸
('~
o
18
,
o
,8
4
o
18
,
o
18
4
0
,8
4
o
18
,
,
o
,
o
o
,~
z
[m]
18
z
[rn]
,8
z
4
[m]
Figure 6. Dynamic evolution of the concentration profile obtained with the method of lines solution procedure (with nz -- 76 and n2 = 76).