Chemical Engineering Science 56 (2001) 45374553

[Link]/locate/ces
Phenomenological model of ltration processes:
1. Cake formation and expression
R. B urger
a,
, F. Concha
b
, K. H. Karlsen
c
a
Institute of Mathematics A, University of Stuttgart, Pfaenwaldring 57, D-70569 Stuttgart, Germany
b
Department of Metallurgical Engineering, University of Concepci on, Casilla 53-C, Correo 3, Concepci on, Chile
c
Department of Mathematics, University of Bergen, Johs. Brunsgt. 12, N-5008 Bergen, Norway
Received 26 September 2000; received in revised form 26 February 2001; accepted 27 March 2001
Abstract
The phenomenological theory of sedimentationconsolidation processes of occulated suspensions is extended to pressure
ltration processes. The local mass and linear momentum balances for the solid and liquid component together with appropriate
constitutive assumptions lead to a strongly degenerate (mixed hyperbolicparabolic) nonlinear partial dierential equation for the
local solids fraction, which together with initial and boundary conditions determines a dynamic cake ltration process. In the case
of a prescribed applied pressure function, we obtain a free boundary problem, in which the piston height has to be determined
simultaneously with the solids concentration. A numerical algorithm approximating the physically correct solution, with possible
discontinuities such as the cake}suspension interface, is presented and employed to simulate various cake ltration processes.
? 2001 Elsevier Science Ltd. All rights reserved.
Keywords: Filtration; Numerical analysis; Sedimentation; Simulation; Suspension; Free boundary problem
1. Introduction
Filtration is a process widely used in industry, espe-
cially in medicine and in chemical and mining as well
as in food and paper companies. The variety of equip-
ment and procedures that are utilized to lter the several
materials in these industries is very large. For example,
vacuum and pressure ltration is used in batch and con-
tinuous manner. Nevertheless, the basic principles under-
lying ltration in all its diversity are the same. Filtration,
as part of the more general eld of solidliquid separa-
tion, has its fundamentals in the ow of uids in a porous
medium and its quantication makes use, and abuse, of
the famous Darcy equation (Darcy, 1856).
The scientic foundation of ltration was accom-
plished in 1970s by several groups of researchers in
the Americas and in Europe. The development of the

Corresponding author. Tel.: +49-711-6857647; fax: +49-711-

6855599.
E-mail addresses: buerger@[Link] (R. B ur-
ger), fconcha@[Link] (F. Concha), [Link]@[Link]
(K. H. Karlsen).
theory of ltration has its parallel in the formulation and
application of the theory of thickening. For thickening,
the phenomenological foundation was laid in 1970s and
1980s (Concha & Bascur, 1977; Concha & Bustos, 1986;
Buscall & White, 1987; Auzerais, Jackson, & Russel,
1988), but the mathematical modelling and analysis and
the design and application of suitable numerical meth-
ods have been performed only recently (Bustos, Concha,
B urger, & Tory, 1999; B urger & Karlsen, 2001; B urger,
Evje, & Karlsen, 2000b). In ltration, the ow in a rigid
lter cake is well established and the case of a compress-
ible cake is approximated by interpreting appropriately
the same equations. For this reason, ltration theory deals
mainly with the conservation equations for the uid, pla-
cing less attention to the behaviour of the solid. If a the-
ory of ltration is aimed at, the deformation and motion
of the solid should be included.
A second, and most important issue, in ltration theory
is to recognize that the process is performed in stages,
which do not necessarily involve the same fundamental
principles. Take as an example the vacuum ltration in
a rotary lter. Four stages can be observed: (1) a cake
0009-2509/01/$ - see front matter ? 2001 Elsevier Science Ltd. All rights reserved.
PII: S0009- 2509( 01)00115- 4
4538 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
is formed by suctioning a suspension through a lter
medium that retains the solid and lets the uid pass; (2)
the cake is dewatered by suctioning air through the pores,
which replaces the liquid; (3) the cake is washed by suc-
tioning water through the cake and (4) the cake is dehu-
midied by again suctioning air. In the case of pressure
ltration, ve stages can be distinguished: (1) a cake is
formed by pressing a suspension through a lter medium;
(2) the lter cake is compressed, reducing its porosity,
by the direct contact with a solid surface, a mechanism
known as expression; (3) dewatering is accomplished by
blowing air through the cake; (4) water is pressed through
the cake to wash it and (5) dehumidication is obtained
by again blowing air through the cake.
In both cases, which represent the general trend in
industrial solidliquid ltration (see Massarani, 1997;
Wakeman & Tarleton, 1999; Rushton, Ward, & Holdich,
2000; Concha, 2001), the cake formation and the expres-
sion eliminate most of the liquid from the suspension,
leaving a saturated lter cake. These two steps are the
simplest of the whole process because they involve the
ow of only one uid, the liquid, through the porous body
formed by the consolidated particles. The principal vari-
ables are pressure, suspension concentration, temperature
and time. These stages of pressure ltration can be ana-
lyzed by using a simple theory of ow through a porous
medium. The dewatering and the washing stages in the
ltration cycle are more complex. They should be ana-
lyzed as drainage or embedding of a liquid in a porous
medium and, therefore, involve the ow of two uids, a
liquid and a gas, through the medium. In these stages, in
addition to the previous variables, the saturation and the
capillary pressure should be added. Saturation is the frac-
tion of liquid lling the pores of the porous medium and
capillary pressure is the pressure dierence established
between the air and the liquid when the rst replaces the
second from the lter cake. Capillary pressure and satu-
ration are microscopic phenomena occurring within the
porous medium and therefore, they rarely appear in the
macroscopic formulations of the ltration process we are
accustomed to. It is evident that the foundation of these
two last stages consists of the theory of two-phase ow
through a porous medium, which will be the topic of a
future paper.
The aim of this rst paper is to develop a theory for
the cake formation and expression in the ltration pro-
cess. It is organized as follows: In Section 2, we briey
specify the type of ltration process considered here. Sec-
tion 3 is devoted to the development of the mathematical
model. The local mass and linear momentum balances of
the solid and of one uid are recalled. Appropriate con-
stitutive equations for the solid and uid stresses and for
the soliduid interaction force are established. Bound-
ary conditions complement the resulting governing eld
equation for the solids concentration or equivalently, the
porosity. The entire formulation leads to the concept of
a dynamic ltration process, which from a mathematical
point of view can be regarded as an initial-boundary value
problem with a free boundary for one scalar strongly
degenerate (or mixed hyperbolicparabolic) convection
diusion equation. In Section 4 we rst mention some
of the particular features of this highly nonlinear partial
dierential equation, whose solutions can be found only
by numerical means. The mixed hyperbolicparabolic na-
ture of this equation implies that its solutions will typi-
cally be discontinuous and hence one has to be careful
when designing numerical methods for computing them.
Previous studies (Evje & Karlsen, 2000) have shown
that numerical methods based on naive nite dierence
discretization may fail to produce correct (i.e. physi-
cally relevant) solutions of strongly degenerate nonlinear
partial dierential equations. Asuitable (i.e. working) nu-
merical method for solving the mathematical model stud-
ied herein is introduced and numerical simulations of the
cake ltration process are shown. A discussion of these
results and an appropriate embedding of our model into
the existing literature is provided in Section 5.
2. Physical model
Consider the pressure ltration of a suspension con-
tained in the vessel with impervious and frictionless lat-
eral walls shown in Fig. 1. Pressure is exerted on the
suspension from above, for example by a piston mov-
ing downwards, while the uid is expelled from the bot-
tom through a semi-permeable membrane, called lter
medium. This medium can be a lter cloth or a membrane
which has the property of retaining the solid while let-
ting the uid ow through. During pressure ltration ve
steps can be distinguished: the cake formation, the ex-
pression, dewatering, washing and the dehumidication
of the lter cake by air blowing. The last three steps are
not considered in this paper because they involve vari-
ables (saturation and capillary pressure) not yet consid-
ered in this formulation.
In the sequel, let [ denote the volumetric solids con-
centration, which is assumed to depend on the height :
and on time t only. Before the ltration begins, the sus-
pension is assumed to have a homogeneous volumetric
solids concentration,
[(:, 0) = [
0
for 0 6: 6h
0
, (1)
where h
0
= h(0) is the initial height of the piston (see
Fig. 1(a)). During the step of cake formation two simul-
taneous processes may occur, the sedimentation of the
suspension and the ow through the porous cake and l-
ter medium at the bottom. In general, upto ve dierent
zones may be distinguished in the pressure vessel (see
Fig. 1(b)): At the top, a supernatant clear liquid zone
may appear with [ = 0. Below a transition zone where
the concentration varies between zero and [
0
, it may be
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4539
Fig. 1. Steps in the pressure ltration of a occulated suspension.
distinguished. A zone of initial concentration [=[
0
fol-
lows and a new transition zone with [
0
[ 6 [
c
may
exist underneath, where [
c
is the critical concentration at
which the solid ocs begin to touch each other. Finally,
there is a fth zone of consolidating porous cake with
concentration equal or higher than [
c
. Whether sedimen-
tation occurs or not depends, for a given material, on the
initial concentration [
0
and on the pressure exerted by the
piston (or the blowing air that may substitute the piston)
and consequently on the time involved on the ltration
process. The cake formation time is usually so short that
no signicant sedimentation occurs.
3. Mathematical model
Sedimentation, centrifugation, consolidation and l-
tration form part of the more general eld of solidliquid
separation processes of particulate systems. Therefore
each of them should be considered as a special case
of this general area of engineering science and their
governing equations should stem from the same root,
the mass and momentum balances for a particulate
system. For sedimentationconsolidation such models
have been presented by Concha, Bustos, and Barrientos
(1996, Chap. 3), B urger and Concha (1998), B urger
(2000) and B urger, Wendland, and Concha (2000c),
and for centrifugationconsolidation, by B urger and
Concha (2001). Filtration is the typical case where a
one-dimensional equation model of the ow in a porous
medium is adequate. Nevertheless, if sedimentation
ltration or centrifugationltration should be treated in
conjunction, it is necessary to start with the more general
one-dimensional model of a compressible particulate
process.
3.1. Dynamic process for a particulate system
Consider a body of solid particles intimately mixed
with a uid under the following conditions:
1. The solid particles are small with respect to the con-
tainer and of the same density.
2. The individual particles and the uid are incompress-
ible.
3. The particles are completely occulated at the begin-
ning of the process.
4. There is no mass transfer between the solid and the
uid.
5. The only body force is gravity.
6. The particles are contained in a vessel with impervious
frictionless lateral walls and a semi-permeable bottom.
Such a system may be regarded as continuous media with
two superimposed components obeying the following lo-
cal mass and linear momentum equations (Bustos et al.,
1999): the continuity equation of the solid,
c[
ct
+
c
c:
([t
s
) = 0, (2)
where t
s
denotes the solid phase velocity; the momentum
balance of the solid,
co
e
([)
c:
= pq[ +
m
d
1 [
, (3)
4540 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
where o
e
([) is the eective solid stress function, which
is a given constitutive function, p := p
s
p
f
0 is
the soliduid density dierence, q the acceleration of
gravity and m
d
the soliduid dynamic interaction force.
The mass balance equation of the mixture is
cq
c:
= 0, (4)
where q := [t
s
+(1[)t
f
t
s
(1[)t
r
is the volume
average velocity of the mixture if t
f
is the uid phase
velocity and t
r
:= t
s
t
f
denotes the soliduid relative
or drift velocity. The linear momentum balance equation
of the uid is
c
c:
= p
f
q
m
d
1 [
, (5)
where denotes the pore pressure. Eqs. (2)(5) follow
from the full mass and linear momentum balance equa-
tions for both components if advective acceleration terms
are neglected and if the solid and uid phase pressures
are expressed in terms of the pressures o
e
and . For de-
tails concerning the derivation of these equations, we re-
fer to Concha et al. (1996), B urger (2000) and B urger et
al. (2000c). If discontinuities appear, Eqs. (2)(5) have
to be replaced by jump conditions, which can be derived
from the macroscopic balances (Bustos et al., 1999).
3.2. Cake formation and expression
As mentioned above, the cake formation stage consists
of the simultaneous sedimentation of the suspension, con-
solidation of the sediment, and ow of ltrate through
the lter cake and the lter medium. Since the develop-
ment of each of these particular dierent operations has
followed dierent approaches, we will at rst respect the
original formulations, and will unify them afterwards.
3.2.1. Sedimentation of the suspension
In the theory of sedimentation it is customary to use
the Kynch batch ux density function
[
bk
([) :=
pq[
2
(1 [)
2
:([)
, (6)
where the resistance coecient :([) is assumed to be a
given function appearing in the constitutive relationship
m
d
= :([)t
r
. (7)
In view of Eq. (3), we obtain from (6)
[t
s
= q[ + [
bk
([)
_
1 +
o

e
([)
pq[
c[
c:
_
. (8)
Substituting this into (2) and taking into account (4),
we obtain the known strongly degenerate parabolic
equation describing the settling of occulated suspensions
(B urger & Concha, 1998):
c[
ct
+
c
c:
(q(t)[+[
bk
([)) =
c
c:
_

[
bk
([)o

e
([)
pq[
c[
c:
_
.
(9)
Introducing the diusion coecient a([) and its primitive
A([) by
a([) :=
[
bk
([)o

e
([)
pq[
, A([) :=
_
[
0
a(s) ds, (10)
Eq. (9) can be shortly rewritten as
c[
ct
+
c
c:
(q(t)[ + [
bk
([)) =
c
2
c:
2
A([). (11)
In addition, we obtain the following linear momentum
balance of the uid, which can be employed to calculate
the pore pressure a posteriori from the concentration
distribution:
c
c:
= p
f
q
pq[
2
(1 [)
[
bk
([)
t
r
. (12)
To illustrate that (9) is of strongly degenerate parabolic
type, we recall that it is usually assumed that o
e
([) van-
ishes for concentration values [ 6[
c
, and that
o

e
([) :=
d
d[
o
e
([) 0 for [[
c
. (13)
A common equation for o
e
that has been used to t ex-
perimental data for [[
c
is the power law
o
e
([) = o
0
(([}[
c
)
K
1) (14)
with parameters o
0
, K 0. Moreover, it is assumed in
general (Kynch, 1952) that the function [
bk
([) satises
[
bk
(0) = [
bk
([
max
) = 0,
[
bk
([) 0 for 0 [[
max
,
[

bk
(0) 0, [

bk
([
max
) 0,
(15)
where [
max
is the maximum solids concentration for the
material considered. A widely used ux density function
is the Michaels and Bolger modication (Michaels &
Bolger, 1962) of the RichardsonZaki ux density func-
tion (Richardson & Zaki, 1954)
[
bk
([) = u

[(1 ([}[
max
))
C
,
u

0, 0 [
max
61, C 1. (16)
Combining the properties (15) of [
bk
([) and those of
o
e
([), we see that the diusion coecient a([) given by
(10) vanishes for [ 6[
c
and [=[
max
and is positive for
[
c
[[
max
. In other words, Eq. (9) is hyperbolic for
[ 6[
c
and [ = [
max
and parabolic for [
c
[[
max
,
where the location of the type changes interface, at which
[ = [
c
is valid, is unknown beforehand. Since Eq. (9)
degenerates from parabolic to hyperbolic type on the
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4541
interval of solution values [0, [
c
] of positive length,
it is called strongly degenerate or mixed hyperbolic
parabolic.
3.2.2. Flow through the lter cake
To analyze a ow through a porous medium, the con-
cept of porosity c 1[ is frequently used instead of the
volumetric solids concentration. In terms of c, the linear
momentum equation of the uid (5) takes the form
c
c:
= p
f
q
m
d
c
. (17)
We assume that the following Darcy equation is valid,
which acts as a constitutive equation for the resistive
force:
m
d
=
j
f
k(c)
H(c)t
r
, (18)
where j
f
is the uid dynamic viscosity and k(c) is the
permeability of the porous medium. The function H(c)
can be chosen in such a way that the resulting linear
momentum equation for the very well known case of
ow through a rigid porous medium is satised. In this
case, we have t
s
0 and hence t
r
= t
f
and q = ct
f
.
Substituting (18) into (17) yields
c
c:
= p
f
q
j
f
k(c)
H(c)
c
2
q.
Since the classic result is
c
c:
= p
f
q
j
f
k(c)
q, (19)
we see that H(c) = c
2
. Consequently, from Eq. (5) we
have
m
d
=
j
f
c
2
k(c)
t
r
(20)
and therefore
c
c:
= p
f
q
j
f
c
k(c)
t
r
. (21)
In ltration, Eq. (21) is frequently referred to (Tiller &
Yeh, 1987; Tiller & Hsyung, 1993) as DarcyShirato
equation (Shirato, Sambuichi, Kato, & Aragaki, 1969).
It is usually stated in the form
c
e
c:
=
j
f
c
k(c)
t
r
, (22)
where

e
:= p
f
q(h(t) :) (23)
is the excess pore pressure, that is the pore pressure in
excess of the hydrostatic pressure.
Note that in view of t
r
= (t
s
q)}c, Eq. (21) can be
rearranged to yield
c
c:
= p
f
q +
j
f
c
k(c)
_
1
1 c
[ q
_
, (24)
where
[ = [t
s
(1 c)t
s
(25)
denotes the total solids ux density. Rearranging Eq. (24),
we obtain
q =
1
1 c
[
k(c)
j
f
c
_
c
c:
+ p
f
q
_
. (26)
Eq. (26) shows that the ltrate ow depends on the solid
ux density in addition to the pressure gradient.
3.2.3. Consistency between sedimentation and
ltration formulations
For the purpose of sedimentation, being in part respon-
sible for the formation of the lter cake, the linear mo-
mentum balance for the uid should be given by (12). In
the lter cake, where [[
c
is valid, this equation should
be equivalent to Eq. (24). Comparing the right-hand sides
of both equations, we see that the permeability function
k(c) is related to the Kynch batch ux density function
[
bk
([) by (B urger, Concha, & Tiller, 2000a)
[
bk
([) =
k(c)
j
f
pq[
2
=
k(1 [)
j
f
pq[
2
. (27)
Consequently, a constitutive equation k = k(c) for a
porous medium determines the portion of the Kynch
batch ux density function [
bk
([) for [[
c
. The in-
terrelation (27) is the key ingredient for the formulation
of a unied mathematical model for pressure ltration.
Provided that [
bk
([) satises (27), the eld equations
(9) and (12) provide all the information necessary for the
calculation of the ltration process. Moreover, inserting
t
r
=
[
bk
([)
[(1 [)
_
1 +
o

e
([)
pq[
c[
c:
_
, (28)
into Eq. (12) we obtain
c
c:
= p([)q
co
e
([)
c:
, (29)
where p([) = p
s
[ + p
f
(1 [) is the local density of
the mixture. Obviously, the pore pressure distribution
=(:, t) can be calculated a posteriori from the concen-
tration distribution. Consequently, the only eld equation
which actually has to be solved is Eq. (9). The unknown
ow rate q = q(t) will be determined from a boundary
condition, as shall be demonstrated below.
In view of this result, it should be mentioned that the
constitutive equation for o
e
given by (14) is equivalent to
4542 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
the frequent constitutive assumption (Tiller & Leu, 1980;
Stamatakis & Tien, 1991) that the local solids concentra-
tion in the cake is given by
[ = [
c
(1 + (o
e
}
A
))
[
(30)
with material specic parameters
A
and [ 0 if we
choose o
0
=
A
and k = 1}[. Eq. (30) can be used to
transform a similar approach for the local permeability
k = k
0
(1 + (o
e
}
A
))
o
(31)
with parameters k
0
, o 0 into values of [
bk
([), which
yields
[
bk
([) =
pqk
0
[
o}[
c
j
f
[
2(o}[)
. (32)
Under most circumstances, however, we have 2
(o}[) 0, so that the expression (32) is appropriate only
where [ is bounded away from zero, which is the case
under the restriction of formula (27) to [[
c
.
3.3. Boundary conditions and ltration rate
3.3.1. Kinematic boundary conditions
To obtain suitable boundary conditions, rst note that
(4) implies that
q = q(t) = h

(t) :=
dh(t)
dt
. (33)
In each mode of operation, the presence of the lter
medium at : =0 implies that no solids leave the ltration
equipment, i.e.
t
s
|
:=0
= 0. (34)
Recalling from (10) the denitions of a([) and A([), we
obtain from Eq. (8) the following boundary condition at
: = 0:
_
[
bk
([)
cA([)
c:
_
(0, t) = h

(t)[(0, t), t 0.
(35)
The second boundary condition has to be prescribed at
: = h(t). Therefore, we note that if the ltration rate is
faster than the settling velocity,
t
s
|
:=h(t)
= h

(t).
Inserting this into (8) yields
_
[
bk
([)
cA([)
c:
_
(h(t), t) = 0, t 0. (36)
3.3.2. Coupling between applied pressure and ltrate
ow rate
Combined with the initial concentration distribu-
tion given by (1), conditions (35) and (36) yield an
initial-boundary value problem (IBVP) for solving
Eq. (9). However, the IBVP (1), (9), (35), (36) is of
dierent character for both modes of operation.
To elucidate this argument, we now derive the cou-
pling equation between the piston height h(t) and the ap-
plied pressure o(t). To this end, rst note that the sum of
the momentum balance equations (3) and (5) yields the
equation
c
t
c:
= p([)q, (37)
where

t
(:, t) = (:, t) + o
e
([(:, t)) (38)
denotes the total pressure. Integrating (37) with respect
to : and bearing in mind that the total pressure at piston
height : = h(t) is the applied pressure, i.e.

t
(h(t), t) = o(t), (39)
we see that at : =0, the total pressure equals the applied
pressure plus the weight of the solid and the uid:

t
(0, t) =o(t) + q
_
h(t)
0
p([(, t)) d
=o(t) + q(m
0
p
f
(h
0
h(t))), (40)
where
m
0
=
_
h
0
0
p([
0
()) d
denotes the total mass of suspension per unit cross-
sectional area of the lter at t = 0. Inserting (40) into
(38), we obtain
o(t) = o
e
([(0, t)) q(m
0
+ p
f
(h(t) h
0
)) + (0, t).
(41)
3.3.3. Filter medium resistance and nal form of the
coupling equation
To determine the term (0, t) in Eq. (41), we consider
the ow of the uid through the lter medium. The thick-
ness of the lter medium, h
m
, is assumed to be very small.
Therefore we denote by 0
+
values of a variable just
above and by 0

just below the lter medium (Land-

man & Russel, 1993). Darcys law for the ow through
the lter medium at : = 0 gives
t
f
(0

, t) =
k
m
j
f
c
c:
(0

, t), (42)
where k
m
is the lter permeability. We may now approx-
imate the derivative on the right-hand side of (42) by a
dierence quotient in order to obtain
t
f
(0

, t) =
k
m
j
f
h
m
((0
+
, t)
atm
). (43)
We can assume that the pressures are scaled in such a
way that the atmospheric pressure
atm
is zero (Landman,
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4543
Sirako, & White, 1991). We assume that the space im-
mediately below the lter medium is entirely lled with
uid, i.e.
[(0

, t) = 0. (44)
To ensure continuity of the liquid ow, we require that
(1 [(0
+
, t))t
f
(0
+
, t) = (1 [(0

, t))t
f
(0

, t). (45)
In view of
(1 [(0
+
, t))t
f
(0
+
, t) = h

(t), (46)
we obtain from (44) to (46)
t
f
(0

, t) = h

(t). (47)
If we replace the ratio h
m
}k
m
by the medium resistance
R
m
(Tiller & Yeh, 1987; Font & Hern andez, 2000), we
obtain by inserting (47) into (43) the desired coupling
equation
o(t) =o
e
([(0, t)) q(m
0
+ p
f
(h(t) h
0
)) j
f
R
m
h

(t).
(48)
We nally mention that Eq. (48) coincides with the
formulations by both Landman and Russel (1993) and
Font and Hern andez (2000). Eq. (48) establishes the
relationship between the pressure applied to the piston,
o(t), and the ltrate ow rate or liquid expression rate
dh(t)}dt. Since Eq. (48) involves the unknown solution
value [(0, t) and is essentially based on a force balance,
we may refer to it as dynamic boundary condition.
3.4. Dynamic cake ltration process
Collecting the results, we can say that the sought vari-
ables, the concentration [=[(:, t) (or, equivalently, the
porosity c =c(:, t) in a saturated medium), the pore pres-
sure = (:, t) and one of the two variables applied
pressure o=o(t) or ltrate ow rate h

(t) q(t), are so-

lution of a dynamic cake ltration process if they satisfy
the eld equations
c[
ct
+
c
c:
(h

(t)[ + [
bk
([)) =
c
2
A([)
c:
2
,
0 : h(t), t 0, (49)
c
c:
= p([)q
co
e
([)
c:
, 0 : h(t), t 0, (50)
where the function A([) is dened in (10). These eld
equations are considered together with the initial condi-
tions
h(0) = h
0
, (51)
[(0, :) = [
0
(:), 0 6: 6h
0
, (52)
the kinematic boundary conditions
_
[
bk
([)
cA([)
c:
_
(0, t) = h

(t)[(0, t), t 0,
(53)
_
[
bk
([)
cA([)
c:
_
(h(t), t) = 0, t 0, (54)
and the dynamic boundary condition
o(t) =o
e
([(0, t))
q(m
0
+ p
f
(h(t) h
0
)) j
f
R
m
h

(t). (55)
The material properties of the suspension and the lter
cake are represented by the Kynch batch ux density
function [
bk
([), the eective solid stress function o
e
([)
and the constants p and j
f
. The ltration equipment is
characterized by the initial height h
0
and the lter medium
resistance. The control functions prescribed externally
are either the applied pressure o(t) or the ltrate ow
rate h

(t).
3.5. Modes of operation
There are principally two dierent ways to operate the
lter. Eqs. (53) and (55) play dierent roles in each of
these modes of operation.
1. The rate of ltrate is a specied function of time.
In this case, the signed ltrate ow rate dh(t)}dt is
specied at all times, and the applied pressure o(t)
has to be determined together with the solution of
Eq. (49).
In this case, Eqs. (49)(54) forman initial-boundary
value problem with a moving boundary : = h(t),
whose location is given a priori. This initial-boundary
value problem can be solved for the concentration [
in a straightforward fashion, and the applied stress
o(t) can be determined a posteriori by evaluating Eq.
(55) using the calculated values of [(0, t).
2. The applied piston pressure is a function of time. In
this more common case, the applied pressure o(t) is
prescribed (for example, is assumed to be constant or
is given by the pump ow rate), and the liquid ex-
pression rate dh(t)}dt is determined with the solution.
This case is more dicult since the unknowns [(0, t)
and h

(t) appear in a nonlinear fashion in both con-

ditions (53) and (55), and both quantities have to be
determined simultaneously.
Transforming (53) into an ordinary dierential
equation for h(t) and (55) into a boundary condition
for [(, t), we obtain
h

(t) =
1
[(0, t)
_
[
bk
([(0, t))
cA([)
c:
(0, t)
_
,
t 0 (56)
4544 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
and the Dirichlet boundary condition
[(0, t) = o
1
e
(o(t) + q(m
0
+ p
f
(h(t) h
0
))
+j
f
R
m
h

(t)), (57)
where o
1
e
is the inverse function of o
e
([) for
[[
c
. As mentioned above, the boundary h(t) de-
pends on o(t) and on the solution itself; therefore,
the initial-boundary value problem (49)(55) is a
problem with a free boundary, i.e. the upper bound-
ary h(t), the location of the piston at time t, is part
of the solution. The second case is, of course, more
dicult.
4. Numerical simulation of cake ltration processes
The equation describing the settling of occulated
suspensions (see (9)) is a highly non-standard par-
tial dierential equation. In particular, because of its
nonlinear nature and mixed hyperbolicparabolic type,
solutions may become discontinuous even if the initial
concentration is smooth. Consequently, the solutions
have to be interpreted in a weak sense. However, weak
solutions are not uniquely determined by their initial
and boundary data, and additional selection criteria, or
the so-called entropy conditions, are needed to single
out the unique physically relevant solution. We refer to
B urger and Karlsen (2001), B urger et al. (2000a, b) and
to the lecture notes by Espedal and Karlsen (2000) (as
well as the references cited in these papers) for details
on the mathematical background of strongly degenerate
nonlinear partial dierential equations.
The nonlinear nature of mathematical model derived
herein rules out analytical solution techniques and one
has no other choice than resorting to numerical methods.
When designing numerical methods for (9) one should be
a bit careful. It was observed in Evje and Karlsen (2000)
that numerical methods based on naive nite dierence
discretization may produce stable but physically wrong
solutions. This feature is intimately connected to the fact
that one has to impose the so-called entropy conditions in
the continuous model to single out its unique physically
relevant solution. Consequently, when designing numeri-
cal methods one has to make sure that the proposed meth-
ods has these entropy conditions or at least discrete
versions of them built in. As it is outside the scope of
this paper to explain in detail why the numerical method
described below produces correct solutions, let us simply
say that this is (roughly speaking) related to the conser-
vative form and the upwind nature of the method. In par-
ticular, the upwind nature implies that the method obeys
a sort of discrete entropy condition. We refer to B urger
and Karlsen (2001) and B urger et al. (2000a,b), Evje
and Karlsen (2000) and Karlsen and Risebro (2001) for
details concerning the mathematical justication of the
numerical method described below. For an overview of
numerical methods for general degenerate partial dier-
ential equations, we refer again to the lecture notes by
Espedal and Karlsen (2000).
4.1. Numerical algorithm
For the numerical treatment it is convenient to trans-
form the initial-boundary value problems to xed height.
To this end, we dene the new coordinate w = :}h(t).
Noting that
c[
ct
(:, t) +
c[
c:
(:, t)
d:
dt
=
c[
ct
(w, t)
wh

(t)
h(t)
c[
cw
(w, t), (58)
we obtain from (9)
c[
ct
+
h

(t)
h(t)
(1 w)
c[
cw
+
1
h(t)
c[
bk
([)
cw
=
1
(h(t))
2
c
2
A([)
cw
2
, 0 6w 61, t 0. (59)
The boundary conditions (35) and (36) take the respec-
tive forms
_
[
bk
([)
1
h(t)
cA([)
cw
_
(0, t) = h

(t)[(0, t),
t 0, (60)
_
[
bk
([)
1
h(t)
cA([)
cw
_
(1, t) = 0, t 0 (61)
and the initial condition (1) reads
[(w, 0) = [
0
(h(0)w, 0), 0 6w 61. (62)
It is now straightforward to adapt the upwind nite dif-
ference method developed for batch and continuous sedi-
mentation of occulated suspensions (B urger & Karlsen,
2001) to the initial-boundary value problem given by
Eqs. (59)(62). We refer to that paper (see also B urger
et al., 2000a,b) for details on the derivation of the method,
its properties, and on the underlying mathematical
background of degenerate hyperbolicparabolic partial
dierential equations, and state here the resulting nite
dierence formulas concisely. To this end, divide the
interval [0,1] into J subintervals [w
)
, w
)+1
] of length
w = 1}J, where w
)
:= )w, ) = 0, . . . , J.
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4545
Assume that an initial time step t
0
is given such that
the stability condition
1
h(t
n
)
_
|h

(t
n
)| + max
06[6[
max
|[

bk
([)|
_
t
n
w
+
2
(h(t
n
))
2
max
06[6[
max
|a(u)|
t
n
w
2
61 c,
n = 0, 1, 2, . . . , c 0 (63)
is satised for n = 0 with t
0
:= 0.
The calculation starts by setting
[
0
)
:=[
0
(h(0)w
)
), ) = 0, . . . , J, h
0
:= h
0
, h
1
:= h
0
,
(64)
where [
n
)
denotes the approximate value of [at w=w
)
and
t =t
n
, and similarly h
n
is the approximate location of the
free boundary, i.e. the piston height, at time t
n
. Assume
now that approximate solution values [
n
)
, ) = 0, . . . , J
and a value h
n
for time t
n
are known. If t
n
1, the
calculation is stopped. Otherwise, a new time step t
n
is
determined from Eq. (63). We set t
n+1
:= t
n
+ t
n
and
calculate [
n+1
)
, ) = 0, . . . , J, and a new height h
n+1
from
the following formulas:
The solution value [
n+1
0
is directly prescribed by setting
[
n+1
0
=o
1
e
(o(t
n+1
) + q(m
0
+ p
f
(h
n
h
0
))
+j
f
R
m
(h
n
h
n1
)). (65)
Calculate a new height h
n+1
from
h
n+1
=h
n

2t
n+1
[
n
0
+ [
n
1
_
[
EO
bk
([
n
0
, [
n
1
)

1
t
n+1
(A([
n
1
) A([
n
0
))
_
, (66)
where the function [
EO
bk
(, ) is the numerical ux func-
tion used by the scheme, namely the EngquistOsher
ux function (Engquist & Osher, 1981) dened by
[
EO
bk
(u, t) :=[
bk
(0) +
_
u
0
max{0, [

bk
(s)} ds
+
_
t
0
min{0, [

bk
(s)} ds. (67)
To ensure second order accuracy in space of the
method, we use the extrapolated variables
[
n, L
)
:= [
n
)

w
2
s
n
)
, [
n, R
)
:= [
n
)
+
w
2
s
n
)
,
) = 0, . . . , J, (68)
where the slopes s
n
)
are given by s
n
0
=s
n
1
=s
n
J1
=s
n
J
=0
and
s
n
)
:=
1
w
MM([
n
)
[
n
)1
, ([
n
)+1
[
n
)
)}2, [
n
)+1
[
n
)
),
) = 2, . . . , J 2, (69)
where the so-called minmod limiter function MM(, , )
is given by
MM(a, b, c) =
_

_
min{a, b, c} if a, b, c 0,
max{a, b, c} if a, b, c 0,
0 otherwise.
(70)
The values [
n+1
)
, ) = 2, . . . , J 1, are then calculated
from the following interior scheme, which is a conser-
vative nite dierence discretization of Eq. (59):
[
n+1
)
=
h
n
h
n+1
_
[
n
)

t
n+1
h
n+1
w
_
h
n+1
h
n
Jt
n+1
((J () + 1))[
n, L
)+1
(J ))[
n, L
)
)
+[
EO
bk
([
n, R
)
, [
n, L
)+1
) [
EO
bk
([
n, R
)1
, [
n, L
)
)
_
+
t
n+1
(h
n+1
)
2
w
2
[A([
n
)+1
) 2A([
n
)
)
+A([
n
)1
)]
_
, ) = 2, . . . , J 1. (71)
Inserting the appropriate discrete analogue of boundary
condition (60) into formula (71) for ) = 1, we obtain
the update formula of [
n
1
:
[
n+1
1
=
h
n
h
n+1
_
[
n
)

t
n+1
h
n+1
w

_
h
n+1
h
n
Jt
n+1
((J 2)[
n, L
2
(J 1)[
n, L
1
)
+[
EO
bk
([
n, R
1
, [
n, L
2
) [
EO
bk
([
n, R
0
, [
n, L
1
)
_
+
t
n+1
(h
n+1
)
2
w
2
[A([
n
2
) A([
n
1
)]
_
.
Similarly, the appropriate discrete version of boundary
condition (61) turns formula (71) for ) = J into an
update formula for [
n
J
:
[
n+1
J
=
h
n
h
n+1
_
[
n
J
+
t
n+1
h
n+1
w
[
EO
bk
([
n, R
J1
, [
n, L
J
)

t
n+1
(h
n+1
)
2
w
2
[A([
n
J
) A([
n
J1
)]
_
.
4.2. Numerical examples
The numerical algorithm was employed to simu-
late several dynamic cake ltration processes of a
4546 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
Fig. 2. The constitutive functions [
bk
(left) and o
e
(right) of the model suspension considered in the numerical examples.
hypothetical suspension represented by the constitutive
functions (see Fig. 2)
[
bk
([) =
_

_
0 if [ 60,
[
1
bk
([) := 5.9303 10
4
[(1 [}0.3)
22
m}s if 0 [ 60.03,
S
3
([) if 0.03 6[ 60.21,
[
2
bk
([) := 4.9758 10
12
[
5.

1
m}s if [0.21,
(72)
o
e
([) =
_
_
_
0 if [ 6[
c
:= 0.25,
19 000(([}[
c
)
11.

1
1)Pa if [[
c
,
(73)
where S
3
is a piecewise cubic spline function connecting
the segments [
1
bk
and [
2
bk
, whose properties are discussed
below.
The remaining parameters assumed to be con-
stant in each example are q = 9.81 m}s
2
, p =
1618.18 kg}m
3
, j
f
= 0.000978 Pa s, h
0
= 0.094 m and
R
m
= 3 10
10
m
1
. The values of h
0
, j
f
and R
m
have
been chosen in such a way that results can be compared
with an experiment of pressure ltration of a kaolin sus-
pension of initial concentration [
0
= 0.03 (see Fig. 8 of
Tiller, Hsyung, & Cong, 1995).
The rst segment [
1
bk
of [
bk
has been chosen accord-
ing to the Michaels and Bolger (1962) formula (16). The
coecient u

, the parameter [
max
and the exponent C
have been chosen in such a way that the settling rate ob-
served with a kaolin suspension at the dilute initial con-
centration 0.03 (Tiller et al., 1995) is well reproduced.
Moreover, in order to produce a lower rarefaction wave
as a smooth transition of concentrations (with respect
to height : at xed time t) between 0.03 and (roughly)
0.05, as suggested by the measurements given by
Fig. 8 of Tiller et al. (1995), these parameters were chosen
such that the local minimum of [
1
bk
is located between
zero and 0.03.
The actual structure of the rarefaction wave depends
on the spline segment S
3
, which in turn is determined by
the choice of a nite number of control points. Similar
constructions of portions of the Kynch batch ux density
function [
bk
, aimed at producing rarefaction waves in the
numerical solution of the corresponding nonlinear partial
dierential equation for [, were performed, for example,
by Diplas and Papanicolaou (1997) (see also Garrido,
Concha, & B urger, 2001) for data obtained from an atta-
pulgite settling experiment by Tiller and Khatib (1984).
The functions [
2
bk
and o
e
have been obtained from
Eqs. (30)(32) using the parameters
A
=19000 Pa, [=
0.09, [
c
=0.25, k
0
=5.85910
15
m
2
and o=0.64. Ex-
cept for the values of k
0
and [
c
, these parameters were de-
termined by Tiller and Horng (1983) for kaolin (see also
Tiller & Kwon, 1998; Tiller, Lu, Kwon, & Lee, 1999).
Our (corrected) choice of [
c
is based on the observation
that the lter cake solids concentrations reported in the
cited gure of Tiller et al. (1995) are lower than 0.32.
Using [
c
= 0.25, we rst performed a calculation with
the original value

k
0
=1.0910
14
m
2
published in these
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4547
references. The result was that the simulated ltration
rate was roughly twice as large as in the cited experiment
of Tiller et al. (1995). Therefore, we reduced the perme-
ability of the sediment at given porosity, and thereby the
ltration rate, by multiplying

k
0
with a constant smaller
than one. We found that the simulation reproduced the
evolution of the measured piston height satisfactorily by
choosing k
0
=0.538

k
0
. It should be pointed out that Tiller
et al. (1995) did not claim explicitly that the kaolin of
their experiment was the same as that examined previ-
ously by Tiller and Horng (1983). Since published con-
stitutive parameters for kaolin vary over large intervals
(see e.g. Tiller, Crump, & Ville, 1980; Wells, 1999), the
constitutive equations we nally utilized for this materi-
als are not unreasonable.
We mention that in absence of experimental data for
that range, the choice of the spline segment S
3
between
[=0.05 and 0.21 is essentially based on the consideration
that no new extremum of the function [
bk
should be
introduced, and that the number of inection points of
the resulting function [
bk
should be kept at the minimal
number of three. The shape of the ux density function
depicted in Fig. 2 is moreover very similar to ux density
functions for occulated suspensions determined by Scott
(1968) and Font, Garca, and P erez (1998).
The formation of and ow through the lter cake and
the evolution of the free boundary h(t) depend essentially
on the portion of [
bk
for [[
c
, so results obtained from
choosing dierent interpolations instead of that given by
S
3
would look very similar to those presented below.
The function A([), as a primitive of a([), used in the
calculations was determined numerically.
Our simulation of the experiment of Tiller et al. (1995),
where [
0
= 0.03, is shown in Fig. 3. Fig. 3(a) displays
a ltration plot (similar to a settling plot) with mea-
sured and simulated iso-concentration lines, while Fig.
3(b) shows the same numerical result as concentration
proles at selected times, together with concentrations
at

h(t), which are plotted as fat dots in time intervals
of length t. As a consequence of the chosen parame-
ters, the phenomenological model correctly predicts the
movement of the piston. The formation and height of
the lter cake are also well reproduced. The values of
the iso-concentration lines have been choosen in such a
way that the formation of a rarefaction wave, which in
part separated the bulk suspension of initial concentration
[
0
from the rising lter cake, where [[
c
, becomes
visible.
One could attempt to further modify the ux density
function [
bk
in such a manner that a concentration gra-
dient forms between the piston (or a possible clear liq-
uid supernate) and the bulk suspension, as suggested by
the concept of upper rarefaction waves (B urger & Tory,
2000). In other words, the main cause of the remaining
essential discrepancies between the measured data and
the numerical simulation lies in the choice of the ux
Fig. 3. Simulation of a ltration experiment by Tiller et al.
(1995): (a) ltration plot showing measured (symbols) and simu-
lated iso-concentration lines corresponding to the annotated values,
(b) simulated concentration proles at selected times and concentra-
tions (dots) at

h(t), plotted in time intervals of length t.
density function [
bk
in the dilute region 0 6 [ 6 0.03,
but not in the advanced theory of cake ltration.
In the remaining three examples, we consider some hy-
pothetical experiments by varying the initial concentra-
tion [
0
and the applied pressure function o(t), but leave
all other parameters unchanged.
We rst choose the initial concentration [
0
= 0.09.
Fig. 4(a) shows the simulated iso-concentration lines, and
Fig. 4(b) displays the corresponding concentration pro-
les. In addition, and following the presentation of nu-
merical results by Tiller and Yeh (1987), the concentra-
tion prole diagram (Fig. 4(b)) displays the calculated
concentration at piston height

h(t), plotted in time inter-
vals of length t (fat dots).
Fig. 5 shows the simulated proles of the excess pore
pressure
e
, where this quantity has been normalized with
respect to the maximumvalue o
max
of the applied pressure
o(t). The proles of
e
have been obtained from the
calculated concentration distributions by the equation
c
e
c:
= pq[
co
e
([)
c:
, 0 6: 6h(t)
4548 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
Fig. 4. Simulation of pressure ltration of a concentrated suspension
with [
0
= 0.09 and constant applied pressure o 102.1 kPa: (a)
ltration plot with simulated iso-concentration lines, (b) concentration
proles. The dots denote the concentration at height

h(t) and have
been plotted in time intervals of length t.
Fig. 5. Excess pore pressure proles corresponding to the simulation
of Fig. 4. The dots denote the excess pore pressure at height

h(t)
and have been plotted in time intervals of length t.
Fig. 6. Simulation of pressure ltration of a concentrated suspension
with [
0
= 0.09 and piecewise constant, increasing applied pressure
function o(t): (a) ltration plot with simulated iso-concentration lines,
(b) concentration proles. The dashed lines correspond (from top to
bottom) to piston heights calculated without increasing o(t) at all,
after t =1000 s, and after t =2000 s, respectively. The dots denote the
concentration at height

h(t) and have been plotted in time intervals
of length t.
obtained from inserting (23) into (50). Furthermore, ex-
cess pore pressure values measured obtained in height

h(t) have been plotted into the

e
diagram (Fig. 5).
Fig. 4(a) shows that the ltration rate is decreasing
while the lter cake is accumulating, and remains (nearly)
constant when the cake is formed. At about t = 4100 s,
the ltration stage is terminated, and the piston gets into
contact with the lter cake. The cake is compressed in a
very short time interval and assumes constant porosity.
The excess pore pressure diagram, Fig. 5, illustrates
that the excess pore pressure distribution remains constant
while the ltration rate is constant, as expressed by the
practical coincidence of the simulated
e
proles for t =
2000, 2500, 3000, 3500 and 4000 s. We also observe that
the maximum excess pore pressure value, attained at the
top of the lter cake, diminishes very rapidly once the
piston height has attained the cake surface level.
Figs. 6 and 7 show a simulation with the same ini-
tial concentration, [
0
0.09, but with the non-constant,
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4549
Fig. 7. Excess pore pressure proles corresponding to the simulation
of Fig. 6. The dots denote the excess pore pressure at height

h(t)
and have been plotted in time intervals of length t.
increasing applied pressure function
o(t) =
_

_
102.1 kPa for 0 t 61000 s,
200 kPa for 1000 t 62000 s,
300 kPa for 2000 t 63000 s,
400 kPa for t 3000 s.
(74)
To demonstrate the eect of successively increasing the
applied pressure, the ltration plot, Fig. 6(a), also shows
(as dashed lines, from top to bottom) the piston level
curves that would have been obtained if the applied pres-
sure had been kept constant at 102.1 kPa, if it had been
changed at t = 1000 s only, and if it had been changed
at t = 1000 and 2000 s only, respectively. Comparison
with the calculated curve

h(t) produced by the function
given by (74) clearly shows that each increment of the
applied pressure function produces an increase of the
ltrate ow rate. In addition, we observe that each of
these steps produces an immediate increase of the con-
centration at : =0, which is, of course, a consequence of
Eq. (57). A similar direct consequences of (74) is ex-
pressed by the excess pore pressure prole plot for this
example as shown in Fig. 7.
In the last numerical example (Figs. 8 and 9), we
choose [
0
=[
c
=0.25, i.e. the suspension is initially net-
worked, and let o 200 kPa. With these parameters, no
hindered settling or clear liquid zones form, and the l-
tration process already starts with the expression stage.
Although our parameters were chosen dierently, our nu-
merical results look similar to those displayed in Figs.
14 and 15 of Tiller and Yeh (1987), Fig. 9 of La Heij,
Kerkhof, Herwijn, and Coumans (1996) and Figs. 6 and
7 of SHrensen, Moldrup, and Hansen (1996).
Fig. 8. Simulation of pressure ltration of a concentrated suspension
with [
0
= [
c
= 0.25 and constant applied pressure o 200 kPa: (a)
ltration plot with simulated iso-concentration lines, (b) concentration
proles. The fat curve in (b) displays the concentrations obtained
at height

h(t) and consists of dots that have been plotted in time
intervals of length t.
Fig. 9. Excess pore pressure proles corresponding to the simulation
of Fig. 8. The fat curve shows the excess pore pressure value at height

h(t) and consists of dots that have been plotted in time intervals of
length t.
4550 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
5. Conclusions
We have shown in this paper how the existing theory
of sedimentationconsolidation processes of occulated
suspensions can be extended to a phenomenological the-
ory of cake ltration. The essential ingredients of this
extension are an additional transport term in the known
parabolic degenerate sedimentationconsolidation equa-
tion plus a coupling equation relating to the piston move-
ment, or equivalently the volumetric ow rate, to the
applied pressure. This coupling condition is essentially
based on a balance of the total pressure.
We point out that theories of sedimentation, which are
mostly based on the pioneering paper by Kynch (1952),
rely on the solids concentration as independent variable,
while the traditional ltration theory is based on the solids
pressure, which is represented here by the concept of ef-
fective solid stress. However, since it is usually assumed
that the concentration (or porosity) is a monotonous func-
tion of the solids stress wherever this quantity does not
vanish, it is easily possible to refer ltration theory to
the solids concentration as independent variable. As a
consequence of this reformulation, the permeability now
becomes a porosity-dependent function, and comparison
with the sedimentation equations shows how the perme-
ability equation has to be combined with known constants
in order to form the portion of the ux density function
for the compression zone.
The present formulation and the corresponding nu-
merical method have some decisive advantages as com-
pared to alternate treatments. It is possible to solve the
free-boundary problem, and thereby to simulate the l-
tration process, explicity without having to trace any of
the interfaces (such as the cake}suspension and possi-
bly the suspension}supernate interfaces) that may form
in the interior of the computational domain. This is due
to the fact that the numerical method approximates dis-
continuities in their correct locations, although they may
be smeared out due to numerical diusion, and is in
contrast to more complicated treatments such as that by
SHrensen et al. (1996), in which the cake}suspension in-
terface gives rise to a moving boundary condition. The
conservative form of the numerical algorithm ensures
that the total amount of solids in the lter will remain
constant during simulation (this is true, of course, up to
roundo errors). Moreover, due to the upwind nature of
the algorithm, it can be shown that only such discontinu-
ities are approximated which are physically admissible,
i.e. which satisfy a selection principle or entropy condi-
tion (Evje & Karlsen, 2000; Karlsen & Risebro, 2001).
The related concept of entropy solutions and existence
and uniqueness theory of degenerate parabolic equations
has been fully developed for the cases of gravity settling,
continuous thickening and centrifugation of suspensions
(B urger et al., 2000b; B urger & Karlsen, 2000, 2001),
and is currently being extended to the present free bound-
ary problem. The most obvious benets for the practi-
tioner that accrue from the mathematical analysis include
here the well posedness of the dynamical cake ltration
process and the design of a numerical method for its
computation.
To put the present work in the proper perspective, we
mention that one of the rst investigations of the eect
on sedimentation on ltration of suspensions was under-
taken by Straumann (1963), who, however, considered
incompressible lter cakes only, which correspond here
to the assumption that o
e
0 and hence a 0, A
0. The present model is comparable to those formulated
by Stamatakis and Tien (1991) and Font and Hern andez
(2000). However, Stamatakis and Tien (1991) did not in-
clude the eect of sedimentation and they assumed that
it is the suspensioncake interface where the prescribed
pressure is applied, while in our work, this interface is a
result of the computation. Since the initial cake thickness
was zero, Stamatakis and Tien (1991) had to formulate
a separate procedure to initiate the numerical solution.
Since it is assumed that the suspension concentration is
uniform, as well as that the solids concentration always
equals the critical, it is not clear to us how their approach
can be extended to the expression stage. Similarly, the
cake}suspension interface boundary condition derived by
Atsumi and Akiyama (1975) is based on the assumption
that the void ratio (or, equivalently, the concentration)
of the suspension above the lter cake level is given by
the initial value. To our view, this assumption is overly
restrictive since it precludes that the slurry concentration
directly above the cake may vary, which may occur due
to prescribing non-constant initial data or due to the for-
mation of rarefaction waves, as evident from the charac-
teristics (iso-concentration lines) emerging tangentially
from the cake surface in a ltration plot. Such tangential
characteristics may be produced by considering dierent
shapes of the ux density function [
bk
and selecting pa-
rameters appropriately, and their existence is even explic-
itly postulated in the Font and Hern andez (2000) basic
assumptions.
In future, the numerical method should be employed to
compare the predictions of the phenomenological model
with experimental data, following the papers of B urger
et al. (2000a) and Garrido, B urger, and Concha (2000)
for the case of gravity settling in a closed column. The
published experimental studies that can be utilized for
this purpose include the papers by Sambuichi, Nakakura,
and Osasa (1982), Zydney, Saltzman, and Colton (1989),
Shen, Russel, and Auzerais (1994), Lu, Tung, Pan, and
Hwang (1998), Alles and Anlauf (1998), Wells (1999),
Sis and Chander (2000), Font and Hern andez (2000) and
Fran ca (2000).
Combining the present work with that on centrifuga-
tion of occulated suspensions (B urger & Concha, 2001),
it should be possible to extend the present phenomeno-
logical theory to ltering centrifuges. This will lead to a
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4551
unied theory of solidliquid separation processes, com-
pleting the work by Tiller and Hsyung (1993) and provid-
ing a formal mathematical background to that of Land-
man and White (1994), which can be employed for ac-
curate design calculations.
Notation
a diusion coecient
A integrated diusion coecient
C exponent of the RichardsonZaki ux
density function
[ total solids ux density
[
bk
Kynch batch ux density function
[
1
bk
, [
2
bk
segments of [
bk
[
EO
bk
EngquistOsher ux function
q acceleration of gravity
h piston height

h numerically calculated value of h

h
0
initial piston height
h
m
lter medium thickness
H auxiliary function introduced in Eq. (18)
J number of spatial subintervals
k permeability
k
0
coecient in the permeability equation
(31)

k
0
published value of k
0
for kaolin
k
m
lter medium permeability
K exponent in eective solid stress func-
tion
m
0
total mass of suspension per unit cross-
sectional area of the lter at t = 0
m
d
soliduid dynamic interaction force
MM minmod limiter function
n time step counter
pore pressure

3
third-order interpolation polynomial

A
parameter in the eective solid stress
equation (30)

atm
atmospheric pressure

e
excess pore pressure

t
total pressure
q volume average velocity of the mixture
R
m
lter medium resistance
s
n
)
slopes in numerical method
S
3
piecewise cubic spline function
t time
t
n
time after n iterations
u

coecient of the RichardsonZaki ux

density function
t
f
, t
s
solid and uid phase volocities
t
r
soliduid relative velocity
w scaled height variable
: height
Greek letters
: resistance coecient
[ exponent in the eective solid stress
equation (30)
o exponent in the permeability equation
(31)
t
n
time step
w spatial meshwidth in numerical method
p soliduid density dierence
t plot interval in prole diagrams
c porosity
c positive constant in inequality (63)
[ volumetric solids concentration
[
0
initial concentration
[
c
critical concentration
[
n
)
approximate value of [
[
n, L
)
, [
n, R
)
extrapolated variables
[
max
maximum concentration
j
f
uid dynamic viscosity
p local density of the mixture
p
s
, p
f
solid and uid mass densities
o applied pressure
o
0
coecient in eective solid stress func-
tion
o
e
eective solid stress
o
1
e
inverse function of o
e
([) for [[
c
Acknowledgements
Preparation of this work was made possible through
the Applied Mathematics in Industrial Flow Problems
(AMIF) programme of the European Science Founda-
tion (ESF); through the Sonderforschungsbereich 404 at
the University of Stuttgart; and through Fondef Project
D97-I2042 at the University of Concepci on.
References
Atsumi, K., & Akiyama, T. (1975). A study of cake ltration:
Formulation as a Stefan problem. Journal of Chemical
Engineering of Japan, 8, 487492.
Auzerais, F. M., Jackson, R., & Russel, W. B. (1988). The resolution
of shocks and the eects of compressible sediments in transient
settling. Journal of Fluid Mechanics, 195, 437462.
B urger, R. (2000). Phenomenological foundation and mathematical
theory of sedimentationconsolidation processes. Chemical
Engineering Journal, 80, 177188.
B urger, R., & Concha, F. (1998). Mathematical model and numerical
simulation of the settling of occulated suspensions. International
Journal of Multiphase Flow, 24, 10051023.
B urger, R., & Concha, F. (2001). Settling velocities of particulate
systems: 12. Batch centrifugation of occulated supensions.
International Journal of Mineral Processing, in press.
B urger, R., Concha, F., & Tiller, F. M. (2000a). Applications of the
phenomenological theory to several published experimental cases
of sedimentation processes. Chemical Engineering Journal, 80,
105117.
4552 R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553
B urger, R., Evje, S., & Karlsen, K. H. (2000b). On strongly degene-
rate convectiondiusion problems modeling sedimentation
consolidation processes. Journal of Mathematical Analysis and
Applications, 247, 517556.
B urger, R., & Karlsen, K. H. (2000). A strongly degenerate
convectiondiusion problem modeling centrifugation of
occulated suspensions. Proceedings of the Eighth International
Conference on Hyperbolic Problems (Hyp 2000), Magdeburg,
Germany, February 28March 3, 2000, in press.
B urger, R., & Karlsen, K. H. (2001). On some upwind nite dierence
schemes for the phenomenological sedimentationconsolidation
model. Journal of Engineering Mathematics, in press.
B urger, R., & Tory, E. M. (2000). On upper rarefaction waves in
batch settling. Powder Technology, 108, 7487.
B urger, R., Wendland, W. L., & Concha, F. (2000c). Model
equations for gravitational sedimentationconsolidation processes.
Zeitschrift f ur Angewandte Mathematik und Mechanik, 80, 79
92.
Buscall, R., & White, L. R. (1987). On the consolidation of
concentrated suspensions, Part 1: The theory of sedimentation.
Journal of the Chemical Society, Faraday Transactions I, 83,
873891.
Bustos, M. C., Concha, F., B urger, R., & Tory, E. M. (1999).
Sedimentation and thickening: Phenomenological foundation and
mathematical Theory. Dordrecht, Netherlands: Kluwer Academic
publishers.
Concha, F. (2001). Manual de Filtraci on & Separaci on, Centro de
Tecnologa Mineral, Universidad de Concepci on}Fundaci on Chile,
Concepci on, Chile.
Concha, F., & Bascur, O. (1977). Phenomenological model of
sedimentation. Proceedings of the 12th International Mineral
Processing Congress (XII IMPC), Vol. 4 (pp. 2946). S ao Paulo,
Brazil.
Concha, F., & Bustos, M. C. (1986). Theory of sedimentation of
occulated ne particles. In B. M. Moudgil, & P. Somasundaran
(Eds.), Proceedings of the engineering foundation conference on
occulation, sedimentation and consolidation (pp. 275284). Sea
Island, GA; New York: AIChE.
Concha, F., Bustos, M. C., & Barrientos, A. (1996).
Phenomenological theory of sedimentation. In E. M. Tory (Ed.),
Sedimentation of small particles in a viscous uid (pp. 5196).
Southampton, UK: Computational Mechanics Publications.
Darcy, H. (1856). Les Fontaines Publiques de la Ville de Dijon.
Paris: Dalmont.
Diplas, P., & Papanicolaou, A. N. (1997). Batch analysis of slurries
in zone settling regime. Journal of Environmental Engineering,
123, 659667.
Engquist, B., & Osher, S. (1981). One-sided dierence
approximations for non-linear conservation laws. Mathematics of
Computation, 36, 321351.
Espedal, M. S., & Karlsen, K. H. (2000). Numerical solution of
reservoir ow models based on large time step operator splitting
algorithms. In A. Fasano, & H. van Duijn (Eds.), Filtration
in porous media and industrial applications. Lecture Notes in
Mathematics, Vol. 1734 (pp. 977). Berlin: Springer.
Evje, S., & Karlsen, K. H. (2000). Monotone dierence
approximations of BJ solutions to degenerate convectiondiusion
equations. SIAM Journal of Numerical Analysis, 37, 18381860.
Font, R., Garca, P., & P erez, M. (1998). Analysis of the variation
of the upper discontinuity in sedimentation batch test. Separation
Science Technology, 33, 14871510.
Font, R., & Hern andez, A. (2000). Filtration with sedimentation:
application of Kynchs theorems. Separation Science Technology,
35, 183210.
Fran ca, S. C. A. (2000). Equac oes constitutivas para a sedimentac ao
de suspens oes oculentas. Doctoral Thesis, Federal University of
Rio de Janeiro, Brazil.
Garrido, P., B urger, R., & Concha, F. (2000). Settling velocities
of particulate systems: 11. Comparison of the phenomenological
sedimentationconsolidation model with experimental results.
International Journal of Mineral Processing, 60,
213227.
Garrido, P., Concha, F., & B urger, R. (2001). Application of
the unied model of solidliquid separation of occulated
suspensions to experimental results. Proceedings of the VI
Southern Hemisphere Meeting on Minerals Technology, Rio de
Janeiro, Brazil, May 2731, 2001, in press.
Karlsen, K. H., & Risebro, N. H. (2001). Convergence of nite
dierence schemes for viscous and inviscid conservation laws with
rough coecients. Preprint, University of Bergen by Mathematical
Modelling and Numerical Analysis, 35, 239270.
Kynch, G. S. (1952). A theory of sedimentation. Transactions of the
Faraday Society, 48, 239270.
La Heij, E. J., Kerkhof, P. J. A. M., Herwijn, A. J. M., & Coumans, W.
J. (1996). Fundamental aspects of sludge ltration and expression.
Water Research, 30, 697703.
Landman, K. A., & Russel, W. B. (1993). Filtration at large pressures
for strongly occulated suspensions. Physics of Fluids A, 5, 550
560.
Landman, K. A., Sirako, C., & White, L. R. (1991). Dewaterung
of occulated suspensions by pressure ltration. Physics of Fluids
A, 3, 14951509.
Landman, K. A., & White, L. R. (1994). Solid}liquid separation
of occulated suspensions. Advances in Colloid and Interface
Science, 51, 175246.
Lu, W. M., Tung, K. L., Pan, C. H., & Hwang, K. J. (1998). The
eect of particle sedimentation on gravity ltration. Separation
Science Technology, 33, 17231746.
Massarani, G. (1997). Fluidodin amica em sistemas particulados.
Editora UFRJ, Federal University of Rio de Janeiro, Brazil.
Michaels, A. S., & Bolger, J. C. (1962). Settling rates and
sediment volumes of occulated Kaolin suspensions. Industrial
and Engineering Chemistry Fundamentals, 1, 2433.
Richardson, J. F., & Zaki, W. N. (1954). Sedimentation and
uidization I. Transactions of the Institution of Chemical
Engineers (London), 32, 3553.
Rushton, A., Ward, A. S., & Holdich, R. G. (2000). Solid
liquid ltration and separation technology (2nd ed.). Weinheim,
Germany: Wiley-VCH.
Sambuichi, M., Nakakura, H., & Osasa, K. (1982). The eect of
gravity settling on constant pressure ltration. Memoirs of the
Faculty of Engineering, Yamaguchi University, 33, 6570.
Scott, K. J. (1968). Experimental study of continuous thickening of
a occulated silica slurry. Industrial and Engineering Chemistry
Fundamentals, 7, 582595.
Shen, C., Russel, W. B., & Auzerais, F. (1994). Colloidal gel ltration:
experiment and theory. [Link].E. Journal, 40, 18761891.
Shirato, M., Sambuichi, M., Kato, H., & Aragaki, H. (1969).
Internal ow mechanism in lter cakes. [Link].E. Journal, 15,
405409.
Sis, H., & Chander, S. (2000). Pressure ltration of dispersed
and occulated alumina slurries. Minerals & Metallurgical
Processing, 17, 4148.
SHrensen, P. B., Moldrup, P., & Hansen, J. A. (1996). Filtration and
expression of compressible cakes. Chemical Engineering Science,
51, 967979.
Stamatakis, K., & Tien, C. (1991). Cake formation and growth in
cake ltration. Chemical Engineering Science, 46, 19171933.
Straumann, R. (1963). Einuss de Sedimentation auf die Filtration,
ChemieIngenieurTechnik, 35, 715720.
Tiller, F. M., Crump, J. R., & Ville, F. (1980). A revised approach
to the theory of cake ltration. In P. Somasundaran (Ed.),
Proceedings of the international symposium on ne particles
processing, Vol. 2 (pp. 15491582). New York: American Institute
of Mining Metallurgical and Petroleum Engineers.
R. B urger et al. / Chemical Engineering Science 56 (2001) 45374553 4553
Tiller, F. M., & Horng, L. (1983). Hydraulic deliquoring of
compressible lter cakes. [Link].E. Journal, 29, 297305.
Tiller, F. M., & Hsyung, N. B. (1993). Unifying the theory of
thickening, ltration and centrifugation. Water Science Tech., 28,
19.
Tiller, F. M., Hsyung, N. B., & Cong, D. Z. (1995). Role of porosity
in ltration: XII. Filtration with sedimentation. [Link].E. Journal,
41, 11531164.
Tiller, F. M., & Khatib, Z. (1984). The theory of sediment volumes
of compressible, particulate structures. Journal of Colloid and
Interface Science, 100, 5567.
Tiller, F. M., & Kwon, J. H. (1998). Role of porosity in ltration:
XIII. Behavior of highly compactible cakes. [Link].E. Journal,
44, 21592167.
Tiller, F. M., & Leu, W. F. (1980). Basic data tting in ltration.
Journal of the Chinese Institute of Chemical Engineers, 11,
6170.
Tiller, F. M., Lu, R., Kwon, J. H., & Lee, D. J. (1999). Variable
ow rate in compactible lter cakes. Water Research, 33,
1522.
Tiller, F. M., & Yeh, C. S. (1987). The role of porosity in ltration;
Part XI: Filtration followed by expression. [Link].E. Journal, 33,
12411256.
Wakeman, R. J., & Tarleton, E. S. (1999). Filtration: Equipment
selection, modelling and simulation. Oxford, UK: Elsevier
Science.
Wells, S. (1999). Analytical modeling of cake ltration. In W.
Leung (Ed.), Advances in ltration and separation technology,
Vol. 12 (pp. 158165). American Filtration Separation Society,
Tuscaloosa, AL.
Zydney, A. L., Saltzman, W. M., & Colton, C. K. (1989). Hydraulic
resistance of red cell beds in an unstirred ltration cell. Chemical
Engineering Science, 44, 147159.