a generalized kinetic model of sedimentation of ... - World Scientific

Comment

Report 7 Downloads 45 Views

October 14, 2008 11:9 WSPC/103-M3AS

00318

Mathematical Models and Methods in Applied Sciences Vol. 18, No. 10 (2008) 1741–1785 c World Scientiﬁc Publishing Company

A GENERALIZED KINETIC MODEL OF SEDIMENTATION OF POLYDISPERSE SUSPENSIONS WITH A CONTINUOUS PARTICLE SIZE DISTRIBUTION

¨ RAIMUND BURGER Departamento de Ingenier´ıa Matem´ atica, Facultad de Ciencias F´ısicas y Matem´ aticas, Universidad de Concepci´ on, Casilla 160-C, Concepci´ on, Chile [email protected] ANTONIO GARCIA Departamento de Ingenier´ıa Metal´ urgica, Facultad de Ingenier´ıa y Ciencias Geol´ ogicas, Universidad Cat´ olica del Norte, Avenida Angamos 0610, Antofagasta, Chile [email protected] MATTHIAS KUNIK Institut f¨ ur Analysis und Numerik, Fakult¨ at f¨ ur Mathematik, Universit¨ at Magdeburg, Universit¨ atsplatz 2, D-39106 Magdeburg, Germany [email protected] Received 31 October 2007 Revised 21 December 2007 Communicated by N. Bellomo

Polydisperse suspensions with particles of a ﬁnite number N of size classes have been widely studied in laboratory experiments. However, in most real-world applications the particle sizes are distributed continuously. In this paper, a well-studied one-dimensional kinematic model for batch sedimentation of polydisperse suspensions of small equaldensity spheres is extended to suspensions with a continuous particle size distribution. For this purpose, the phase density function Φ = Φ(t, x, ξ), where ξ ∈ [0, 1] is the normalized squared size of the particles, is introduced, whose integral with respect to ξ on an interval [ξ1 , ξ2 ] is equivalent to the volume fraction at (t, x) occupied by particles of that size range. Combining the Masliyah–Lockett–Bassoon (MLB) model for the solidﬂuid relative velocity for each solids species with the concept of phase density function yields a scalar, ﬁrst-order equation for Φ, namely the equation of the generalized kinetic theory. Three numerical schemes for the solution of this equation are introduced, and a numerical example and an L1 error study show that one of these schemes introduces less numerical diﬀusion and less spurious oscillations near discontinuities than the others. Several numerical examples illustrate the simulated behavior of this kind of suspensions. Numerical results also illustrate the solution of an eigenvalue problem associated with the equation of the generalized kinetic theory. 1741

October 14, 2008 11:9 WSPC/103-M3AS

1742

00318

R. B¨ urger, A. Garc´ıa & M. Kunik Keywords: Sedimentation; polydisperse suspension; continuous particle size distribution; conservation law; kinetic model; ﬁnite diﬀerence method; numerical simulation. AMS Subject Classiﬁcation: 76T20, 35L65, 65M06

1. Introduction 1.1. Scope Numerous engineering applications involve the sedimentation of small solid particles dispersed in a viscous ﬂuid. In polydisperse suspensions, the particles diﬀer in size or density and segregate, creating areas of diﬀerent composition, which is frequently described by spatially one-dimensional models. We assume that the particle diameters are small compared with that of the settling vessel, which justiﬁes identifying each species with a continuous phase. In general, we distinguish between N species of particles, where particles of species i, associated with the volume fraction φi , have size di and density ρi , and di = dj or ρi = ρj for i = j. If vi is the phase velocity of species i, then the continuity equations of the N species in diﬀerential form are ∂t φi + ∂x (φi vi ) = 0 for i = 1, . . . , N , where t is time and x is depth. In kinematic models, the velocities v1 , . . . , vN are assumed to be given functions of the vector φ := φ(t, x) := (φ1 (t, x), . . . , φN (t, x))T of local concentrations. This yields systems of conservation laws of the type ∂t φ + ∂x f (φ) = 0, fi (φ) = φi vi (φ),

f (φ) = (f1 (φ), . . . , fN (φ))T , i = 1, . . . , N.

(1.1)

The algebraic expressions of the functions f1 (φ), . . . , fN (φ) are introduced in Sec. 2.1. If φmax is a maximum solids concentration, then the solution φ is assumed to take values in D := {φ = (φ1 , . . . , φN )T ∈ RN | φ1 ≥ 0, . . . , φN ≥ 0, φ1 + · · · + φN ≤ φmax }. In this paper, we study a model of sedimentation of polydisperse suspensions that emerges from the kinematic model (1.1) if we consider equal-density particles, but let N tend to inﬁnity, i.e. if the suspension has a continuous particle size distribution. This property describes many real suspensions more realistically than a ﬁnite number of species. The resulting mathematical model is an equation for the phase density function Φ(t, x, ξ), where the normalized squared particle size ξ is the kinetic parameter, such that integrating Φ(t, x, ξ) over a ξ-interval yields the total solids concentration within that size range. This equation, to which we refer as the equation of the generalized kinetic theory, replaces (1.1). (Note that we distinguish between a kinematic model and a generalized kinetic equation.) We ﬁrst analyze an eigenvalue problem associated with the equation of the generalized kinetic theory. Then we propose three alternative ﬁnite volume schemes for the discretization of this equation, two of which are adapted to the “concentration times velocity” structure of the ﬂux functions of the kinematic model, and were advanced recently.19 A series of numerical examples made with the best-performing

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1743

scheme illustrates the settling behavior predicted by the kinetic model. Numerical results suggest that solutions to a succession of problems with increasing N converge to a solution of the equation of the generalized kinetic theory. Finally, two numerical experiments illustrate a result of the eigenvalue analysis, which states that the spectrum of the eigenvalue problem associated with the equation of the generalized kinetic theory is contained in a particular bounded interval. 1.2. Related work The model of sedimentation of polydisperse suspensions used herein goes back to Masliyah37 and Lockett and Bassoon36 (“MLB model”). General overviews of polydisperse sedimentation models that give rise to (1.1), but with ﬂuxes f1 , . . . , fN that possibly diﬀer from those of the MLB model, are given in Refs. 32, 35 and 52 and the introductory parts of Refs. 16, 17, 21 and 22. Solutions of (1.1) are in general discontinuous, and include kinematic shocks that separate areas of diﬀerent composition. The strongly coupled, nonlinear nature of (1.1) makes it diﬃcult to determine exact solutions. In fact, piecewise constant solutions with areas of constancy separated by kinematic shocks are constructed in Refs. 31 and 47; however, these solutions may fail to satisfy an admissibility condition.14 Numerical techniques that have been applied to (1.1) include universal ﬁnite volume schemes for conservation laws,16–18 schemes designed for kinematic ﬂow models,19 multiresolution methods that concentrate computational eﬀort near kinematic shocks,22 and the front tracking method.14 The papers by Xue and Sun50 and Zeidan et al.51 report the implementation of the Kurganov–Tadmor34 and Godunov ﬁnite volume methods for the sedimentation model. Alternatively, the MLB model explicitly provides the settling velocities of individual particles in a particle-based approach.3–8 Similar kinematic models describe multi-species traﬃc ﬂow12, 19, 22, 49, 53–55 and the settling of oil-in-water dispersions.46 We recall that for a given vector φ, (1.1) is called hyperbolic if the Jacobian Jf (φ) := (∂fi (φ)/∂φj )i,j=1,...,N has real eigenvalues only, and strictly hyperbolic if these are moreover pairwise distinct. Hyperbolicity ensures that concentration changes, and solution information in general, propagate at ﬁnite speed. Recent advances were made in the hyperbolicity analysis and characterization of eigenvalues of multi-species kinematic models. For the model of settling of oil-in-water dispersions, Rosso and Sona46 proved for arbitrary N strict hyperbolicity in the interior of D. The proof is based on deriving an explicit closed formula of the characteristic polynomial p(λ; φ) = det(Jf (φ) − λI) by exploiting elimination possibilities in the determinant. Then p(λ; φ) is evaluated at N suitable λ-arguments that produce values of alternating sign, which along with a discussion of p(λ; φ) for λ → ±∞ implies that p(λ; φ) must have N distinct zeros. Berres et al.16 proved in a similar way that the MLB model16, 36, 37 with equal-density particles is strictly hyperbolic for an arbitrary number N of species (size classes). The basic idea was also used by Zhang et al.54 to prove strict hyperbolicity of the multi-class traﬃc model proposed in Refs. 12 and 49.

October 14, 2008 11:9 WSPC/103-M3AS

1744

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

Relatively few experimental and theoretical studies have explicitly been concerned with the settling of a suspension with a continuous particle size distribution. Greenspan and Ungarish31 proposed a model which is roughly equivalent to our equation of the generalized kinetic theory (set out in Sec. 2.5) for the cumulative phase density. They essentially suggested an iterative procedure to construct approximate solutions of their model either by successions of kinematic shock constructions for ﬁnite N or, for the particular case of Riemann-type problems, by iteratively solving an integral equation. The formation of sediment with a continuous particle size distribution is, however, not considered. Austin et al.2 studied suspensions with particles having two diﬀerent densities and a continuous size distribution due to Rosin and Rammler,45 which is widespread in mineral processing (Sect. 20 of Ref. 39). The continuous distribution is approximated by a discrete one, and their model features a regularizing diﬀusion term. Experimental results, though no complete mathematical models for their interpretation, for suspensions with continuous particle size distribution are given in Refs. 24 and 25. Xue and Sun50 present numerical solutions of the MLB model with N = 35, and obtain good agreement of the numerically predicted total solids concentration with light extinction measurements made at several heights (but well above the level of the forming sediment). Substantial eﬀorts have been made to derive collective (macroscopic) descriptions from individual (microscopic or mesoscopic) models; standard examples are the convergence of the Boltzmann equation to the Euler equations of gas dynamics or the derivation of macroscopic traﬃc models; see Refs. 1, 11, 13, 26–28 for overviews. Moreover, considerable eﬀorts have been made in mathematical biology in deriving macroscopic equations from the underlying microscopic description delivered by the kinetic theory, see for example Filbet et al.30 or Bellomo et al.10 The latter approaches are of particular interest for the application considered herein, since they can be generalized to large systems of dispersed elements. In a sense, we address the opposite limit, since in our kinetic system the convective term is nonlinear and nonlocal (with respect to the kinetic parameter ξ), but the collision terms are zero. Our model is similar to (but more involved than) a model of multi-species traﬃc ﬂow, which was proposed for a ﬁnite number of species in Refs. 12 and 49, further studied in Refs. 19, 22, 53–55 and analyzed in the kinetic limit of an inﬁnite number of species in Ref. 13. We brieﬂy relate our ﬁndings to these works in Sec. 6. 1.3. Outline of this paper The mathematical model is introduced in Sec. 2. To this end, we consider in Sec. 2.1 a suspension of particles having the same density s , and belonging to a ﬁnite number N of size classes or species. We restate the MLB model describing the settling of this mixture, which yields a system of conservation laws of type (1.1) with a particular nonlinear ﬂux vector. In Sec. 2.2, we recall a theorem from Ref. 16, which roughly states that for concentration vectors φ from the interior of the phase

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1745

space, the system is strictly hyperbolic. In addition, inclusions for the eigenvalues are stated. The model for ﬁnite N is extended in Sec. 2.3 to a kinetic model for a continuous particle size distribution, where the normalized squared size ξ is the kinetic parameter and, as stated above, the phase density function Φ is sought. The ingredients of the kinetic model are summarized in Sec. 2.4. In Sec. 2.5, we state a separate equation of the generalized kinetic theory for the cumulative phase density Φ∗ , which is the primitive of Φ with respect to ξ, and has a direct physical interpretation (namely, Φ∗ is a volume fraction) and better regularity than Φ. Section 3 is devoted to the eigenvalue problem associated with the equation of the generalized kinetic theory. In Sec. 3.1 we pass from the eigenvalue problems of the hyperbolic systems with ﬁnite N , which are discussed in Sec. 2.2, to the eigenvalue problem for Φ∗ ; both give rise to the same eigenvalues. In Sec. 3.2, we prove that the Jacobians of the cumulative problem, evaluated at a ﬁxed function Φ∗ , converge in L1 to a continuous operator when N → ∞. (The proof of this statement, Lemma 3.2, is deferred to the Appendix.) We then characterize the spectrum of the continuous operator, which is contained in a bounded interval. In Sec. 3.3 we determine the eigenfunctions associated with the eigenvalue problem, under the assumption that both Φ∗ and the eigenfunctions are smooth. (This assumption will, however, not hold in general.) In Sec. 4, we deﬁne three numerical schemes to approximate solutions of the kinetic model. The discretization parameters and boundary conditions are introduced in Sec. 4.1, while three alternative discretizations (Schemes 1, 2 and 3) are introduced in Sec. 4.2. All schemes are based on an exact integration with respect to ξ, so that conservativity is ensured. While Scheme 1 is a version of the well-known Lax–Friedrichs scheme for conservation laws, Schemes 2 and 3 are versions of a family of schemes introduced in Ref. 19 that explicitly utilize the algebraic structure of the numerical ﬂux for kinematic models. We propose in Sec. 4.3 a numerical test for mass conservation. Numerical examples are presented in Sec. 5. In Sec. 5.1 (Examples 1, 2 and 3), we deal with problems of batch settling. We ﬁrst compare the performance of the three schemes, with the conclusion that Scheme 1 suﬀers from excessive numerical viscosity, while Scheme 2 produces slightly oscillatory solutions. Scheme 3 is chosen for the remainder of numerical experiments. In Sec. 5.2, we study the Riemann problem for the equation of the generalized kinetic theory (Examples 4 and 5), and display numerical results combined with an eigenvalue analysis of the numerical solution to illustrate the ﬁndings of Sec. 3. Some conclusions of this paper are collected in Sec. 6. 2. The Mathematical Model 2.1. Polydisperse suspensions with a discrete particle size distribution We assume that particles of species i have size (diameter) di , where d1 > d2 > . . . > dN > 0. According to the MLB model,21, 36, 37 the sedimentation of this mixture is

October 14, 2008 11:9 WSPC/103-M3AS

1746

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

modeled by the ﬁrst-order system of N conservation laws ∂t φi + ∂x (φi [q + ui − (φ1 u1 + . . . + φN uN )]) = 0,

i = 1, . . . , N,

(2.1)

where t is time, x is depth, φi is the volume fraction of particle species i, ui is the solid-ﬂuid relative velocity of particle species i, and q is the volume average velocity of the mixture, which can be controlled externally. We limit ourselves to a closed settling column, for which q = 0. The speciﬁc equation of the MLB model is ui =

¯s g 2 ˜ d V (φ)(1 − φ), 18µf i

(2.2)

where ¯s := s −f , g is the acceleration of gravity, f is the density of the ﬂuid, µf is its viscosity, and φ = φ1 + . . . + φN is the total solids volume fraction. The function V (φ) := V˜ (φ)(1 − φ) is the hindered settling factor. We choose the continuous function  γ  for φ ∈ [0, φq ),  (1 − φ) γ ˜ V (φ) = (1 − φq ) (φmax − φ)/(φmax − φq ) for φ ∈ [φq , φmax ], (2.3)   0 otherwise, with γ ≥ 1 and the maximum solids concentration φmax ∈ (0, 1]. We assume that φ˜q := (γφmax − 1)/(γ − 1) satisﬁes 0 < φ˜q < φmax ; in this case, setting φq := φ˜q ensures that (2.3), and therefore V (φ), are diﬀerentiable across φ = φq . (Note that (2.3) is a slight modiﬁcation of the common Richardson–Zaki formula.44 ) In the sequel, we will always assume that V ∈ C 1 [0, φmax ]. For φ = 0, (2.2) recovers the well-known formula ui = ¯s gd2i /(18µf ) for the Stokes velocity, i.e. the settling velocity in an unbounded ﬂuid, of a particle of size di . For detailed derivations of (2.1)–(2.3), statements of the underlying assumptions and discussions of the limitations of the model, we refer to Refs. 15, 16, 21, 36 and 37. To further simplify the model, let v1 := ¯s gd21 /(18µf ) be the settling velocity of the largest species. We consider a settling column of height L and introduce the dimensionless space variable x ˜ = x/L, the time variable t˜ = (v1 /L)t, and 2 2 the parameters δi := di /d1 , i = 1, . . . , N . Inserting (2.2) into (2.1) for q = 0 and switching back to (t, x) instead of (t˜, x ˜), we obtain the governing system of conservation laws ∂t φ + ∂x f (φ) = 0,

(2.4)

where the components f1 (φ), . . . , fN (φ) of the ﬂux vector f (φ) are deﬁned by fi (φ) = φi (δi − (δ1 φ1 + . . . + δN φN ))V (φ),

i = 1, . . . , N.

(2.5)

The zero-ﬂux boundary conditions corresponding to a closed column can then be stated as fi |x=0 = fi |x=1 = 0,

i = 1, . . . , N.

(2.6)

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1747

Finally, we prescribe an initial concentration distribution φ0 (x) = (φ01 (x), . . . , φ0N (x))T , i.e. φi (t, x)|t=0 = φ0i (x),

x ∈ [0, 1],

i = 1, . . . , N ;

φ0 (x) ∈ D for x ∈ [0, 1]. (2.7)

2.2. Properties of the model for a finite number of particle sizes The system (2.5) is strictly hyperbolic for all N ∈ N, size distributions δ1 = 1 > δ2 > . . . > δN , and vectors φ with φ1 , . . . , φN > 0 and φ < φmax . More precisely, the following theorem was proved in Ref. 16. Theorem 2.1. (Ref. 16) For any vector φ = (φ1 , . . . , φN )T ∈ D we recall the abbreviation φ = φ1 + · · · + φN and define the quantities δ¯ := δ1 φ1 + . . . + δN φN ,

¯ (φ) + V (φ)(δ¯ + φ). λ− := −2δV

(2.8) ˜

Assume that the hindered settling function V˜ (φ) satisfies V˜ (φ) > 0 and V (φ) < 0 for φ ∈ [0, φmax ). If δ1 = 1 > δ2 > . . . > δN and φ is chosen from the interior of D, then the eigenvalues λN (φ) ≤ λN −1 (φ) ≤ . . . ≤ λ1 (φ) of the Jacobian Jf (φ) satisfy ¯ V (φ)(δi − δ)), ¯ λi (φ) ∈ (V (φ)(δi+1 − δ),

i = 1, . . . , N − 1,

¯ λN (φ) ∈ (λ− , V (φ)(δN − δ)).

(2.9) (2.10)

Remark 2.1. Note that V (φ) > 0 is guaranteed in the previous theorem due to (2.3) since we have assumed that φ is chosen in the interior of D, i.e. φ < φmax . Thus the ﬁrst N − 1 eigenvalues λi (φ) in (2.9) indicate a continuous distribution in the limit N → ∞, whereas the last (smallest) eigenvalue λN (φ) plays a separate role in the theory. Corollary 2.1. Under the conditions of Theorem 2.1, for any vector φ ∈ D we have that ¯ V (φ)(δi − δ)], ¯ λi (φ) ∈ [V (φ)(δi+1 − δ), λN (φ) ∈ [λ− , V (φ)(δN

i = 1, . . . , N − 1, ¯ − δ)].

(2.11) (2.12)

In particular, if V (φ) = 0, then at most two of the eigenvalues coincide. More precisely, if λi (φ) = λi+1 (φ) for one index i ∈ {1, . . . , N − 1}, then λi (φ) = ¯ and λj (φ) < λi (φ) for j = i + 1, . . . , N and λj (φ) > λi+1 (φ) = V (φ)(δi+1 − δ), λi (φ) for j = 1, . . . , i − 1. If V (φ) = 0, then λi (φ) = 0 for i = 1, . . . , N − 1 and λN (φ) ∈ [V (φ)(δ¯ + φ), 0]. Proof. Corollary 2.1 follows from Theorem 2.1 if we take into account that for 0 ≤ φ ≤ φmax , Jf is a diﬀerentiable function of φ, and use the well-known fact that the eigenvalues of an N × N matrix are H¨ older continuous functions of its 38 coeﬃcients with index 1/N .

October 14, 2008 11:9 WSPC/103-M3AS

1748

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

2.3. Polydisperse suspensions with a continuous particle size distribution Polydisperse suspensions with particles of N distinct size classes have been utilized in many laboratory experiments. However, in most applications, for example in mineral processing, the sizes of particles are continuously distributed. Since most real suspensions have passed through a sieve with a determined mesh-width, we can assume that there exists a maximum particle size, which is normalized to one as in the discrete case. Consequently, the local composition of the polydisperse suspension as a function of time t and depth x is no longer determined by a vector φ ξ2= φ(t, x), but by a phase density function Φ = Φ(t, x, ξ), ξ ∈ [0, 1], where ξ1 Φ(t, x, ξ) dξ denotes the volume fraction at (t, x) occupied by all particles with normalized squared size ξ1 < ξ < ξ2 . To derive an equation for the evolution of Φ, we may approximate the continuous size distribution by a discrete one, multiply the ith equation in (2.4) by δi − δi−1 , sum the equations from i1 to i2 , and then require that the summed equation holds for all choices 1 ≤ i1 ≤ i2 ≤ N . Passing to N → ∞ and replacing the sum by an integration, we obtain that the following equation holds for all 0 ≤ ξ1 < ξ2 ≤ 1: ξ2 ¯ x)))} dξ = 0, ¯ x))(ξ − ξ(t, {∂t Φ(t, x, ξ) + ∂x (Φ(t, x, ξ)V (Φ(t, (2.13) ξ1

¯ and the local average value ξ¯ of ξ with respect where the total volume fraction Φ to Φ are given by the respective expressions 1 1 ¯ x) := ¯ x) := Φ(t, Φ(t, x, ξ) dξ, ξ(t, ξ Φ(t, x, ξ) dξ. (2.14) 0

0

Since (2.13) is assumed to hold for all suitable ξ1 < ξ2 , the integrand in (2.13) should vanish, which leads to the following equation of the generalized kinetic theory, ¯ and ξ: ¯ written down compactly by dropping the arguments t, x and ξ of Φ, Φ ¯ = 0. ¯ − ξ)) ∂t Φ + ∂x (ΦV (Φ)(ξ

(2.15)

2.4. Kinetic problem in final form For easy reference, we collect here the ingredients of the kinetic model in ﬁnal and dimensionless form, where we have put L = 1. For a given ﬁnal time T > 0, we seek a function Φ : [0, T ] × [0, 1] × [0, 1] → [0, ∞) which satisﬁes the equation of the generalized kinetic theory ¯ x))) = 0, ¯ x))(ξ − ξ(t, ∂t Φ(t, x, ξ) + ∂x (Φ(t, x, ξ)V (Φ(t,

t ∈ (0, T ],

x ∈ (0, 1) (2.16)

subject to the initial condition Φ(0, x, ξ) = Φ0 (x, ξ),

x ∈ [0, 1],

Φ0 : [0, 1] × [0, 1] → [0, ∞),

(2.17)

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1749

and zero-ﬂux boundary conditions, i.e. F (t, 0, ξ) = F (t, 1, ξ) = 0

for all ξ ∈ [0, 1] and all t ∈ [0, T ]

(2.18)

with the “kinetic ﬂux function” ¯ x)). ¯ x))(ξ − ξ(t, F (t, x, ξ) := Φ(t, x, ξ)V (Φ(t,

(2.19)

2.5. Cumulative quantities for the conservation of mass In this subsection, we introduce a concept that will be important for the eigenvalue analysis of the hyperbolic system (2.16) in the asymptotic case N → ∞ of inﬁnitely many species. We ﬁrst deﬁne the primitive Φ∗ of the phase density Φ with respect to the kinetic parameter ξ by ξ ∗ Φ(t, x, ξ ) dξ (2.20) Φ (t, x, ξ) := 0

and observe that Φ∗ (t, x, ξ2 ) − Φ∗ (t, x, ξ1 ) =

ξ2

Φ(t, x, ξ ) dξ

(2.21)

ξ1

determines the volumetric concentration at (t, x) occupied by all particles with dimensionless size 0 ≤ ξ1 ≤ ξ2 ≤ 1. Then we can also rewrite the moments (2.14) of Φ in terms of Φ∗ as 1 ∗ ∗ ¯ ¯ Φ∗ (t, x, ξ) dξ, (2.22) Φ(t, x) = Φ (t, x, 1), ξ(t, x) = Φ (t, x, 1) − 0

¯ x). Applying an integration by where we have applied integration by parts for ξ(t, parts to (2.16) yields the conservation law of mass for all species for all t ∈ (0, T ] and x ∈ (0, L):

ξ ∗ ∗ ∗ ∗ ¯ Φ (t, x, ξ) dξ = 0. ∂t Φ (t, x, ξ) + ∂x V (Φ (t, x, 1)) Φ (t, x, ξ)(ξ − ξ(t, x)) − 0

(2.23) This may be regarded as a separate equation of the generalized kinetic theory for the so-called cumulative phase density Φ∗ , but we will not use (2.22) and (2.23) for numerical purposes. Equation (2.23) contains exactly the same information as the original equation (2.16), but has two advantages: The quantity Φ∗ is a standard volume fraction, and thus has a direct physical meaning. The second advantage is important for the eigenvalue analysis in the asymptotic case N → ∞, namely the better regularity properties of the cumulative quantity Φ∗ in the case that Φ is a distributional or singular solution, which cannot be excluded a priori. We remark that the kinetic model (2.16)–(2.19) includes the model (2.4)–(2.7) for ﬁnite N as a special case. For example, if the initial condition (2.7) is speciﬁed

October 14, 2008 11:9 WSPC/103-M3AS

1750

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

for the normalized squared particle sizes δ1 = 1 > δ2 > · · · > δN , then the model for ﬁnite N is recovered if the initial distribution Φ0 in (2.17) is deﬁned by Φ0 (x, ξ) =

N

φ0i (x)δ(ξ − δi ),

ξ, x ∈ [0, 1],

(2.24)

i=1

where δ(·) is the Dirac δ-distribution. For the cumulative version (2.23) of (2.16), (2.24) is replaced by the following initial condition, where H(·) is the Heaviside function: N φ0i (x)H(ξ − δi ), ξ, x ∈ [0, 1]. (2.25) Φ∗ (0, x, ξ) = Φ∗0 (x, ξ) = i=1

3. The Eigenvalue Problem for the Equation of the Generalized Kinetic Theory 3.1. The eigenvalue problem for finite N For the subsequent analysis, the arguments t and x in Φ∗ are ﬁxed parameters, and we simply write Φ∗ = Φ∗ (ξ) instead of Φ∗ (t, x, ξ). We start with the hyperbolic system (2.5) corresponding to N diﬀerent particle species, but with a slight change of notation which is better adapted for the continuum limit N → ∞. We introduce for simplicity the N equidistant numbers ξi := i/N,

i = 1, . . . , N,

(3.1)

where ξi = δN +1−i , since the ﬁnite sequence of the δi is starting with δ1 = ξN = 1 and decreasing. Adapted to the deﬁnition of the quantities ξi , we will also reverse the numbering of the basic variables φi in the hyperbolic system (2.5) by introducing the new basic variables Φi := φN +1−i . Now we may rewrite (2.5) as ∂t Φi + ∂x Fi = 0,

i = 1, . . . , N,

(3.2)

where the ﬂuxes Fi are given functions of the new basic variables Φi according to Fi := Φi (ξi − (Φ1 ξ1 + · · · + ΦN ξN ))V (Φ1 + · · · + ΦN ),

i = 1, . . . , N.

(3.3)

Φ∗n

Next, we deﬁne the cumulative volume fractions along with the cumulative ∗ ∗ ∗ ﬂuxes Fn by Φn := Φ1 + · · · + Φn and Fn := F1 + · · · + Fn for n = 1, . . . , N . By using the identities n

Φi ξi = Φ∗n ξn −

i=1

n−1

Φ∗i (ξi+1 − ξi ),

n = 1, . . . , N,

(3.4)

i=1

we may express the cumulative ﬂuxes in terms of the cumulative volume fractions: N −1 ∗ ∗ ∗ ∗ Fn = − Φn ΦN ξN − Φi (ξi+1 − ξi )

− Φ∗n ξn −

i=1 n−1 i=1

Φ∗i (ξi+1 − ξi )

V (Φ∗N ),

n = 1, . . . , N.

(3.5)

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1751

From these ﬂuxes, we form the N × N Jacobian with the coeﬃcients ∂Fn∗ ∂Fn∗ ∗ = (Φ , . . . , Φ∗N ), ∂Φ∗k ∂Φ∗k 1 and obtain

∂Fn∗ = − δnk ∂Φ∗k

−

Φ∗N ξN

−

N −1

Φ∗i (ξi+1

δN k ξN − (ξk+1 − ξk )

N

δnk ξn − (ξk+1 − ξk ) Φ∗N ξN −

Φ∗n ξn

−

V (Φ∗N )

δi(n−1)

N −1

Φ∗i (ξi+1 − ξi )

i=1

i=k

− δN k Φ∗n

N

Φ∗n V (Φ∗N )

δi(N −1)

i=k

−

− ξi ) V (Φ∗N )

i=1

+

k, n = 1, . . . , N,

n−1

Φ∗i (ξi+1

V (Φ∗N ).

− ξi )

(3.6)

i=1

Next we prescribe a vector with N real components Ψ∗1 , . . . , Ψ∗N and apply the Jacobian matrix in (3.6) on it to obtain the resulting vector with components

N N −1 ∂Fn∗ ∗ ∗ ∗ ∗ Ψ = − Ψn ΦN ξN − Φi (ξi+1 − ξi ) V (Φ∗N ) ∂Φ∗k k i=1 k=1

N −1 ∗ ∗ − ΨN ξN − (ξk+1 − ξk )Ψk Φ∗n V (Φ∗N ) k=1

+

Ψ∗n ξn

−

(ξk+1 −

ξk )Ψ∗k

− Ψ∗N Φ∗n

Φ∗N ξN −

N −1 i=1

Φ∗n ξn

V (Φ∗N )

k=1

−

n−1

−

n−1

Φ∗i (ξi+1

Φ∗i (ξi+1 − ξi )

− ξi )

V (Φ∗N ).

(3.7)

i=1

The eigenvalue problem for the cumulative N × N hyperbolic system ∂t Φ∗n + ∂x Fn∗ = 0,

n = 1, . . . , N,

(3.8)

with an eigenvalue λ and an eigenvector (Ψ∗1 , . . . , Ψ∗N )T is N ∂F ∗ n

k=1

∂Φ∗k

Ψ∗k = λΨ∗n ,

n = 1, . . . , N,

(3.9)

October 14, 2008 11:9 WSPC/103-M3AS

1752

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

whereas the eigenvalue problem for the original system reads N ∂Fn Ψk = λΨn , ∂Φk

n = 1, . . . , N.

(3.10)

k=1

The relation between (3.9) and (3.10) is explained in the following lemma. Lemma 3.1. For Φ = (Φ1 ,. . . ,ΦN )T ∈ RN , let Ψ = (Ψ1 , . . . , ΨN )T ∈ RN be an eigenvector for the real eigenvalue λ such that (3.10) is satisfied. Define the cumulative quantities Ψ∗n := Ψ1 + · · · + Ψn for n = 1, . . . , N . Then Ψ∗ = (Ψ∗1 , . . . , Ψ∗N )T is eigenvector of (3.9) with the same eigenvalue λ. Proof. We deﬁne the N × N Jacobians ∇∗ F ∗ and ∇F with the respective coeﬃcients ∂Fn∗ ∂Fn∗ ∗ ∂Fn ∂Fn = (Φ , . . . , Φ∗N ), = (Φ1 , . . . , ΦN ), k, n = 1, . . . , N, ∗ ∂Φk ∂Φ∗k 1 ∂Φk ∂Φk and the regular N × N matrix

TN

  1 0 ··· 0 .  1 1 . . . ..  , = .  . . . . . 0 1 ··· ··· 1

such that Ψ∗ = TN Ψ. Then we have ∗ −1 −1 ∗ −1 (∇F )Ψ = λΨ ⇔ (∇F )T−1 N (TN Ψ) = TN (λTN Ψ) ⇔ (∇F )TN Ψ = TN (λΨ ) ∗ ∗ ∗ ∗ ∗ ⇔ (TN (∇F )T−1 N )Ψ = λΨ ⇔ (∇∗ F )Ψ = λΨ ,

which is the statement of the lemma. Corollary 3.1. Consider a vector Φ∗ = (Φ∗1 , . . . , Φ∗N )T ∈ D∗ := {Φ∗ = TN Φ | Φ ∈ D} and define ξ¯N := ξ1 Φ1 + · · · + ξN ΦN ,

λ− := −2ξ¯N V (Φ∗N ) + V (Φ∗N )(ξ¯N + Φ∗N ).

(3.11)

If V (Φ∗N ) > 0, then the eigenvalues λ1 (Φ∗ ) ≥ λ2 (Φ∗ ) ≥ · · · ≥ λN (Φ∗ ) of the eigenvalue problem (3.9) satisfy λN +1−i (Φ∗ ) ∈ [ξi−1 , ξi ], i = 2, . . . , N, ξ¯N + V (Φ∗N ) ∗ V (ΦN ) ¯N λN (Φ∗ ) λ− N N ∗ N ¯ ¯ ¯ ∈ ξ + , ξ1 = (ξ + ΦN ) − ξ , ξ1 . ξ + V (Φ∗N ) V (Φ∗N ) V (Φ∗N )

(3.12) (3.13)

If V (Φ∗N ) = 0, then we have that λi (Φ∗ ) = 0 for i = 1, . . . , N − 1, and λN (Φ∗ ) ∈ [V (Φ∗N )(ξ¯N + φ∗N ), 0].

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1753

Proof. The corollary is a consequence of Theorem 2.1, Corollary 2.1 and Lemma 3.1. 3.2. Convergence of the linearized operators for N → ∞ and spectrum of the “continuous” operator Lemma 3.1 motivates why we consider the cumulative problem (3.9) instead of the original problem (3.10), because in the limit N → ∞ the cumulative eigenvectors will be replaced by integrals. We ﬁrst state precisely in which sense the ﬁnitedimensional operators converge to a “continuous” operator. For a given vector Φ∗ = (Φ∗1 , . . . , Φ∗N )T ∈ [0, φmax ]N , we deﬁne the discrete (ﬁnite-dimensional) linear operator A˜N (Φ∗ ) : RN → RN by (A˜N (Φ∗ )Ψ∗ )n =

N ∂F ∗ n

k=1

∂Φ∗k

(Φ∗ )Ψ∗k ,

n = 1, . . . , N,

Ψ∗ = (Ψ∗1 , . . . , Ψ∗N )T ∈ RN . (3.14)

We associate with A˜N the “continuous” operator AN (Φ∗ ) : F → L1 [0, 1],

AN (Φ∗ )Ψ∗ = (PN ◦ A˜N (Φ∗ ) ◦ ΠN )(Ψ∗ ),

(3.15)

where we deﬁne F := BV [0, 1] ⊂ L1 [0, 1] associated with the norm f F := f 1 + TV(f ), and the projection operator ΠN : L1 [0, 1] → RN , which maps a function Ψ∗ = Ψ∗ (ξ) onto its cell averages with respect to the partition ZN := {ξi = i/N, i = 0, . . . , N }, such that ξi ∗ Ψ∗ (ξ) dξ, i = 1, . . . , N ; ΠN Ψ∗ 1 = N Ψ∗ L1 [0,1] , (ΠN Ψ )i = N ξi−1

and the inverse of ΠN , the prolongation operator PN : RN → L1 [0, 1], which assigns to a vector y = (y1 , . . . , yN )T ∈ RN the function PN y deﬁned by (PN y)(ξ) =

N

χ[ξi−1 ,ξi ) (ξ)yi .

i=1

In other words, taking into account that ξi+1 − ξi = 1/N for i = 0, . . . , N − 1, AN (Φ∗ ) maps Ψ∗ ∈ F onto the piecewise constant function AN (Φ∗ )Ψ∗ deﬁned by N ∂Fi∗ ∗ (AN (Φ )Ψ )(ξ) = χ[ξi−1 ,ξi ) (ξ) Ψ . ∂Φ∗m m m=1 i=1 ∗

∗

N

(3.16)

We deﬁne an operator norm AN (Φ∗ )1 :=

sup

Ψ∗ 1 ≤1

AN (Φ∗ )Ψ∗ 1 .

In what follows, we will compare operators AN (Φ∗ ) with diﬀerent values of N , and in the notation, we replace Φ∗ by Φ∗,N if Φ∗ ∈ RN . Speciﬁcally, we assume that Φ∗,N is the vector of cumulative concentration, deﬁned with respect to ZN , of a

October 14, 2008 11:9 WSPC/103-M3AS

1754

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

non-negative function Φ ∈ L1 (0, 1), such that Φ1 ≤ φmax , i.e. ξi i ξj ∗,N ∗,N T ∗,N ∗,N Φ = (Φ1 , . . . , ΦN ) , Φi = Φ(ξ)dξ = Φ(ξ) dξ, j=1

ξj−1

0

(3.17)

i = 1, . . . , N. Lemma 3.2. For a fixed function Φ ∈ L1 [0, 1], the sequence of operators {AN (Φ∗,N )}N ∈N is a Cauchy sequence, which converges in · 1 to a bounded linear operator A = A[Φ∗ ] : A[Φ∗ ] : L1 [0, 1] ⊃ F Ψ∗ → (A[Φ∗ ])Ψ∗ ∈ L1 [0, 1].

(3.18)

Note that A = A[Φ∗ ] denotes a functional dependence. The operator A depends on the limit of vectors Φ∗,N for N → ∞, i.e. on the primitive of Φ, the cumulative phase density deﬁned by ξ ˜ dξ, ˜ ξ ∈ [0, 1]. Φ(ξ) (3.19) Φ∗ (ξ) := 0

1

Theorem 3.1. Let Φ ∈ L [0, 1] be a given phase density function and Φ∗ ∈ C 0 [0, 1] ¯ = be the associated cumulative phase density (3.19), where it is assumed that Φ ∗ ∗ Φ (1) ≤ φmax . Then the operator A = A[Φ ] has a bounded spectrum Λ(A) = Λ− (A) ∪ Λ+ (A),

(3.20)

¯ := [V (Φ)( ¯ −V (Φ) ¯ ¯ ξ) ¯ ξ¯ + Φ) ¯ − 2V (Φ) ¯ ξ, ¯ ξ), Λ− (A) ⊂ I1 (Φ, 1 ¯− ξ¯ := Φ Φ∗ (ξ) dξ,

(3.21)

¯ := [−V (Φ) ¯ V (Φ)(1 ¯ ¯ ξ) ¯ ξ, ¯ Λ+ (A) = I(Φ, − ξ)].

(3.23)

where

(3.22)

0

Proof. We select Φ ∈ L1 [0, 1] and assume that Φ∗ , its primitive, satisﬁes Φ∗ (1) < φmax . Assume that Φ∗,N ∈ RN is determined by (3.17), and let the operators A˜N (Φ∗,N ) and AN (Φ∗,N ) be deﬁned by (3.14) and (3.15), respectively. To emphasize the dependence on N , we denote the eigenvalues of A˜N (Φ∗,N ) by ∗,N N ), i = 1, . . . , N , where λN λN 1 ≥ · · · ≥ λN , and recall that these eigenvali = λi (Φ ues satisfy the inclusions given by Corollary 3.1. If we use that ¯ ¯ ≤ Φ1 = Φ ≤ 1 , (3.24) |ξ¯N − ξ| N N N ¯ is the same for all values of N , then (3.11)–(3.13) imply along with the fact that Φ that

i − 1 ¯N i N N ¯ ¯ ¯ λN +1−i ∈ V (Φ) −ξ −ξ , V (Φ) N N

i−2 ¯ i+1 ¯ N ¯ ¯ ¯ ¯ − ξ , V (Φ) − ξ , i = 2, . . . , N, ⊂ Ii (Φ, ξ) := V (Φ) N N (3.25)

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1 N N N ¯ ¯ ¯ ¯ ¯ ¯ ¯ −ξ ∈ V (Φ)(ξ + Φ) − 2V (Φ)ξ , V (Φ) N

2 1 1 N ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ −ξ . ⊂ I1 (Φ, ξ) := V (Φ) ξ + Φ + − 2V (Φ) ξ + , V (Φ) N N N (3.26)

λN N

1755

Observing that

N ∞ N ¯ ¯ ¯ ξ) ¯ ∪ I(Φ, ¯ ξ) ¯ Ii (Φ, ξ) = I1 (Φ, N =1

i=1

¯ V (Φ)(1 ¯ =: I( ¯ ¯ ξ¯ + Φ ¯ − 2V (Φ) ¯ ξ, ¯ ˜ Φ, ¯ ξ), = [V (Φ) − ξ)]

which is a compact interval, we conclude that every sequence {λN i(N ) }N ∈N of eigen∗,N values (of the corresponding sequence of operators AN (Φ )) has a subsequence, ¯ for ˜ Φ, ¯ ξ) which we do not bother to relabel, that converges to a number λ ∈ I( ∗,N N N N → ∞. Associated with an eigenvalue λi is an eigenvector Ψi ∈ R , or a piecewise constant eigenfunction Ψ∗,N (ξ). We assume the normalization ΨN i F = 1, so i ∗,N {Ψi(N ) }N ∈N is the sequence of corresponding (normalized) eigenfunctions. Since these eigenfunctions are normalized and uniformly bounded in L1 [0, 1], we may select a subsequence which converges in L1 [0, 1] to a function Ψ∗ with Ψ∗ F = 1. We will not bother to relabel this sequence, so we assume that ∗ Ψ∗,N i(N ) − Ψ 1 → 0

as N → ∞.

To see that Ψ∗ coincides a.e. with an eigenfunction of the operator A, we note that ∗ (A − λI)Ψ∗ 1 ≤ (A − AN (Φ∗,N ))Ψ∗ 1 + (AN (Φ∗,N ) − λN i(N ) )Ψ 1 ∗ + |λN i(N ) − λ|Ψ 1 ∗,N ∗ = (A − AN (Φ∗,N ))Ψ∗ 1 + (AN (Φ∗,N ) − λN i(N ) )(Ψ − Ψi(N ) )1 ∗ + |λN i(N ) − λ|Ψ 1

≤ (A − AN (Φ∗,N ))Ψ∗ 1 + (AN (Φ∗,N ) − λ)(Ψ∗ − Ψ∗,N i(N ) )1 ∗,N ∗ ∗ + |λN i(N ) − λ|(Ψ 1 + Ψ − Ψi(N ) 1 ).

(3.27)

Since the right-hand side of (3.27) converges to zero for N → ∞, while the left-hand side does not depend on N , we conclude that AΨ∗ − λΨ∗ = 0

a.e., while Ψ∗ = 0,

(3.28)

so λ is an eigenvalue of A, and Ψ∗ is the corresponding eigenfunction. ¯ contains no eigenvalue of A, let us assume that λ is ˜ Φ, ¯ ξ) To show that C\I( an eigenvalue of A, with a normalized eigenfunction Ψ∗ , Ψ∗ F = 1, such that AΨ∗ − λΨ∗ 1 = 0. On the other hand, let {λN i(N ) }N ∈N be an arbitrary sequence of

October 14, 2008 11:9 WSPC/103-M3AS

1756

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

eigenvalues of the corresponding operators AN (Φ∗,N ) with the normalized eigenfunctions Ψ∗,N i(N ) . We then have that 0 = λΨ∗ − AΨ∗ ∗,N ∗ = (λ − λN ))Ψ∗ i(N ) )Ψ − (A − AN (Φ ∗,N ∗,N ∗,N ∗ − (AN (Φ∗,N ) − λN ) − λN i(N ) )Ψi(N ) − (AN (Φ i(N ) )(Ψ − Ψi(N ) ).

(3.29)

Now we select a subsequence of {Ψ∗,N i(N ) }N ∈N , which we do not relabel, which converges in L1 to Ψ∗ . Since the second and fourth terms on the right-hand side of (3.29) vanish in L1 [0, 1] as N → ∞, and the third term is zero, we conclude that for a given constant δ > 0, we may choose N0 ∈ N such that ∗ N ∗ δ ≥ (λ − λN i(N ) )Ψ = |λ − λi(N ) |Ψ for N ≥ N0 .

Since Ψ∗ > 0, this implies that the corresponding sequence {λN i(N ) }N ∈N converges ¯ ˜ ¯ to λ for N → ∞, which implies that λ ∈ I(Φ, ξ), i.e. all eigenvalues of A[Φ∗ ] belong ¯ ˜ Φ, ¯ ξ). to I( ¯ con¯ ξ) Finally, to prove (3.20), we recall that each of the intervals IiN (Φ, ¯ ξ) ¯ can be approximated tains at least one eigenvalue, such that each λ ∈ I(Φ, N N ¯ ¯ N ¯ ¯ by a sequence of eigenvalues {λN i(N ) }N ∈N , λi(N ) ∈ I2 (Φ, ξ) ∪ . . . ∪ IN (Φ, ξ), and the proof of (3.28) ensures the existence of a corresponding eigenfunction. Conse¯ is indeed an eigenvalue of A. ¯ ξ) quently, every λ ∈ I(Φ, Remark 3.1. The smallest eigenvalues {λN N }N ∈N are bounded in the sequence ¯ but these intervals do not shrink as N → ∞, so it is not ¯ ξ), of intervals I1N (Φ, clear at present whether the smallest eigenvalues converge to an eigenvalue of A. Though numerical experiments support this conjecture, a rigorous proof would involve additional arguments to those of the proof of Theorem 3.1, and is beyond the scope of this paper. 3.3. The eigenvalue problem under smoothness assumptions Even though the preceding analysis implies that eigenfunctions of the cumulative, continuous problem will in general not be smooth, which, for example, occurs if Φ represents a distribution corresponding to a ﬁnite number of sizes, it is still of interest to analyze the eigenvalue problem under the assumption that Φ∗ : [0, 1] → [0, φmax ] is a continuously diﬀerentiable and monotone function such that Φ∗ (0) = 0 and Φ∗ (1) ≤ φmax . We require that ¯ = 0 V (Φ)

¯ < φmax . and Φ

(3.30)

The ﬁrst condition in (3.30) implies the second due to the constitutive law (2.3) for V˜ . Furthermore, for Ψ∗ ∈ C 1 [0, 1] we may consider the operator A = A[Φ∗ ] arising in the limit N → ∞ as an integral operator A : C 1 [0, 1] → C 1 [0, 1], which

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1757

according to (3.7) is given by ∗ ¯ ¯ (AΨ )(ξ) := − Ψ (ξ)ξV (Φ) − Ψ (1) − ∗

∗

+ Ψ∗ (ξ)ξ −

1

∗

¯ Ψ (η) dη Φ∗ (ξ)V (Φ)

0

ξ

¯ Ψ∗ (η) dη V (Φ)

0

− Ψ∗ (1) Φ∗ (ξ)ξ¯ − Φ∗ (ξ)ξ −

ξ

Φ∗ (η)dη

¯ V (Φ).

(3.31)

0

The analogue of the discrete eigenvalue problem (3.9) will now be considered in the limit N → ∞ and reads AΨ∗ = λΨ∗ , where we deﬁne ξλ := ξ¯ +

λ ¯ V (Φ)

(3.32)

and have to require for all ξ ∈ [0, 1] that ξ ∗ ∗ (ξ − ξλ )Ψ (ξ) = Ψ (η) dη + Ψ∗ (1) − 0

1

∗

Ψ (η) dη Φ∗ (ξ)

0

+ Ψ∗ (1)

ξ ¯ V (Φ) ∗ ∗ ¯ − ξ) + (ξ)( ξ Φ (η) dη . Φ ¯ V (Φ) 0

(3.33)

We ﬁrst deal with the characterization of the eigenvalues of the operator A for a given function Φ(ξ) or equivalently, a given cumulative phase density function Φ∗ (ξ), and then turn to the corresponding eigenfunctions. For the solution we distinguish two cases. Case A: Ψ∗ (1) = 0. This case is important for the discussion of the uniqueness problem of the eigenfunctions corresponding to a given eigenvalue λ. We put ξ := 1 in (3.33) to obtain ﬁrst with (3.30) 1 Ψ∗ (η) dη = 0. 0

Here the eigenvalue problem reduces to ∗

(ξ − ξλ )Ψ (ξ) =

ξ

Ψ∗ (η) dη.

(3.34)

0

We ﬁrst assume that ξλ ∈ / [0, 1]. Then (3.34) implies that Ψ∗ is arbitrarily smooth with derivative zero, such that Ψ∗ (ξ) = 0 for all ξ ∈ [0, 1]. Next, we assume that ξλ ∈ [0, 1] and obtain from (3.34) in the limit ξ → ξλ ξλ Ψ∗ (η) dη = 0, 0

October 14, 2008 11:9 WSPC/103-M3AS

1758

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

which implies ∗

(ξ − ξλ )Ψ (ξ) =

ξ

Ψ∗ (η) dη.

ξλ

In particular, in this subcase Ψ∗ is arbitrarily smooth with zero derivative and hence zero for all ξ ∈ [0, 1]. We ﬁnally conclude that in Case A we have only the zero- solution for Ψ∗ , which is not an eigensolution. Case B: Ψ∗ (1) = 1. This is a normalization condition for the eigenfunctions which can be assumed without loss of generality if Case A does not hold. Assume ﬁrst that we have two eigenfunctions Ψ∗1 and Ψ∗2 with the same real eigenvalue λ, which both satisfy the normalization condition Ψ∗1 (1) = Ψ∗2 (1) = 1. Then we conclude that Ψ∗ = Ψ∗1 − Ψ∗2 satisﬁes (3.34) in Case A with Ψ∗ (1) = 0, such that Ψ∗ ≡ 0 and hence Ψ∗1 = Ψ∗2 . This means that in Case B, the eigenvalue problem admits at most one eigenfunction with respect to a ﬁxed eigenvalue λ. If we put ξ := 1 in (3.33), then we can calculate from (3.22) that 1 ¯ ξλ ¯V (Φ) . (3.35) − ξ Ψ∗ (η) dη = 1− ¯ ¯ 1−Φ V (Φ) 0 Thus the eigenvalue problem may be rewritten in the form ξ (ξ − ξλ )Ψ∗ (ξ) = Ψ∗ (η) dη + F (ξ) 0 ξ ¯ ξλ V (Φ) ∗ ∗ ∗ with F (ξ) = Φ (ξ) − (ξ) − Φ (η) dη . ξΦ ¯ ¯ 1−Φ V (Φ) 0

(3.36)

But F is continuously diﬀerentiable on ξ ∈ (0, 1) due to our assumptions, and therefore so is the expression Ψ∗ (ξ) on the left-hand side of (3.36), except at most for ξ = ξλ ∈ (0, 1). Thus we deﬁne at least for ξ ∈ (0, 1), ξ = ξλ , the quantities Φ(ξ) :=

dΦ∗ (ξ), dξ

Ψ(ξ) :=

dΨ∗ (ξ), dξ

(3.37)

and obtain from the ﬁrst equation in (3.36) for ξ ∈ (0, 1), ξ = ξλ , the necessary condition ¯ ξλ V (Φ) Φ(ξ) − (3.38) Ψ(ξ) = ¯ ¯ ξ . ξ − ξλ 1 − Φ V (Φ) In order to obtain a solution of the eigenvalue problem, also the condition (3.35) has to be satisﬁed. This requires that the right-hand side in (3.38) is continuous for 0 ≤ ξ ≤ 1 and hence integrable. We summarize our results and interpret Eq. (3.38) as an eigensolution of the original eigenvalue problem (3.10) in the continuum limit

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1759

N → ∞, provided that (a) in the case ξλ ∈ [0, 1] there exists for ξ ∈ [0, 1] \ {ξλ } the limit ¯ ξλ V (Φ) Φ(ξ) Ψ(ξλ ) := lim ¯ − V (Φ) ¯ ξ , ξ→ξλ ξ − ξλ 1 − Φ

(3.39)

(b) the integration constant C ∈ R of the primitive function Ψ∗ (ξ) := C +

ξ

0

¯ ξλ V (Φ) Φ(η) − ¯ ¯ η dη, η − ξλ 1 − Φ V (Φ)

ξ ∈ [0, 1],

(3.40)

is chosen such that the normalization condition Ψ∗ (1) = 1 is satisﬁed. If Ψ∗ is a normalized eigensolution of the eigenvalue problem (3.33) with eigenvalue ξλ , then it is necessarily given by (3.40). For ξλ < 0, an eigenfunction Ψ∗ may be readily determined from (3.40). This occurs, for example, in rarefaction-type solutions of Riemann problems with upwards-propagating solution information; a numerical example of such a case is given by Example 5 in Sec. 5.2. On the other hand, for ξλ ∈ [0, 1] satisfaction of (3.39) ensures that we can evaluate (3.40). Obviously, (3.39) is satisﬁed if ¯ V (Φ) 1 = ¯ ¯ , 1−Φ V (Φ)

(3.41)

or Φ(ξλ ) = 0. We performed numerical tests (see Sec. 5.2) to decide whether the discrete analogues of ξλ , namely the transformed eigenvalues ξλ1 , . . . , ξλN deﬁned by N

λi λi (Φ) = , ξλi = ξ¯ + ξi Φi + V (φ) V (Φ + · · · + ΦN ) 1 i=1

i = 1, . . . , N,

(3.42)

belong to the interval [0, 1] or not. In agreement with Corollary 3.1, it turns out that almost all discrete values for ξλ belong to [0, 1], with the possible exception of the smallest value.

4. Numerical Schemes 4.1. Discretization parameters and boundary conditions We discretize time by setting tm := mδt for m = 0, . . . , M , where δt := T /M with M ∈ N, and space by xk := kδx for k = 0, . . . , K, where δx := L/K and K ∈ N is given. We set ξn := nδξ for n = 0, . . . , N , where N ∈ N and δξ := 1/N . The numbers M , K and N are given discretization parameters. Finally, we denote by φm k,n the ﬁnite diﬀerence approximation of Φ(tm , xk−1/2 , ξ) for ξn−1 ≤ ξ < ξn . At a ﬁxed time t = tm , we discretize Φ = Φ(tm , x, ξ) in the x-ξ phase space by the

October 14, 2008 11:9 WSPC/103-M3AS

1760

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

piecewise constant function ΦD (tm , x, ξ), where ΦD (tm , x, ξ) = φm k,n on the phase space cell Ωk,n := {(x, ξ) | xk−1 ≤ x < xk , ξn−1 ≤ ξ < ξn },

k = 1, . . . , K,

n = 1, . . . , N .

Thus, ΦD (tm , ·, ·) is given by a K × N matrix. However, for the numerical computation, we only store the updated current initial data, starting at the current time step with the initial matrix (Φ0k,n ) ∈ RK×N . Boundary conditions for a kinetic scheme are in general a subtle issue. Below we discuss two cases of numerical boundary conditions, which can be described suﬃciently in every time step tm by two given vectors m m m Φm − := (φ0,1 , φ0,2 , . . . , φ0,N ),

m m m N Φm + := (φK+1,1 , φK+1,2 , . . . , φK+1,N ) ∈ R .

They describe essentially initial value problems, but from the numerical point of view they are still important, since the computation can only be performed on a ﬁnite domain. For Riemann initial data, which provide the simplest test case of interest, we can choose the boundary data as m φm 0,n = φ1,n ,

m φm K+1,n = φK,n ,

n = 1, . . . , N .

(4.1)

This is suﬃcient for the numerical study of certain Riemann initial value problems on the spatial computational domain [0, L], provided that up to the ﬁnal time t = T , no signal coming from the interior of the spatial domain has reached the boundaries. We deﬁne the boundary phase cells Ω0,n := {(x, ξ) |−δx ≤ x < 0, ξn−1 ≤ ξ < ξn }, ΩK+1,n := {(x, ξ) | L ≤ x < L + δx , ξn−1 ≤ ξ < ξn }, on which ΦD (tm , x, ξ) assumes the piecewise constant boundary values of (4.1). Thus, ΦD (tm , x, ξ) is given on the extended domain −δx ≤ x ≤ L + δx , 0 ≤ ξ ≤ 1 in a piecewise constant manner by the extended matrix (φm k,n )k=0,...,K+1,n=1,...,N . Finally, we deﬁne for k = −1, . . . , K the quantities xk+1/2 := (k + 1/2)δx . For k = 0, . . . , K and n = 1, . . . , N , we deﬁne the intervals J˜k := [xk−1/2 , xk+1/2 ] and ˜ k,n := J˜k × [ξn−1 , ξn ]. The latter the rectangles, also-called dual phase space cells, Ω ˜ k,n | = δx δξ , as have the cells Ωk,n . have the area |Ω However, for the sedimentation process in a closed vessel we require a discrete version of the zero-ﬂux boundary condition (2.18), which cannot be captured appropriately by this approach. Rather, this boundary condition has to supplement the discretization of the equation of the generalized kinetic theory, which will be described ﬁrst. 4.2. Discretization of the equation of the generalized kinetic theory Using (2.16), we can exactly calculate at time t = tm the time derivative of the integral average 1 ΦD (t, x, ξ) dx dξ, k = 0, . . . , K, n = 1, . . . , N . (4.2) φ˜k,n (t) := δx δξ ˜ k,n Ω

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1761

To this end, we consider for ξn−1 < ξ < ξn the relationships ΦD (tm , xk−1/2 , ξ) = φm k,n ,

k = 0, . . . , K + 1,

n = 1, . . . , N .

(4.3)

Using the following deﬁnitions, which are valid for every x with xk−1 < x < xk , 1 N := Φ (t , x, ξ) dξ = δ φm Am D m ξ k k,n , 0

Bkm :=

0

n=1 1

ξΦD (tm , x, ξ) dξ = δξ

N n=1

φm k,n

ξn + ξn−1 , 2

we obtain by taking the integral average of (2.16) with respect to the phase space ˜ k,n : cell Ω ξn dφ˜k,n 1 m (tm ) + (φm V (Am k+1 )(ξ − Bk+1 ) dt δx δξ ξn−1 k+1,n m m − φm k,n V (Ak )(ξ − Bk )) dξ = 0,

(4.4)

which implies δx

ξn + ξn−1 dφ˜k,n m m (tm ) =−φm − B V (A ) k+1,n k+1 k+1 dt 2

ξn + ξn−1 m − Bkm , k = 0, . . . , K, + φm k,n V (Ak ) 2

n = 1, . . . , N .

Note that for the integration with respect to x, we have used here the fundamental theorem of inﬁnitesimal calculus, and the integration with respect to ξ has been performed exactly as well, so conservativity is ensured. Consequently, for the time evolution, the ﬁrst proposed numerical scheme (Scheme 1) is the following:

dφ˜k−1,n 1 m δt dφ˜k,n m (φ + = + φ ) − φm+1 k−1,n k,n 2 k+1,n 2 dt dt

1 1 δt m ξn + ξn−1 m m m − B = (φm + φ ) − V (A ) φ k−1,n k+1 k+1 2 k+1,n 2 δx k+1,n 2

ξn + ξn−1 m m − B − φm V (A ) , k = 1, . . . , K, n = 1, . . . , N . k−1,n k−1 k−1 2 If we use the following marching formula in conservation form: m m ˜ m φm+1 k,n = φk,n − λ(hk+1/2,n − hk−1/2,n ),

˜ := δt /δx , λ

k = 1, . . . , K,

n = 1, . . . , N ,

(4.5)

where hm k+1/2,n is the numerical ﬂux, then Scheme 1 represents the well-known conservative ﬁrst-order Lax–Friedrichs central diﬀerence scheme with numerical ﬂuxes 1 m 1 ¯ m m m m ¯ m hm k+1/2,n = hn (Φk , Φk+1 ) = ˜ (φk,n − φk+1,n ) + (fn (Φk ) + fn (Φk+1 )), 2 2λ k = 0, . . . , K, n = 1, . . . , N ,

October 14, 2008 11:9 WSPC/103-M3AS

1762

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

m m T m m ¯ m with Φm k := (φk,1 , . . . , φk,N ) and fn (Φk ) := φk,n vn (Φk ), where we deﬁne

ξn + ξn−1 m m m vn (Φk ) := V (Ak ) − Bk . (4.6) 2 m T ¯ m ¯ We also deﬁne the vector ¯f (Φm k ) := (f1 (Φk ), . . . , fN (Φk )) . This scheme is stable provided the following CFL condition holds

˜ max ρ(J¯ (Φm )) ≤ 1, λ f k

(4.7)

¯ m where ρ(·) denotes the spectral radius, and J¯f (Φm k ) the N × N Jacobian of f (Φk ). m We approximate max ρ(J¯f (Φk )) by α := max |vn (Φm k )|. k,m,n

The second proposed numerical scheme (Scheme 2), which is introduced in Ref. 19, is given by (4.5) with the following numerical ﬂuxes, where vn (Φm k ) is deﬁned by (4.6): m m m m hm k+1/2,n = φk,n max{0, vn (Φk+1 )} + φk+1,n min{0, vn (Φk+1 )},

k = 0, . . . , K,

n = 1, . . . , N .

(4.8)

If negative velocities are present, which is possible in our model, Scheme 2 produces sharply resolved interfaces, but with overshoots or oscillations in some situations.19 To overcome this shortcoming we propose a more viscous version of Scheme 2 that provides a good compromise between sharply resolved interfaces and suppression of overshoots. Scheme 3 is the conservative formula (4.5) with numerical ﬂuxes hm k+1/2,n =

m Ek+1 1 m m m m {φk+1,n vn (Φm (φm ) + φ v (Φ )} − n k+1 k,n k k+1,n − φk,n ) 2 2 φm k,n m m m |vn (Φm − k ) − vn (Φk+1 )|sgn(φk+1,n − φk,n ), 2 k = 0, . . . , K, n = 1, . . . , N ,

m m := max{|v1 (Φm where Ek+1 k+1 )|, . . . , |vN (Φk+1 )|}, and which is also introduced in Ref. 19. For Schemes 2 and 3, according to Ref. 19, there is not a proof of stability, but based on an analysis of a model with non-negative velocities, the following CFL condition is given: ˜ max |vn (Φm )| ≤ 1 . (4.9) λ k k,m,n 2 It is easy to prove that for our model, the estimate maxk,m,n |vn (Φm k )| = 1. m and φ , we utilize (4.1) for the Riemann initial value probTo compute φm 0,n K+1,n lem. However, as mentioned above, for the description of a sedimentation process in a closed vessel, we have to formulate zero-ﬂux boundary conditions as a discrete version of (2.18). To this end, we modify the scheme by setting m hm 1/2,n = hK+1/2,n = 0 ,

in the conservative formula (4.5).

m = 0, . . . , M,

n = 1, . . . , N

(4.10)

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1763

4.3. Numerical test of mass conservation We compute the original phase density Φ in terms of the matrix with coeﬃcients φm k,n according to the numerical scheme described above. In order to verify the numerical mass conservation for a sedimentation process in a closed vessel, we only have to check that the “species vector” with components m Snm := φm 1,n + · · · + φK,n ,

n = 1, . . . , N ,

(4.11)

is conserved in time to machine accuracy, i.e. is independent of the time step tm . If we deﬁne at the current time step tm the cumulative K × N phase density matrix with coeﬃcients m m φ∗m k,n := φk,1 + · · · + φk,n ,

(4.12)

then we can also compute at each time step tm the cumulative N -vector with coeﬃcients Sn∗m :=

K k=1

φ∗m k,n =

n K

φm k,ν ,

n = 1, . . . , N .

(4.13)

k=1 ν=1

Instead of using (4.11), we can also check numerically whether this cumulative vector is independent on the time step tm in order to justify mass conservation for the batch sedimentation process. Of course, for the implementation we need not calculate and store the cumulative matrix with the coeﬃcients given in (4.12) because we can directly calculate the double sum on the right-hand side of (4.13). The conservation property of Skm may now be interpreted as the fact that none of the species will exchange mass with other species or with the ﬂuid, so that the total mass of each species is conserved during the sedimentation process. 5. Numerical Examples The equi-distribution of particle size is represented by the distribution function Feq (ξ) := ξ for ξ ∈ [0, 1]. On the other hand, the Rosin–Rammler distribution function45 is given by FRR (ξ) := 1 − exp(−( ξd1 /l)m ), ξ ∈ [0, 1], where l is a characteristic size, m is a uniformity coeﬃcient, and d1 is the diameter of the largest species. Since FRR (1) < 1, we use the normalized version FRRn (ξ) := FRR (ξ)/FRR (1) for ξ ∈ [0, 1]. If φtot (x) denotes the initial proﬁle of total solids volume fraction, then the initial solids phase density of a suspension with particle size distribution F (ξ) is Φ0 (x, ξ) = φtot (x)F (ξ),

x ∈ [0, 1],

ξ ∈ [0, 1],

F (ξ) :=

d F (ξ). dξ

(5.1)

In all examples, we employ V˜ (φ) given by (2.3) with γ = 2.7, φmax = 0.8 and φq = φ˜q = 0.6824, so the transition between the linear and nonlinear portions of V˜ (φ) is smooth.

October 14, 2008 11:9 WSPC/103-M3AS

1764

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

5.1. Examples with zero-flux boundary conditions 5.1.1. Example 1: Sedimentation of a suspension with equi-distribution of particle size As an example with zero-ﬂux boundary conditions, we simulate the sedimentation of a suspension with equi-distributed particle sizes in a closed vessel. The initial condition is (5.1) with φtot (x) ≡ 0.1, and the solution is calculated by Scheme 3 with ˜ = 0.5 and δξ = 1/128. (Note that in our dimensionless variables, δx = 1/3600, λ we measure time in multiples of the time a particle of the largest species needs to settle through the vessel in pure ﬂuid.) Figure 1 shows the simulated phase density Φ as a function of ξ and x at three diﬀerent times. Figures 1(a) and (b) illustrate that larger particles settle faster than smaller ones, and therefore larger ones ﬁll the bottom of the vessel and form a sediment. We also see that the smaller particles are partially removed from the bottom and, combined with those that were initially located in upper regions of the suspension, form a very thin and concentrated sediment layer above the larger ones.

(a)

(b)

(c) Fig. 1.

Example 1: Simulated phase density Φ at (a) t = 1, (b) t = 5, (c) t = 20.

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1765

The formation of this sediment layer is explained by the form of the ﬂux function (2.19). We observe for ﬁxed x that if the normalized squared size ξ of a particle ¯ then the particle will move upwards due to the is smaller than the average size ξ, negative value of the ﬂux, yielding an accumulation of particles of the same size in the upper zone of the suspension. Figure 2 shows the simulated phase density Φ as a function of ξ at four diﬀerent positions and three diﬀerent times. Figures 2(a) and (b) reveal that at x = 0.5 ¯ the average normalized squared size ξ, ¯ and x = 0.75 the total volume fraction Φ, and the width of the size range, which we denote by wr , decrease as time advances. Moreover, ξ¯ is on the right half of the size range. In Fig. 2(c) we see that at x = 0.9, ¯ is not evident, both ξ¯ and wr decrease as time advances, while the evolution of Φ ¯ increases. Likewise, ξ¯ is on the right half of the size but Fig. 5 indicates that Φ ¯ ξ¯ and wr do not range. According to Fig. 2(d), at the bottom of the vessel Φ, ¯ change, and ξ is on the right half of the size range. Figure 3 shows the simulated phase density Φ as a function of x for three values of ξ at three diﬀerent times. In Fig. 3(a) we see that for the smallest species the peak of concentration moves upwards. In Figs. 3(b) and (c) we observe that the

(a)

(b)

(c)

(d)

Fig. 2. Example 1: Simulated phase density Φ as a function of ξ for three diﬀerent times. (a) x = 0.5, (b) x = 0.75, (c) x = 0.9, (d) x = 1.0.

October 14, 2008 11:9 WSPC/103-M3AS

1766

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

(a)

(b)

(c)

(d)

Fig. 3. Example 1: Simulated phase density Φ as a function of x for three diﬀerent times: (a) ξ = 0.0039602, (b) ξ = 0.33203, (c) ξ = 0.66797, (d) ξ = 0.99609.

peak of concentration for ξ = 0.33203 and ξ = 0.66797 moves downwards, and at t = 1 there is a peak of concentration in the upper part of the suspension. Figure 3 indicates that for the largest species there is a peak of concentration at the bottom, and at t = 1 there is a smooth change of concentration at the upper part of the suspension. Figure 4 shows the simulated total solids volume fraction until t = 20 and a zoom on the same solution for t ∈ [0, 4]. In this example, we also record approximate L1 errors deﬁned with respect to a reference solution, and convergence rates to evaluate the performance of the numerical schemes. The approximate L1 errors e1 and e2 are deﬁned by e1 := δξ δx

µ MR N n=1 j=1 k=ML

m |φ˜m µ(n−1)+j,k − φn,k |,

M µ N R m e2 := δξ δx (φ˜m µ(n−1)+j,k − φn,k ) , n=1 j=1 k=ML

m where φ˜m ˜ and the approximate n ˜ ,l and φn,l are the reference solution for ξ = ξn solution for ξ = ξn , respectively, both at x = xl and t = tm ; µ is the value of N of

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

(a) Fig. 4.

1767

(b) Example 1: Simulated total solids volume fraction φ.

the reference solution divided by that of the approximate solution; ML and MR are the indices of the positions between which we calculate the errors of the numerical approximation; and δ˜ξ is the size discretization parameter of the reference solution. The reference solution is calculated with Scheme 3 with δ˜ξ = 1/128, and the other solutions for the error study are calculated with δξ ∈ {1/2, 1/4, 1/8, 1/16, 1/32}. For the reference solution and all other computations of the error study, we use ˜ = 0.5. δx = 1/3600 and λ Figures 5 and 6 show the simulated total solids volume fraction φ as a function of x for t = 5 and t = 20, and the simulated phase density Φ as a function of ξ for t = 1, t = 5 and t = 20, respectively, produced by Schemes 1, 2 and 3, while Tables 1 and 2 list the approximate L1 errors for Φ and φ, respectively, measured over the domain [0, 1] × [0, 1] at t = 1 and t = 20. Figure 5 shows that Scheme 2 introduces less numerical diﬀusion than Schemes 1 and 3, but with spurious oscillations, while Tables 1 and 2 indicate that Scheme 3 produces smaller L1 errors and

(a)

(b)

Fig. 5. Example 1: Comparison of numerical schemes for the simulated total solids volume fraction φ as a function of x: (a) t = 5, (b) t = 20.

October 14, 2008 11:9 WSPC/103-M3AS

1768

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

(a)

(b)

(c) Fig. 6. Example 1: Comparison of numerical schemes for the simulated phase density Φ as a function of ξ at x = 1 for three diﬀerent times. (a) Scheme 1, (b) Scheme 2, (c) Scheme 3.

better convergence rates. Moreover, Fig. 6 alerts to the increasing numerical viscosity of Scheme 1 as time advances, while solutions at x = 1 generated by Schemes 2 and 3 do not change. Therefore, we select Scheme 3 for the remaining examples.

5.1.2. Example 2: Steady state of a sedimentation of a suspension with equi-distribution of particle size We present a simulation up to steady state of the settling of a suspension with equi-distributed particle sizes. The initial condition is (5.1) with φtot (x) ≡ 0.1. To ˜ = 0.5 save computing time, we use the discretization parameters δx = 1/1800, λ and δξ = 1/32. Figure 7 shows the simulated phase density Φ as a function of ξ and x at a suﬃciently large time (t = 200), which represents the steady state. We observe that the steady state consists of a thick layer of large particles with a small amount of small particles above the bottom of the vessel, and a very thin and concentrated layer of small particles above the layer of large ones. Figures 8(a) and (b) show the simulated phase density Φ as a function of ξ at three diﬀerent positions, and the simulated phase density Φ as a function of x for three values of ξ, respectively.

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions Table 1. N

Example 1: Approximate L1 errors for the phase density Φ. t=1

e1 10−3 2 4 8 16 32

24.020 14.686 9.922 7.935 7.163

2 4 8 16 32

20.946 11.506 6.284 3.461 2.089

2 4 8 16 32

20.909 11.506 6.037 3.087 1.554

Table 2.

t = 20

Conv. rate

e2 10−6

Conv. rate

e2 10−6

Conv. rate

0.710 0.566 0.322 0.148

0.1082 0.4834 0.1678 0.2264 0.1965

Scheme 1 91.733 −2.160 89.446 1.526 88.188 −0.432 85.811 0.204 84.528

0.0364 0.0204 0.0394 0.0217

69.424 69.513 69.189 67.669 67.739

−0.0018 0.0067 0.0320 −0.0015

0.864 0.873 0.861 0.728

0.1893 0.3310 0.2335 0.1050 0.0795

Scheme 2 73.260 −0.807 49.957 0.503 31.725 1.153 26.695 0.402 24.806

0.552 0.655 0.249 1.059

189.835 91.324 65.466 39.736 21.648

1.056 0.480 0.720 0.876

0.862 0.930 0.968 0.990

0.1907 0.3310 0.0824 0.0358 0.0197

Scheme 3 73.140 −0.796 49.957 2.007 29.633 1.204 16.386 0.864 8.172

0.550 0.754 0.855 1.004

186.852 91.324 36.368 14.537 6.722

1.033 1.328 1.323 1.113

Conv. rate

e1 10−3

Example 1: Approximate L1 errors for the solids volume fraction φ.

N

t=1 e1 10−3

2 4 8 16 32

49.281 24.240 11.723 5.470 2.344

2 4 8 16 32

49.292 24.257 11.737 5.476 2.348

2 4 8 16 32

49.290 24.257 11.730 5.471 2.344

Conv. rate

e2 10−3

1.0236 1.0481 1.0996 1.2224

49.232 24.226 11.722 5.470 2.344

t = 20 Conv. rate

e2 10−3

Conv. rate

Scheme 1 49.265 1.0231 24.259 1.0473 11.754 1.0995 5.491 1.2224 2.492

1.0221 1.0453 1.0979 1.1397

49.231 24.225 11.722 5.471 2.345

1.0230 1.0472 1.0995 1.2223

1.0230 1.0474 1.0998 1.2219

49.232 24.225 11.722 5.470 2.344

Scheme 2 49.194 1.0231 24.281 1.0473 11.810 1.0995 5.559 1.2224 2.439

1.0187 1.0398 1.0872 1.1886

49.118 24.194 11.709 5.466 2.343

1.0216 1.0470 1.0992 1.2222

1.0229 1.0482 1.1002 1.2226

49.232 24.225 11.722 5.470 2.344

Scheme 3 49.189 1.0231 24.281 1.0473 11.757 1.0995 5.479 1.2224 2.347

1.0185 1.0463 1.1016 1.2228

49.119 24.194 11.713 5.467 2.343

1.0217 1.0465 1.0992 1.2223

Conv. rate

e1 10−3

1769

October 14, 2008 11:9 WSPC/103-M3AS

1770

Fig. 7.

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

Example 2: Sediment composition near steady state (t = 200): simulated phase density Φ.

(a)

(b)

(c) Fig. 8. Example 2: Sediment composition near steady state (t = 200): simulated phase density Φ as a function (a) of ξ at three diﬀerent positions x and (b) of x for three diﬀerent values of ξ; and (c) total solids volume fraction φ.

Figure 8(c) shows the simulated total solids volume fraction φ as a function of x at steady state (t = 200). In Fig. 8(c) we clearly observe that above around x = 0.875, we have φ = 0, while φ = φmax below this position (both equalities hold to within machine accuracy). Figure 8(c) conﬁrms that there is conservation of

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

Fig. 9.

1771

Example 2: Simulated total solids volume fraction φ.

mass or volume because the area between the graph and φ = 0 is approximately 0.8(1 − 0.875) = 0.1, which is the same value of the area under the graph at t = 0 (initial volume of the solids). Figure 9 shows the simulated total solids volume fraction φ until steady state (t = 200). 5.1.3. Example 3: Comparison between sedimentation of suspensions with equi-distribution and Rosin–Rammler distribution of particle size We now simulate the batch sedimentation of suspensions with the particle size distribution functions Feq and FRRn . The initial total solids volume fraction is φtot (x) ≡ 0.1; the Rosin–Rammler parameters are m = 2 and d1 /l = 2; and the ˜ = 0.5 and δξ = 1/32. discretization parameters are δx = 1/3600, λ Figure 10 shows the simulated phase density Φ for suspensions with equidistribution and Rosin–Rammler distribution of particle sizes at t = 1 and t = 20. In Fig. 10 we observe the eﬀect of a greater ratio of small particles to large particles of the Rosin–Rammler distribution compared with the equi-distribution. Figure 11 shows the corresponding simulated total solids volume fractions φ. In Fig. 11 we observe that the suspension with Rosin–Rammler distribution of particle size settles more slowly than the suspension with equi-distribution of particle size, due to the mentioned greater ratio small particles to large particles of the Rosin– Rammler distribution. Figure 12 shows the simulated total solids volume fraction φ of suspensions with equi-distribution and Rosin–Rammler distributions of particle sizes. 5.2. Examples 4 and 5: Sedimentation of a suspension with Riemann initial data In this pair of examples, we solve numerically two Riemann problems for the sedimentation of a suspension with equi-distributed particle sizes. The respective initial

October 14, 2008 11:9 WSPC/103-M3AS

1772

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

(a)

(b)

(c)

(d)

Fig. 10. Example 3: Simulated phase density Φ at t = 1 and t = 20: F (ξ) = Feq (ξ) for (a), (c), F (ξ) = FRRn (ξ) for (b), (d) .

(a)

(b)

Fig. 11. Example 3: Simulated total solids volume fraction φ of suspensions particle size distributions Feq and FRRn : (a) t = 1, (b) t = 20.

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

(a)

1773

(b)

Fig. 12. Example 3: Simulated total solids volume fraction φ with particle size distributions (a) Feq , (b) FRRn .

total solids volume fractions for Examples 4 and 5 are 0 for x < 0, (Example 4), φtot (x) = 0.1 for x ≥ 0 0.3 for x < 0, (Example 5). φtot (x) = 0 for x ≥ 0 ˜ = 0.5 and δξ = 1/32. We use the discretization parameters δx = 1/3600, λ Figure 13 shows the simulated phase density Φ as a function of the normalized squared size and the spatial position, at t = 2. We see, as in Figs. 3(b) and (c) of Example 1, that there is a considerable accumulation of medium size particles on the back of the wave of solid particles. The explanation of this phenomenon is given in Example 1. Moreover, it is clear that larger particles move faster than smaller ones. Figures 14(a) and (b) show the simulated phase density as a function of ξ at two spatial positions and three diﬀerent times. Figures 14(c) and (d) show the

Fig. 13.

Example 4: Simulated phase density Φ at t = 2.

October 14, 2008 11:9 WSPC/103-M3AS

1774

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

(a)

(b)

(c)

(d)

Fig. 14. Example 4: Simulated phase density Φ (a, b) as a function of ξ for three diﬀerent times: (a) x = 0.25, (b) x = 0.5; and (c,d) as a function of x for three diﬀerent times: (c) ξ = 0.015625, (d) ξ = 0.32812.

Fig. 15.

Example 4: Simulated total solids volume fraction φ at t = 1, t = 2 and t = 20.

simulated phase density as a function of x for ξ = 0.015625 and ξ = 0.32812 at three diﬀerent times. Figure 15 shows the simulated total solids volume fraction as a function of x at three diﬀerent times, and Fig. 16(a) shows the simulated total solids volume fraction until t = 20. The “steps” in the graph are a consequence

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

(a) Fig. 16.

1775

(b)

Examples 4 (a) and 5 (b): simulated total solids volume fraction φ.

of the discretization of the ξ-axis. Figure 16(b) shows the corresponding result for Example 5. In addition, we use this example, and Example 5, to analyze the behavior of the eigenvalues as the number of species N is increased. To this end, we recall ﬁrst that the exact solution of Example 4 depends on ω := x/t only, so it is sufﬁcient to calculate the solution for one ﬁxed time. We chose the cases N = 2, 4, 8 and 32 for close inspection, and plotted the eigenvalues λ1 ≥ λ2 ≥ · · · ≥ λN of Corollary 3.1, evaluated at the numerical solution, as a function of ω; according to Lemma 3.1, these eigenvalues coincide with those of the Jacobian Jf (see Theorem 2.1) if both are evaluated at corresponding vectors φ and Φ, respectively. For Example 4, Figs. 17 and 18 display the eigenvalues λ1 , . . . , λN (solid lines) and corresponding transformed eigenvalues ξλ1 , . . . , ξλN (dashed lines) for N = 2, N = 4 and N = 8 (Fig. 17) and N = 32 (Fig. 18). We recall that the transformed eigenvalues ξλ1 , . . . , ξλN are deﬁned by (3.42). The dashed lines are the bounds established by (2.9) and (2.10) for λ1 , . . . , λN and (3.12) and (3.13) for ξλ1 , . . . , ξλN , respectively. Figures 19 and 20 show the analogous results for Example 5.

6. Conclusions In this paper we propose the equation of the generalized kinetic theory (2.16), along with its “cumulative” version (2.23), for the sedimentation of a suspension with continuously distributed particle sizes. Our main emphasis has been on gaining insight into the solution by numerical experiments. We point out that unlike common kinetic formulations leading to an equation of Boltzmann type, our equation of the generalized kinetic theory (2.16) has a transport term depending nonlinearly ¯ and ξ¯ are involved. on Φ and nonlocally on ξ, since the quantities Φ In the present model the collision term that usually appears in kinetic models is zero. Note that since the kinetic parameter ξ is equivalent to squared particle size, not velocity, we cannot incorporate non-destructive particle–particle contacts by a kinetic collision term. On the other hand, although in the present work the

October 14, 2008 11:9 WSPC/103-M3AS

1776

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 17. Example 4: (a, c, e) eigenvalues λ1 , . . . , λN (solid lines) and (b, d, f) corresponding transformed eigenvalues for (a, b) N = 2, (c, d) N = 4 and (e, f) N = 8. The dashed lines are the bounds given by (a, c, e) (2.9), (2.10) and (b, d, f) (3.12), (3.13).

solid particles are considered to be rigid, one could additionally assume that they undergo breakage, or that the particles are actually small droplets in a liquid–liquid emulsion that aggregate (coagulate or coalesce). These eﬀects can be described by nonlocal terms in population balance models.43 Recent advances in their modeling and discretization include Refs. 20, 29, 33, 40–42. The next step should therefore to combine the present kinetic hindered settling model with the population balance

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

(a)

1777

(b)

Fig. 18. Example 4: (a) eigenvalues λ1 , . . . , λ32 (solid lines) and (b) corresponding transformed eigenvalues ξλ1 , . . . , ξλ32 (solid lines) for N = 32. The dashed lines are the bounds given by (a) (2.9), (2.10) and (b) (3.12), (3.13).

approach. On the other hand, the model could be adapted to particles made of diﬀerent materials, i.e. having diﬀerent densities. However, it is well known9, 21, 48 that such mixtures, even if the particles of each species have one size only, exhibit unstable sedimentation behavior, which is reﬂected by loss of hyperbolicity of the model equations, and the appearance of complex-conjugate pairs of eigenvalues. The analysis of the eigenvalue problem shows that also in the limit N → ∞, we can expect that the solution information travels at ﬁnite speeds. This information is useful when further developing numerical schemes for the equation of the generalized kinetic theory (2.16) (or (2.23)). However, our analysis is limited in that no closed-form expression for the eigenvectors or eigenfunctions is derived, at least not for the general case, in which Φ∗ and Ψ∗ are possibly discontinuous functions. Furthermore, as stated in Remark 3.1, the eigenvalues of A can possibly be located even more precisely, and we conjecture that the spectrum of A actually consists of ¯ plus one isolated point. ¯ ξ) I(Φ, The fact that the eigenfunction Ψ∗ is in general a discontinuous function, which will usually occur if the suspension under study is one with a ﬁnite number of size classes, combined with the observations that the operator A corresponds to a local linearization of the equation of the generalized kinetic theory (2.23) and that the eigenvalue problem appears in the solution of a Riemann problem, suggests that the solution of (2.16) or (2.23) will, in general, be measure-valued, as are the solutions for the kinetic limit of a multi-species traﬃc model considered in Ref. 13. Indeed, it seems interesting to investigate whether the main result of Ref. 13, namely, the equivalence between the well-posedness of the traﬃc model for all ﬁnite numbers of species and that of the corresponding kinetic model, could also be proved to hold for the model (2.4)–(2.7) for ﬁnite N and the kinetic model (2.16)–(2.19). However, for several reasons it is not entirely obvious how to achieve this, since the present model makes sense only if we impose boundary conditions (which do not appear for a traﬃc model); the ﬂuxes of the sedimentation model are algebraically more

October 14, 2008 11:9 WSPC/103-M3AS

1778

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 19. Example 5: (a, c, e) eigenvalues λ1 , . . . , λN (solid lines) and (b, d, f) corresponding transformed eigenvalues ξλ1 , . . . , ξλN (solid lines) for (a, b) N = 2, (c, d) N = 4 and (e, f) N = 8. The dashed lines are the bounds given by (a, c, e) (2.9), (2.10) and (b, d, f) (3.12), (3.13).

¯ x)) in the kinetic model (2.19) has no complicated (for example, the factor (ξ − ξ(t, counterpart in the traﬃc model); and as a consequence of the latter point, a pair of a convex entropy function and a convex entropy ﬂux for (2.4), (2.5) has not yet been found. Finally, let us emphasize that our numerical experiments further support the superiority of Scheme 3 for multi-species kinematic ﬂow models with velocities of

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

(a)

1779

(b)

Fig. 20. Example 5: (a) eigenvalues λ1 , . . . , λ32 (solid lines) and (b) corresponding transformed eigenvalues ξλ1 , . . . , ξλ32 (solid lines) for N = 32. The dashed lines are the bounds given by (a) (2.9), (2.10) and (b) (3.12), (3.13).

variable sign, which conﬁrms the observations made in Ref. 19 (our Scheme 3 is based on Scheme 8 of Ref. 19). We point out that the L1 convergence result of Lemma 3.2 for the ﬁnite-dimensional, linearized operators is consistent with the observed convergence of the L1 error e1 (with respect to δξ ) in Table 1 and of the errors e1 and e2 (for φ) in Table 2. Furthermore, we comment that Figs. 17–20 illustrate that Scheme 3 apparently approximates the correct solution of the Riemann problem satisfying the Liu entropy condition, since a fan of kinematic shocks, as is visible in Example 4 (Figs. 17 and 18), is characterized by stepwise decreases of eigenvalues, while successions of rarefaction fans, which emerge in Example 5 (Figs. 19 and 20), are characterized by successive increases of eigenvalues. See Refs. 14 and 23 for a precise formulation of this statement and further details. Appendix A. Proof of Lemma 3.2 Proof. (Proof of Lemma 3.2) We show that {AN (Φ∗,N )}N ∈N is a Cauchy sequence with respect to · 1 , which converges to a limit operator A. We need to consider that the diﬀerence between functions AN (Φ∗,N )Ψ and AM (Φ∗,M )Ψ is piecewise i−1 i constant, where the intervals of constancy are Ii := [ N M , N M ), i = 1, . . . , NM . To ∗ this end, we select Ψ ∈ F. Taking into account that 1 ¯ for all N ∈ N, Φ∗,N = Φ(ξ) dξ =: Φ (A.1) N 0

we obtain from (3.7) and (3.14) the following expression for i = 1, . . . , N M : (AN (Φ∗,N ) − AM (Φ∗,M ))Ψ∗ |Ii

N −1 M−1 ∗,M 1 1 ∗,N ∗,N ∗,M ¯− ¯− ¯ = −Ψ i Φ Φ Φ + Ψ i Φ V (Φ) M N N j=1 j M j=1 j

October 14, 2008 11:9 WSPC/103-M3AS

1780

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

+ −

N −1 M−1 1 ∗,N 1 ∗,M N,∗ ∗,M ∗,N ¯ − Ψk Ψk Φ i + Ψ M − Φ i V (Φ) M N N M

Ψ∗,N N

k=1

k=1

i i M −1 N −1 1 1 i i ∗,N ∗,N ∗,M ∗,M ¯ Ψk Ψk − V (Φ) + Ψ i − Ψ i − M N N M M N k=1

+

−Ψ∗,N N

Ψ∗,N i M

+ Ψ∗,M M

1 − M

Ψ∗,M i N

k=1

N −1 ¯− 1 Φ Φ∗,N N j=1 j

M−1 ∗,M ¯− 1 Φ Φ M j=1 j

1 − N

i M −1 i ∗,N ∗,N Φj Φ i − M M j=1

i N −1 i ∗,M ∗,M ¯ Φj V (Φ) Φ i − N N j=1

¯ + I 2 V (Φ) ¯ + I 3 V (Φ) ¯ + I 4 V (Φ). ¯ =: Ii1 V (Φ) i i i Our goal is now to estimate (AN (Φ∗,N ) − AM (Φ∗,M ))ΨL1 [0,1] =

NM 1 |(AN (Φ∗,N ) − AM (Φ∗,M ))Ψ|Ii |. N M i=1

(A.2) To this end, we note that

Ii1

=

−∆Ψ M,N,i

N −1 ¯− 1 Φ Φ∗,N N j=1 j

∆Φ , + Ψ∗,M i N

(A.3)

where we deﬁne ∗,N ∗,M ∆Ψ M,N,i := Ψ i − Ψ i , M

N

∆Φ :=

N −1 M−1 1 ∗,N 1 ∗,M Φj − Φ . N j=1 M j=1 j

Now, taking into account that ! i

i

" i i # i ∈ {1, . . . , N M } Ii ⊂ M , N ∪ N , M ≤ max{N, M }, N M M N we see that N M i=1

∗ |∆Ψ M,N,i | ≤ max{N, M }TV(Ψ ).

(A.4)

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1781

On the other hand, for

N M 1 ∗,N 1 ∗,M 1 ∗,N 1 ∗,M Φ Φ Φ − Φ − − ∆ = N j=1 j M j=1 j N N M M

k/N /M N M k/M /N 1 1 ¯ 1 − Φ = Φ(ξ) dξ − Φ(ξ) dξ − NM N M 0 0 k=1

NM 1 1 ¯ 1 k/M /N − Φ (A.5) = Φ(s) ds − NM N M k/N /M Φ

k=1

we may apply similar arguments to deduce that Φ |N − M | + max{N, M } ∆ ≤ Φ1 . NM

(A.6)

The term in parentheses in (A.3) may be estimated by 2Φ1 . Thus, using (A.3), (A.4) and (A.6), we get NM 3 max{N, M } + |N − M | 1 1 Φ1 Ψ∗ F . |I | ≤ N M i=1 i NM

(A.7)

Since the terms Ii2 have the same structure as Ii1 with the roles of Φ and Ψ reversed, we also get NM 3 max{N, M } + |N − M | ∗ 1 2 Ψ F (Φ∗ 1 + Φ1 ). |I | ≤ N M i=1 i NM

(A.8)

Ψ ˜Ψ +D Furthermore, we rewrite Ii3 as I3 = DM,N,i M,N,i , where

i i/M ∗,N Ψ Ψ i − N Ψ∗,M DM,N,i := i N M M

N i i M N i 1 ∗,N − Ψ i + = − Ψ∗,M (Ψ∗,N i i ). M N M N M M N Noting that i i |M − N | N M − , ≤2 M N NM

(A.9)

we obtain, by applying an argument similar to the one which leads to (A.4), that NM 2|N − M | + max{N, M } ∗ 1 Ψ Ψ F . DM,N,i ≤ N M i=1 NM

(A.10)

October 14, 2008 11:9 WSPC/103-M3AS

1782

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

Furthermore, we have that i

−1 1 N 1 Ψ ˜ DM,N,i := Ψ∗,M − k M N

i M −1

k=1

k=1

Ψ∗,N k

1 ∗,M 1 ∗,M 1 ∗,N 1 Ψ i − Ψ∗,N Ψk − Ψk − i M N M N N M k=1 k=1  i  i

N N M M ∗,N 1 ∗,M 1  1 ∗,N Ψ  Ψ i − Ψ i = ∆M,N,k − Ψk/M − . NM M N N M i i N

i M

=

k=1

k=N N

Applying calculus similar to that leading to (A.7) and (A.8), we get NM 1 ˜Ψ 3 max{N, M } + |M − N | ∗ Ψ 1 |DM,N,i | ≤ N M i=1 NM

(A.11)

NM 1 3 4 max{N, M } + 3|N − M | ∗ Ψ 1 . |Ii | ≤ N M i=1 NM

(A.12)

and ﬁnally

In light of (A.1), we obtain by similar arguments the estimate NM 4 max{N, M } + 3|M − N | 1 4 (1 + Φ1 )Ψ∗ 2F . |Ii | ≤ N M i=1 NM

(A.13)

Collecting (A.7), (A.12) and (A.13) and noting that for the purpose of estimate, we may replace |N − M | by max{N, M }, we see that there exists a constant C depending on Φ1 and Ψ∗ F such that ! " 1 1 ∗,N ∗,M (AN (Φ ) − AM (Φ , ))ΨL1 [0,1] ≤ C max (1 + V ∞ ), (A.14) N M whose right-hand side can be made arbitrarily small by choosing N and M sufﬁciently large. Consequently, for a ﬁxed function Φ ∈ L1 [0, 1], the sequence of operators {AN (Φ∗,N )}N ∈N is a Cauchy sequence, which converges to a bounded linear operator A : F → L1 [0, 1]. Acknowledgments RB acknowledges support by Fondecyt project 1050728 and Fondap in Applied Mathematics (project 15000001). Part of this work was performed while MK was visiting the University of Concepci´on. He acknowledges support by Fondecyt project 7050211. AG acknowledges support by Centro de Investigaci´ on Cient´ıfica y Tecnol´ ogica para la Miner´ıa (CICITEM), Antofagasta, Chile.

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1783

References 1. K. Arlotti, N. Bellomo and E. De Angelis, Generalized kinetic (Boltzmann) models: Mathematical structures and applications, Math. Mod. Meth. Appl. Sci. 12 (2002) 567–591. 2. L. G. Austin, C. H. Lee, F. Concha and P. T. Luckie, Hindered settling and classiﬁcation partition curves, Minerals & Metallurgical Process. 9 (1992) 161–168. 3. M. Bargiel, R. A. Ford and E. M. Tory, Simulation of sedimentation of polydisperse suspensions: A particle-based approach, AIChE J. 51 (2005) 2457–2468. 4. M. Bargiel and E. M. Tory, Simulation of sedimentation and ﬂuidization of polydisperse suspensions via a Markov model, Chem. Eng. Sci. 61 (2006) 5575–5589. 5. M. Bargiel and E. M. Tory, An extension of the particle-based approach to simulating the sedimentation of polydisperse suspensions, Int. J. Mineral Process. 79 (2006) 235–252. 6. M. Bargiel and E. M. Tory, A ﬁve-parameter Markov model for simulating the paths of sedimenting particles, Appl. Math. Model. 31 (2007) 2080–2094. 7. M. Bargiel and E. M. Tory, Stability of tridisperse suspensions, Comput. Visual. Sci. 10 (2007) 163–170. 8. M. Bargiel and E. M. Tory, A particle-based approach to sedimentation and ﬂuidization, in Leading-Edge Applied Mathematical Modeling Research, ed. M. P. Alvarez (Nova Science Publishers, 2008), pp. 19–65. 9. G. K. Batchelor and R. W. Janse van Rensburg, Structure formation in bidisperse sedimentation, J. Fluid Mech. 166 (1986) 379–407. 10. N. Bellomo, A. Bellouquid, J. Nieto and J. Soler, Multicellular biological growing systems: Hyperbolic limits towards macroscopic description, Math. Mod. Meth. Appl. Sci. 17 (2007) 1675–1692. 11. N. Bellomo, M. Delitalia and V. Coscia, On the mathematical theory of vehicular traﬃc ﬂow I. Fluid dynamic and kinetic modelling, Math. Mod. Meth. Appl. Sci. 12 (2002) 1801–1843. 12. S. Benzoni-Gavage and R. M. Colombo, An n-populations model for traﬃc ﬂow, Eur. J. Appl. Math. 14 (2003) 587–612. 13. S. Benzoni-Gavage, R. M. Colombo and P. Gwiazda, Measure valued solutions to conservation laws motivated by traﬃc modelling, Proc. Roy. Soc. A 462 (2006) 1791–1803. 14. S. Berres and R. B¨ urger, On Riemann problems and front tracking for a model of sedimentation of polydisperse suspensions, Z. Angew. Math. Mech. 87 (2007) 665–691. 15. S. Berres, R. B¨ urger and H. Frid, Neumann problems for quasi-linear parabolic systems modeling polydisperse suspensions, SIAM J. Math. Anal. 38 (2006) 557–573. 16. S. Berres, R. B¨ urger, K. H. Karlsen and E. M. Tory, Strongly degenerate parabolichyperbolic systems modeling polydisperse sedimentation with compression, SIAM J. Appl. Math. 64 (2003) 41–80. 17. R. B¨ urger, F. Concha, K.-K. Fjelde and K. H. Karlsen, Numerical simulation of the settling of polydisperse suspensions of spheres, Powder Technol. 113 (2000) 30–54. 18. R. B¨ urger, K.-K. Fjelde, K. H¨ oﬂer and K. H. Karlsen, Central diﬀerence solutions of the kinematic model of settling of polydisperse suspensions and three-dimensional particle scale simulations, J. Engrg. Math. 41 (2001) 167–187. 19. R. B¨ urger, A. Garc´ıa, K. H. Karlsen and J. D. Towers, A family of schemes for kinematic ﬂows with discontinuous ﬂux, J. Engrg. Math. 60 (2008) 387–425. 20. R. B¨ urger, K. H. Karlsen and J. D. Towers, Closed-form and ﬁnite diﬀerence solutions to a population balance model of grinding mills, J. Engrg. Math. 51 (2005) 165–195.

October 14, 2008 11:9 WSPC/103-M3AS

1784

00318

R. B¨ urger, A. Garc´ıa & M. Kunik

21. R. B¨ urger, K. H. Karlsen, W. L. Wendland and E. M. Tory, Model equations and instability regions for the sedimentation of polydisperse suspensions of spheres, Z. Angew. Math. Mech. 82 (2002) 699–722. 22. R. B¨ urger and A. Kozakevicius, Adaptive multiresolution WENO schemes for multispecies kinematic ﬂow models, J. Comput. Phys. 224 (2007) 1190–1222. 23. C. M. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, 2nd edn. (Springer-Verlag, 2005). 24. R. H. Davis and H. Gecol, Classiﬁcation of concentrated suspensions using inclined settlers, Int. J. Multiphase Flow 22 (1996) 563–574. 25. R. H. Davis and M. A. Hassen, Spreading of the interface at the top of a slightly polydisperse sedimenting suspension, J. Fluid Mech. 196 (1988) 107–134. 26. E. De Angelis and M. Delitalia, Modelling complex systems in applied sciences; methods and tools of the mathematical kinetic theory for active particles, Math. Comp. Model. 43 (2006) 1310–1328. 27. Eds. P. Degond, L. Pareschi and G. Russo, Modeling and Computational Methods for Kinetic Equations (Birkh¨ auser, 2004). 28. M. Delitala and A. Tosin, Mathematical modeling of vehicular traﬃc: A discrete kinetic theory approach, Math. Mod. Meth. Appl. Sci. 17 (2007) 901–932. 29. F. Filbet and P. Lauren¸cot, Numerical simulation of the Smoluchowski coagulation equation, SIAM J. Sci. Comput. 25 (2004) 2004–2028. 30. F. Filbet, P. Lauren¸cot and B. Perthame, Derivation of hyperbolic models for chemosensitive movement, J. Math. Biol. 50 (2005) 189–207. 31. H. P. Greenspan and M. Ungarish, On hindered settling of particles of diﬀerent sizes, Int. J. Multiphase Flow 8 (1982) 587–604. 32. Z. Ha and S. Liu, Settling velocities of polydisperse concentrated suspensions, Canad. J. Chem. Engrg. 80 (2002) 783–790. 33. J. Kumar, M. Peglow, G. Warnecke, S. Heinrich and L. M¨ orl, Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique, Chem. Engrg. Sci. 61 (2006) 3327–3342. 34. A. Kurganov and E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection-diﬀusion equations, J. Comput. Phys. 160 (2000) 241–282. 35. H.-S. Law, J. H. Masliyah, R. S. MacTaggart and K. Nandakumar, Gravity separation of bidisperse suspensions: Light and heavy particle species, Chem. Engrg. Sci. 42 (1987) 1527–1538. 36. M. J. Lockett and K. S. Bassoon, Sedimentation of binary particle mixtures, Powder Technol. 24 (1979) 1–7. 37. J. H. Masliyah, Hindered settling in a multiple-species particle system, Chem. Engrg. Sci. 34 (1979) 1166–1168. 38. A. M. Ostrowski, Solution of Equations in Euclidian and Banach Spaces, 3rd edn. (Academic Press, 1973). 39. Eds. R. H. Perry, D. W. Green and J. O. Maloney, Perry’s Chemical Engineer’s Handbook, 7th edn. (McGraw-Hill, 1997). 40. S. Qamar, M. P. Elsner, I. A. Angelov, G. Warnecke and A. Seidel-Morgenstern, A comparative study of high resolution schemes for solving population balances in crystallization, Comput. Chem. Engrg. 30 (2006) 1119–1131. 41. S. Qamar, A. Ashfaq, G. Warnecke, I. Angelov, M. P. Elsner and A. SeidelMorgenstern, Adaptive high-resolution schemes for multidimensional population balances in crystallization processes, Comput. Chem. Engrg. 31 (2007) 1296–1311.

October 14, 2008 11:9 WSPC/103-M3AS

00318

Sedimentation of Polydisperse Suspensions

1785

42. S. Qamar and G. Warnecke, Numerical solution of population balance equations for nucleation, growth and aggregation processes, Comput. Chem. Engrg. 31 (2007) 1576– 1589. 43. D. Ramkrishna, Population Balances (Academic Press, 2000). 44. J. F. Richardson and W. N. Zaki, Sedimentation and ﬂuidization I, Trans. Inst. Chem. Engrg. (London) 32 (1954) 35–53. 45. P. Rosin and E. Rammler, The laws governing the ﬁneness of powdered coal, J. Inst. Fuel 7 (1933) 29–36. 46. F. Rosso and G. Sona, Gravity-driven separation of oil-water dispersions, Adv. Math. Sci. Appl. 11 (2001) 127–151. 47. W. Schneider, G. Anestis and U. Schaﬂinger, Sediment composition due to settling of particles of diﬀerent sizes, Int. J. Multiphase Flow 11 (1985) 419–423. 48. R. H. Weiland, Y. P. Fessas and B. V. Ramarao, On instabilities arising during sedimentation of two-component mixtures of solids, J. Fluid Mech. 142 (1984) 383–389. 49. G. C. K. Wong and S. C. Wong, A multi-class traﬃc ﬂow model — an extension of LWR model with heterogeneous drivers, Transp. Res. A 36 (2002) 827–841. 50. B. Xue and Y. Sun, Modeling of sedimentation of polydisperse spherical beads with a broad size distribution, Chem. Engrg. Sci. 58 (2003) 1531–1543. 51. A. Zeidan, S. Rohani and A. Bassi, Dynamic and steady-state sedimentation of polydisperse suspension and prediction of outlets particle-size distribution, Chem. Engrg. Sci. 59 (2004) 2619–2632. 52. A. Zeidan, S. Rohani, A. Bassi and P. Whiting, Review and comparison of solids settling velocity models, Rev. Chem. Engrg. 19 (2003) 473–530. 53. M. Zhang, C. -W. Shu, G. C. K. Wong and S. C. Wong, A weighted essentially nonoscillatory numerical scheme for a multi-class Lighthill–Whitham–Richards traﬃc ﬂow model, J. Comput. Phys. 191 (2003) 639–659. 54. P. Zhang, R.-X. Liu, S. C. Wong and S.-Q. Dai, Hyperbolicity and kinematic waves of a class of multi-population partial diﬀerential equations, Eur. J. Appl. Math. 17 (2006) 171–200. 55. P. Zhang, S. C. Wong and C.-W. Shu, A weighted essentially non-oscillatory numerical scheme for a multi-class traﬃc ﬂow model on an inhomogeneous highway, J. Comput. Phys. 212 (2006) 739–756.

Recommend Documents

A Generalized Model of Commitment

nominal cross recurrence as a generalized lag ... - World Scientific

A Simple Kinetic Model Describes the Processivity of Myosin-V

A neural model of border-ownership from kinetic occlusion, Vision