Mathematical and Numerical Modeling of 1-D and 2-D Consolidation

(1)

and 2-D Consolidation

Katarina Gustavsson

Doctoral Dissertation Royal Institute of Technology

Department of Numerical Analysis and Computer Science

(2)

ges till offentlig granskning för avläggande av teknisk doktorsexamen torsdagen den 24 april 2003 kl 10.15 i Kollegiesalen, Administrationsbyggnaden, Kungl Tekniska Högskolan, Valhallavägen 79, Stockholm.

ISBN 91-7283-481-1 TRITA-NA-0304 ISSN 0348-2952

ISRN KTH/NA/R--03/04--SE

Katarina Gustavsson, April 2003c ,

(3)

A mathematical model for a consolidation process of a highly concentrated, flocculated suspension is developed.The suspension is treated as a mixture of a fluid and solid particles by an Eulerian two-phase fluid model.We characterize the suspension by constitutive relations correlating the stresses, interaction forces, and inter-particle forces to concentration and velocity gradients.This results in three empirically determined material functions: a hysteretic permeability, a non-Newtonian viscosity and a non-reversible particle interaction pressure.Pa- rameters in the models are fitted to experimental data.

A simulation program using finite difference methods both in time and space is applied to one and two dimensional test cases.Numerical experiments are performed to study the effect of different viscosity and permeability models.

The eﬀect of shear on consolidation rate is studied and it is signiﬁcant when the permeability hysteresis model is employed.

ISBN 91-7283-481-1 • TRITA-NA-0304 • ISSN 0348-2952 • ISRN KTH/NA/R--03/04--SE

iii

(4)

(5)

Acknowledgments

I thank my advisors Jesper Oppelstrup and Bj¨orn Engquist for their guidance throughout this project.I would like to express my deepest gratitude to Jesper for his engagement and support.His advises, ideas and comments has been crucial for the ﬁnal result.

I am grateful to my co-author Bj¨orn Sj¨ogreen for his contributions to this work and for showing a great interest in this project.

I would also like to thank Anna-Karin Tornberg and Daniel Appel¨o for read- ing parts of this thesis and making suggestions for improvement.

Furthermore, I would like to thank Gunilla Kreiss for advice and support and Ingrid Melinder for her friendly encouragement.

Many thanks to all my friends and colleagues at NADA.

Jon Eiken and Bent Madsen (Alfa Laval Separation AB) have supplied the experimental data and invaluable engineering, process, and modeling experience and is gratefully acknowledged.

Sist men inte minst vill jag tacka Anders, Henrik och Hanna.Utan er - ingenting...

Financial support from Vinnova (former NUTEK) and Alfa Laval Separation AB through the Parallel Scientiﬁc Computing Institute (PSCI) and NADA is gratefully acknowledged.

vii

(8)

(9)

Introduction

Mathematical modeling and numerical approximation of processes involving suspensions are the main topics of this work.

A suspension is a composite material consisting of two or more immiscible constituents or phases.It can for example be a fluid with small solid particles, fibers or fluid particles.Suspensions are very common in many industrial applications such as

• paper manufacturing

• medical applications

• food industry

• mineral industry

• treatment of industrial and municipal waste water.

Therefore, it has become an important area of research and much effort is devoted to develop physical and mathematical models that describe the behavior of a suspension in different processes.One such process is separation, where by natural gravity or external forces, the suspension is separated into its different components.

In this thesis, we focus on the dewatering of industrial or municipal waste water sludge, modeled as a suspension with agglomerates of small particles, ﬂocs.The dewatering reduces the water contents in the sludge, which is of great importance before disposal.

Two diﬀerent processes are studied, see Figure 1.1. In the one dimensional process, gravitational forces and buoyancy separate the particle phase from the ﬂuid phase.The other process is a two dimensional process where, in addition to gravity, shear is used to improve and speed up the separation.

The goal is to develop and study a mathematical and numerical model of such processes with the aim to improve the understanding and, eventually, the design of separation process machinery.

1

(10)

0 0.1 0

0.1

clear fluid

(a) Separation by gravity.

0 0.1

clear fluid

(b) Separation by gravity and shear.

Shear induced by movement of bottom wall.

Figure 1.1. The two diﬀerent separation processes. Here, contour plots of the volume fraction of solids and mass ﬂux arrows are displayed. Numerical results are presented in more details in sections 4.4 and 5.5.

Mathematical models of a fluid with dispersed particles can be formulated in different ways depending on e.g. the application and the concentration of the suspension.One possibility is a microscopic description of the evolution of all particles and all their interfaces.[13, 31, 48] present results from direct simulations of rigid particles suspended in a liquid.However this is only feasible for a limited number of particles.In many practical situations, rather than considering each individual particle, a macroscopic description of the suspension is more appropriate.Then, the two phases (fluid and particle) are modeled as interacting continua, i.e. both the fluid and the solid phases are considered as continuous phases.To keep track of the different phases, a scalar volume fraction field is introduced indicating the proportion of the total volume occupied by the particle phase.A mathematical model is obtained by formulating effective or mean equations for conservation of mass and momentum for each phase.This view of the process is considered in the present work.

There is one major drawback with the macroscopic description compared to a microscopic description; The averaging process involved in obtaining the mean formulation leads to unspeciﬁed terms.In order to characterize the suspension and close the formulation, constitutive relations are needed.These are often hard to determine both from modeling as well as from experimental point of view.Here the suspension is characterized by three empirical material functions:

(11)

viscosity, permeability and yield pressure.These are ﬁtted to experimental data obtained on waste water sludge.

1.1 One Dimensional Gravity and Buoyancy Sep- aration

The one dimensional separation process is a classical problem and has been the subject of many papers, [3, 8, 25] among others.In the classical formulation, the mathematical model reduces to one partial diﬀerential equation for the volume fraction, here referred to as the classical or inviscid model.

In our study, the approach is to work with a slightly diﬀerent model which includes a term related to the viscous forces on the suspension.This term is small and is usually neglected in the classical way of treating this problem.We choose to keep this term and solve two coupled equations, one equation for the volume fraction and one for the velocity ﬁeld of the particle phase.We refer to this model as the viscous model.

From a numerical point of view, the viscous model is advantageous since a less restrictive condition on the time step (CFL condition) is obtained than with the classical model.This is indicated by a linear stability analysis and exempliﬁed by numerical experiments.

Numerical simulations as well as experimental data show a steep gradient in the volume fraction at the clear fluid interface, see Figure 1.1. The equation for the volume fraction in the classical approach is hyperbolic in the clear fluid region and parabolic otherwise.Mathematically, the regions are separated by a discontinuity.The viscous model gives a volume fraction profile with very steep gradients at the interface but it is not clear whether this is a true shock or not.

Still, the fact that steep gradients are developed has to be considered in the choice of numerical method.A high resolution scheme based on second order central ﬁnite diﬀerences is developed for the viscous problem.It is constructed in such way that it also works in the inviscid limit, i.e., when the viscous terms are neglected, the classical model is obtained, and the solution develops shocks.

1.1.1 Numerical Results

We present computational results for one dimensional gravity separation.The convergence of the high resolution scheme is studied clearly showing second order accuracy and well resolved sharp gradients.

A numerical comparison between the viscous and the inviscid model shows good agreement when the viscous term is small.The less restrictive CFL condition for the viscous model translates into faster simulation.

Furthermore, our numerical results are validated by a comparison to experimental data for alumina suspensions under natural gravity.Experimentally determined concentration proﬁles are compared to computed proﬁles and the agreement is good.

(12)

1.2 Gravity and Shear Induced Separation

The two dimensional process is quite diﬀerent from the one described above.

The separation is irreversible and the micro structure in the suspension is very important for the dewatering.These features have to be incorporated in the mathematical model.A memory function is introduced to treat the irreversibility effects.This function is also used as an additional state variable which describes, in a simplified way, effects of structure in certain constitutive relations.

Speciﬁcally, a new heuristic model for how the permeability function depends on micro structure is developed and tested numerically.

The mathematical model describing the 2-D process consists of an elliptic system of equations for the fluid pressure and the two velocity components of the particle phase coupled to the volume fraction.The equation for the volume fraction is approximated numerically by a two dimensional extension of the high resolution scheme developed for the one dimensional problem.The elliptic system is discretized with second order finite differences.

1.2.1 Numerical Results

Numerical experiments are performed to investigate the properties of the mathematical model and to study its sensitivity to the material functions and their influence on the solution.As it turns out the final result is sensitive to the choice of these functions and to parameters such as wall velocity.Furthermore, the effects of shear on the process are significant when the permeability model with irreversibility is employed.Below we show the time evolution of a typical 2-D computation, Figure 1.2.

1.3 Outline

The outline of this thesis is as follows.An industrial background is given and the physical process is discussed in chapter 2.In section 2.5 we present other work related to the topic.In chapter 3, the mathematical model is introduced together with the constitutive relations and their ﬁtting to experimental data.

One dimensional gravity induced separation models are discussed in chapter 4.

In section 4.3.3 the high resolution method is developed and described, and section 4.4 gives numerical results and comparison to experimental data. In chapter 5, the two dimensional gravity and shear induced process is described with mathematical analysis of the model and the numerical treatment.In section 5.5 we present numerical experiments and parametric studies of the most interesting process parameters.

(13)

0 0.02 0.04 0.06 0.08 0.1 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

(a)t = 10s

0.05

0.1

0.15

0.2

0.25

0.3

0 0.02 0.04 0.06 0.08 0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

(b)t = 50s

0 0.02 0.04 0.06 0.08 0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

(c)t = 100s

0.05

0.1

0.15

0.2

0.25

0.3

0 0.02 0.04 0.06 0.08 0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

(d)t = 150s

0 0.02 0.04 0.06 0.08 0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

(e)t = 200s

0.05

0.1

0.15

0.2

0.25

0.3

0 0.02 0.04 0.06 0.08 0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

(f)t = 250s

Figure 1.2. Numerical 2-D simulation of separation by gravity and shear, volume fraction of the solids phase.

(14)

This thesis is partially based on the following papers:

1.K.Gustavsson, Simulation of Consolidation Processes by Eulerian Two- Fluid Models, Licentiate’s thesis, Department of Numerical Analysis and Computing Science, Royal Institute of Technology (KTH), Stockholm, Sweden (1999).ISBN 91-7170-419-1, TRITA-NA-9907.

2.K.Gustavsson and J.Oppelstrup, A Numerical Study of the Consolida- tion Process of Flocculated Suspensions Using a Two-Fluid Model.Pro- ceedings, The Third European Conference on Numerical Mathematics and Advanced Applications (ENUMATH), World Scientiﬁc (1999)

3. K.Gustavsson and J.Oppelstrup, Consolidation of Concentrated Suspen- sions - Numerical Simulations Using a Two-Phase Fluid Model, Computing and Visualization in Science 3 (2000) 39-45

4.K.Gustavsson and J.Oppelstrup and J.Eiken, Numerical 2D Models of Consolidation of Dense Flocculated Suspensions, Journal of Engineering Mathematics 41(2001) 189-201

5.K.Gustavsson and J.Oppelstrup and Jon Eiken, Consolidation of Concen- trated Suspensions - Shear and Irreversible Floc Structure Rearrangements, Computing and Visualization in Science 4 (2001) 61-66

6.K.Gustavsson and B.Sj¨ogreen, Numerical Study of a Viscous Consoli- dation Model, Proceedings of the Ninth International Conference on Hy- perbolic Problems Theory, Numerics, Applications, (HYP2002), Springer (2002)

(15)

Physical Model

Continuous separation of solids from liquids to achieve thickening of suspensions, clariﬁcation of liquids, etc., is often carried out by decanter centrifuge. One typical application is dewatering of sludge from municipal and industrial waste water.

The decanter is a high speed rotational device with an outer cylindrical bowl and a screw conveyor installed inside the bowl, see Figure 2.1. The suspension is fed into the decanter through the feed pipe A, to the extreme left.“Centrifugal forces” cause the solids in the suspension to sediment and migrate towards the bowl wall producing a “cake” with high volume fraction of solids.The diﬀeren- tial speed between the bowl and the screw conveyor transports the solids cake towards the solid outlet B, to the right.The clariﬁed liquid is discharged from the liquid outlet C to the left.A more detailed description of the dewatering process in a decanter can be found in [56].

The movement of the screw conveyor shears the suspension.Shear has a large inﬂuence on the separation process and makes the process more eﬃcient:

Figure 2.1. A decanter centrifuge. The suspension is fed into the decanter through the feed pipe A, right. Directly to the left of the feed pipe is the solid phase outlet B. The feed pipe empties into the bowl little bit further to the left.

The clariﬁed liquid outlet C is on the extreme left. (Picture courtesy of Alfa Laval.)

7

(16)

A much higher volume fraction of the solid phase (dryer sludge) is achieved in a decanter centrifuge than in an ordinary centrifuge.

The box model introduced in [28] is a simpliﬁed and conceptual two dimen- sional model for how these processes might work, see Figure 2.2. This model serves as a conﬁguration for both experimental and mathematical modeling.

The separation is driven by a gravitational force directed downwards, per- pendicular to the bottom.Shear forces are induced by moving the bottom wall to the right, see Figure 2.2. Due to gravity, the particles fall to the bottom of the box where an increase in concentration is obtained.The moving bottom wall transports the solid phase to the right and a maximum concentration is obtained in the lower right corner of the box.

An experimental device was constructed by Alfa Laval to verify the box model, [47].The device is sketched in Figure 2.3.It is a centrifuge with an outer cylinder that rotates with a diﬀerent speed than the inner cylinder.The space between the two cylinders, the bowl, is divided into a number of chambers to obtain the eﬀect of the walls of the screw conveyor.

In the experiments, the bowl was ﬁlled with waste water, sped up, and stopped after a while.The dry solids content of the cake was determined by removing the remaining water by drying.In Figure 2.4 the concentration of the solid phase is visualized after a short time.The maximum concentration is obtained in the right lower corner where the concentrated part has climbed up along the side wall.This ﬁgure can be compared to Figure 1.2 where a numerical computation of a 2-D consolidation process is displayed.

The box model is also used for the numerical 2-D experiments presented in chapter 5.

2.1 Consolidation of Flocculated Suspensions

The waste water sludge is a two constituent ﬂuid-particle mixture, also called suspension, with many small solid particles dispersed in water.The suspension is dense which means that the volume fraction of particles is large, ( 10%).The density of the particles is higher than the density of the ﬂuid: The particles will settle (fall) towards the bottom and form a layer of sediment, i.e. a separation process will take place.

When certain surface active substances, known as ﬂocculants, are added, the particles attach to each other and form larger agglomerates of particles.

The agglomerates are called ﬂocs and the suspension is said to be ﬂocculated.

Flocculating a suspension speeds up the separation process because large ﬂocs will settle faster than small particles.

Under the influence of external forces, e.g. gravitational or centrifugal forces, the flocs can form ever denser aggregates.For a solid phase volume fraction higher than the so-called gel forming fraction, the structure in the suspension changes from individual flocs to a continuous porous network filled with fluid.

This network can support normal stresses and resists compression until the ap-

(17)

Figure 2.2. The original sketch of the conceptual model, [28], of a dewatering process of a ﬂocculated suspension. The suspension with iso-lines of the concentration is pictured before and after being subjected to shear and gravity. The maximum volume fraction is found in the lower right corner. (Picture courtesy of Alfa Laval.)

(18)

ω Ω

Figure 2.3. Experimental setup. The outer cylinder rotated with a diﬀerent speed than the inner cylinder. The centrifuge is divided into a number of chambers.

Figure 2.4. Experimental result with waste water sludge. Dry solids content in one chamber in the centrifuge sketched in Figure 2.3. The result is visualized after a short spin with wall relative movement. (Picture courtesy of Alfa Laval.)

(19)

plied external force exceeds a yield pressure.Once this critical force is exceeded, there will be a compression of the network resulting in ﬂuid being released and an increase in concentration of solids.This is what we refer to as a consolidation process.A detailed discussion of ﬂocculated suspensions can be found in e.g.[6].

2.2 Characterization of the Suspension

We need constitutive relations to accomplish the characterization of the material as well as to form a closed set of equations.As mentioned earlier, experimental data is correlated to three material functions: permeability, viscosity, and yield pressure.

A detailed discussion of the constitutive relations and the material functions, ﬁtted to experimental data, is found in section 3.2.

Yield Pressure

The compressive response of the suspension is described by the yield pressure.

When stresses in the suspension exceed the volume fraction dependent yield pressure, ps,yield(φ), the suspension consolidates to a higher volume fraction with a higher yield pressure.

If the ﬂocculation of the suspension is strong enough, the compression of the network is irreversible and a plastic deformation takes place.This means that even if the external compression load subsequently decreases, the suspension will remain compressed.

Permeability

The interaction between the solid and the fluid phase is due to drag forces on the particles or flocs when transported through the fluid.This can be described by Darcy’s law and the concept of permeability from the theory of flow through porous media [57].

Viscosity

A Newtonian ﬂuid exhibits a linear relation between the shear force, τ and the velocity gradients (shear rate, ˙γ)

τ = η ˙γ

where the viscosity η is independent of shear rate.Flocculated suspensions are non-Newtonian and viscosity decreases with increased shear rate.This kind of behavior is generally called shear thinning.Flocculated suspensions usually have very high viscosity at low shear rates and are strongly shear thinning.Viscosity also increases with increased volume fraction.A treatment of the rheological behavior of suspensions can be found in e.g. [6] and [38].

(20)

The basis for this work has been that the material data needed should be furnished by techniques which can be applied routinely.Traditional pressure ﬁltration devices and piston experiments in 1-D settings are used for determin- ing permeability and force-deformation relationships.Various types of rotating and oscillatory viscometers can be used to determine the shear behavior.For an overview of a number of techniques developed for mineral applications, see [32].It proposes the yield pressure and the hindered settling function - the lat- ter related to the permeability - as important characteristics.The work gives ample illustration of the experimental diﬃculties and how some of these can be overcome.

Usually the experimental apparatus is limited to a certain range of the parameters.In particular, the industrial sludge studied in this work has very low permeability and standard ﬁlter press devices require too long times to be prac- tically useful.

A more fundamental problem is the determination of volume fraction solids.

This is easy for mineral suspensions where the density of the solids is accurately known, but we have not found any experimental technique which can measure this quantity for organic sludge.Solids weight fraction, however, is easily determined by drying, and it has been necessary to convert this to solids volume fraction by assumptions on intra-cellular water content.

In what follows we assume that the models determined by 1-D devices are also valid in the 2-D case.

For details about the parameter ﬁtting to experimental data, see [36].

2.3 Consolidation Under Gravity

When no shear forces are present the sedimentation is solely driven by gravitational forces.In this gravity consolidation process, no variation occurs in the horizontal direction, and the process can be described as a 1-D process in the vertical direction.

If the suspension is conﬁned to a container with an impermeable bottom, three distinct zones can be recognized during the consolidation process, see Fig- ure 2.5.

At the bottom there is the compression zone with consolidated material where the stress equals the yield pressure and a variable volume fraction of solids, φ.

On top of that is the overburden where the stress is smaller than the critical yield stress, and the concentration is constant.This as also called the hindered settling zone.Finally, on the top there is a zone with clear liquid separated from the overburden with a sharp interface.

The consolidation process will continue until the yield stress balances the gravity and buoyancy forces.

(21)

clear fluid

G

φ0 φ

y

φ > φ 0

φ₀ φ0

Figure 2.5. Consolidation driven solely by gravity. Distribution of the volume fraction of solids, φ at some time t. Initially, the concentration was uniform, φ = φ0. The compression zone (φ > φ0), the overburden where the volume fraction is constant (φ = φ0) and the clear ﬂuid zone are visualized.

2.4 Box Model - Consolidation by Shear and Gra- vitational Forces

The box model introduces a shearing velocity ﬁeld.This process is very diﬀerent compared to the one described above.There is still consolidation due to gravitational forces and the solid phase consolidates on the bottom of the box.But in addition, there are shear forces acting on the suspension by the moving bottom wall, see Figure 2.6.

In this case, two zones can be recognized: A consolidated part on the bottom of the container and clear fluid on top of that.The concentrated part of the suspension moves with the bottom wall.As the concentration increases, the viscosity of the suspension will increase and the suspension behaves more like a solid and starts to climb up along the wall on the right side.There, the diverging velocity field will deform and break the suspension apart into smaller compressed flocs separated by larger channels.The fluid can flow through these channels more easily and the suspension can be further compressed.Such a change in structure of the suspension is incorporated in the mathematical model by the permeability.These structure effects are discussed in section 3.2.

The break-up of the consolidated part of the suspension is a plausible con- tributor to the faster dewatering that has been observed in experiments where shear and gravitational forces are active.

(22)

>

φ ₀

clear fluid

G φ

φ0

τ τ

Figure 2.6. The Box Model: A consolidation process driven by gravity and shear. Closed container with a moving bottom wall to induce shear. Initially, the concentration in container is uniform,φ = φ0.

2.5 Related Work

Consolidation processes subjected to gravitational forces only are well under- stood and have been treated both from physical, mathematical and numerical point of view.

A vast number of references treat consolidation or sedimentation of concentrated suspensions.Here we will only point out a few that have been used in this work.For an extensive review of the contributions to research in sedimentation and thickening made during the 20th century, see B¨urger and Wedland, [16].

In 1952, Kynch, [43], published a paper on mathematical models for the sedimentation of a particle-ﬂuid suspension considered as a continuum.This paper is often considered as the origin of modern sedimentation theory.Kynch’s model makes a correct description of the behavior of a suspension of equally sized, small and rigid particles but it does not however work for ﬂocculated suspensions forming compressible sediments.

A more general model and mathematical theory for the sedimentation of ﬂocculated suspensions in one dimension is presented in e.g Buscall and White, [18], and in Auzerais et al., [3].In the latter, a numerical solution to the prob- lem is also presented.This is the classical or inviscid model, mentioned in the introduction chapter.

This model is also considered in Dorobantu, [25], where two different numerical approaches are used to produce a solution; a fixed grid two-phase method and an interface tracking method.The numerical solution is compared to experimental data published by Bergström, [5].

The thermodynamics of batch consolidation of suspensions is presented by B. Raniecki and J. Eiken in [55].

Bürger and coworkers published a number of works concerning consolidation of flocculated suspension, see e.g.[10, 11, 13, 14, 15, 17].In B¨urger et al., [17], a phenomenological theory of consolidation of flocculated suspension is developed and a general mathematical model is derived following classical theory of continuum mechanics.This model is further analyzed in [15].In one space dimension it reduces to the model describing 1-D gravity consolidation mentioned above.

(23)

Mathematical analysis of this model is reported in [10] and numerical treatment in [11].In [13, 14], the theory for sedimentation of suspensions of small spheres of equal size and density is generalized to polydisperse suspensions.

The multidimensional consolidation problem has been treated by Ystr¨om, [63].A mathematical model is developed, well posedness of that model is considered, and the problem is also solved numerically in two dimensions.

An application to paper pulp is presented in Zahrai, [64], where a two phase ﬂuid model for dewatering of an accelerating ﬁbrous suspension is presented.

The inﬂuence of diﬀerent stress models for the solids phase on the dewatering is investigated.

(24)

(25)

Mathematical Models

A ﬂuid-particle system can either be considered on a meso scale where every particle or ﬂoc is explicitly kept track of or on a macro scale where the system is considered to consist of two interacting continua, see Figure 3.1.

Whether to consider the ﬂuid-particle suspension on a meso or macro scale depends on the application and whether the suspension is dilute or dense.For a dense suspension where the particles are small and the number of particles is very large, the macro scale description is more suitable.

In this chapter, we discuss a mathematical macro scale model of the physical process described in chapter 2.Both the 1-D gravitational consolidation process and the process driven by gravity and shear can be described by this mathematical model.

The starting point is the two-phase model, discussed in e.g Ungarish, [62], Drew and Passman, [27] and Bustos et al., [19].This model has been used in diﬀerent applications such as dewatering of paper pulp, Zahrai, [64], and in ﬂuidization processes of turbulent gas-solid suspensions, Enwald et al., [29].

It consists of two coupled sets of conservation equations for mass and momentum and the ﬁnal set of conservation equations resembles a Navier-Stokes system.This is sometimes referred to as the Eulerian/Eulerian approach, see [29].

This is also the starting point for many of the classical 1-D consolidation studies, treated in e.g Auzerais et al., [3] , and B¨urger and Concha, [8].

The mathematical model is first formulated without any specific applications in mind.To make the model fit to a specific process or material, constitutive relations are needed.This is discussed in section 3.2.

3.1 Eulerian Two-Fluid model

We are concerned with the macroscopic, averaged behavior of the suspension and a continuum approach is used to formulate a mathematical model.The

17

(26)

V ol. frac.: Φ(rrrrr,t) Solid: uuuu(rurrrr,t) Flocs::::: rrrrri(t)

Fluid: vvvvv(rrrrr,t)

Particles Continua

Φ*(rrrrr,t)

Figure 3.1. Diﬀerent levels of description

two phases are treated as separate interpenetrating continuous media and the presence of more than one phase in the suspension is modeled by the concept of volume fractions for each phase.This model is referred to as a two-ﬂuid model.

The general idea is to formulate the integral balance equations for a ﬁxed control volume containing both phases and a moving interface.This yields local instantaneous equations for each phase and local instantaneous jump conditions for the interaction between the phases at the interface.In the Eulerian approach these equations are then averaged in a suitable way.Diﬀerent averaging procedures and closure laws have been employed, [27, 29].

The model that we present below is based on the model discussed in e.g. [27], also treated in [17], where a general phenomenological theory is developed for a consolidation process of a ﬂocculated suspension, which reduces to the classical one dimensional model discussed in chapter 4.

The following is assumed to be valid for the suspension:

• The particles are small in comparison to the dimensions of the conﬁguration containing the suspension;

• The particles have the same density;

• Both phases are incompressible;

• No mass transfer takes place between the phases;

• Surface tension between the phases is neglected.

(27)

The two-ﬂuid model conservation equations for mass and momentum look as follows,

∂φ

∂t +∇ · (φu) = 0 (3.1)

∂

∂t +∇ · (v) = 0 (3.2)

ρsφDu

Dt + φ∇pf− ∇ · Ts= ρsφg + m (3.3) ρfDv

Dt + ∇pf− ∇ · τf = ρfg− m (3.4) where the material derivatives are deﬁned as

Du

Dt = ut+ u· ∇u Dv

Dt = vt+ v· ∇v.

Subscript s denotes the solid phase and f denotes the fluid phase. φ∈ R is the volume fraction of solids and = 1− φ. u, v ∈ R² are the velocities of the solid and fluid phase respectively. T and τ ∈ R^2×2 denote stress tensors and m ∈ R² is a phase interaction term. pf ∈ R is the fluid pressure. ρs and ρf are densities of the two phases.

These equations were postulated in [2] and used in [3] to model consolidation of suspended particles under the inﬂuence of gravity.

Equations (3.1) - (3.4) must be supplemented by constitutive relations, boundary and initial conditions.Boundary and initial conditions are discussed in con- nection with the numerical treatment discussed in chapter 4 and in chapter 5.

The constitutive models are discussed below.

3.2 Constitutive Relations

Constitutive relations specify how the material behaves and how the two phases interact with each other.They are based on experience with the behavior of the material and conﬁrmed by rheological measurements.The modeling is strongly dependent on experimental data.The experiments on industrial waste water sludge, see also section 2.2, have provided material data for the constitutive relations.

Here we consider models for the stress tensors in the ﬂuid and solid phase, T and τ and the phase interaction term m.

3.2.1 Irreversibility and Micro Structure

Consolidation is irreversible, i.e. once the ﬂocs in the suspension are compressed they will remain compressed.

(28)

To make the model mimic such behavior we introduce a memory function φ^∗(x, t), see [63], which describes the maximal φ encountered since t = 0 by a material particle.It satisﬁes the equation

Dφ^∗ Dt =

Dφ

Dt when φ≥ φ^∗, ^Dφ_Dt > 0,

0 otherwise. (3.5)

Since

Dφ

Dt = φt+ u· ∇φ = −φ∇ · u

the sign of ∇ · u determines whether the suspension is under compression (φ is increasing) or dilation (φ is decreasing) along a particle path.

The constitutive models should take into account the micro structure since micro structural eﬀects are believed to be of great importance for the process.

Here micro structure is incorporated in the mathematical model in a simpli- ﬁed way by generalizing φ^∗ to be an additional state descriptor in addition to the volume fraction in the phase interaction model.This is described in details in section 3.2.3.

3.2.2 Stresses: Yield Pressure and Non-Newtonian Vis- cosity

The solid phase is assumed to support shear forces only when velocity gradients are present, i.e., in this respect to behave like a ﬂuid, with stress tensor

Ts=−psI + τs (3.6)

where ps is the solid phase pressure, I is the unit tensor and τs= ηs(∇u + ∇u^T) + λs(∇ · u)I.

We denote the viscosity of the solid phase by ηsand λs=−2/3ηssimilar to the Stokes’ hypothesis for a single phase ﬂuid, [4].

In the material studied the mixture viscosity of the suspension will be orders of magnitude larger than the viscosity of the ﬂuid which we assume to be water.

We will therefore neglect the fluid stress tensor, τf, and assume that the mixture viscosity equals the solid phase viscosity.This simplification will not be valid in the clear fluid limit but its influence in the compression zone, where the concentration φ is large, will most likely be of little importance.With these assumptions, the solid phase stress tensor is

τs= η

∇u + ∇u^T −2

3(∇ · u)I

. (3.7)

where η denotes the mixture viscosity.

(29)

Viscosity

We assume that the mixture viscosity is a function of the volume fraction and the shear rate, η = η(φ,| ˙γ|).

The rate of strain tensor is deﬁned as

γ(u) =˙ 1

2(∇u + ∇u^T) (3.8)

and the shear rate,| ˙γ| is a scalar formed from the invariants of the rate of strain tensor

| ˙γ| =

1

2|(tr ˙γ)²− tr( ˙γ · ˙γ)|

tr ˙γ =

n

( ˙γ)nn.

At low shear rates the ﬂocs are assumed not to break and the viscosity is essentially constant, η₀, with respect to shear rate.This is called the lower Newtonian level.At higher shear rates the ﬂocs deform and begin to break down.

The velocity gradients induce an orientation of the ﬂoc structure to form layers separated by clear ﬂuid.This will cause the apparent viscosity to decrease and the suspension will be shear thinning.For high enough shear rate the viscosity will again reach a region with constant viscosity but at a lower level, the upper Newtonian level, η_∞.

The upper and lower regions are not well determined by available experimental data.The measurements were conducted at shear rates where the viscosity was still dependent on the shear rate.Thus η₀and η_∞are treated as somewhat adjustable parameters for which order of magnitude estimates are available.The viscosity is also strongly increasing with the volume fraction of particles.

Experimental data for the viscosity as a function of shear rate show a power law behavior and have been ﬁt to a Carreau-Yasuda type of model, [7], combined with a volume fraction dependence into (see Figure 3.2)

η(φ,| ˙γ|) = Cφ^m

1 + | ˙γ|

0.05

n

C = 9.15× 10⁷P as m = 3.4

n =−0.9

(3.9)

Since n is close to−1 , the ﬂow behaves almost like a Bingham ﬂuid, which has a shear rate independent critical shear stress.

(30)

10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³ 10⁰

10¹ 10² 10³ 10⁴ 10⁵ 10⁶ 10⁷

η [Pas]

shear rate [1/s]

φ=0.5 φ=0.3 φ=0.2 φ=0.1

(a)η(φ, | ˙γ|) as a function of | ˙γ| for diﬀerent values of φ.

0 0.1 0.2 0.3 0.4 0.5 0.6

10⁻⁶ 10⁻⁴ 10⁻² 10⁰ 10² 10⁴ 10⁶ 10⁸

η [Pas]

φ [−]

0.1 1/s 1 1/s 10 1/s 100 1/s

(b)η(φ, | ˙γ|) as a function of the volume fraction, φ, for diﬀerent values of| ˙γ|

Figure 3.2. The viscosity is not constant but strongly shear thinning for con- centrated suspensions. This means a decrease in viscosity with increasing shear rate. The viscosity also increases withφ.

(31)

Yield Pressure

The interaction between the particles is described by the solids phase pressure ps, related to the yield pressure.

We saw earlier that the suspension supports normal stresses until the force exceeds a volume fraction dependent yield pressure.When the external compression load decreases, the suspension will remain compressed, and the interparticle forces vanish immediately in dilation.The irreversibility eﬀects are modeled us- ing the memory function φ^∗ introduced in section 3.2.1, equation (3.5). The solids pressure, psis assumed to depend on φ and φ^∗ as

ps(φ, φ^∗) =

0 if φ < φ^∗,

ps,yield(φ) if φ≥ φ^∗. (3.10) To ﬁt the yield pressure to experimental data we assume a power law behavior with concentration, see Figure 3.3,

ps,yield(φ) = Cφⁿ C = 5× 10⁶ P a n = 4.846

(3.11)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0 1 2 3 4 5 6x 10⁵

ps [Pa]

φ [−]

(a) Loading

0 0.05 0.1 0.15 0.2 0.25 0.3

0 100 200 300 400 500 600 700 800 900 1000

ps [Pa]

φ [−]

(b) Loading and unloading

Figure 3.3. Yield pressure as a function of φ ﬁtted to experimental data. In the right ﬁgure there is unloading atφ^∗= 0.15 (dotted line).

A typical particle phase pressure - concentration relation is shown in Figure 3.3. In initial compression the pressure follows the “virgin’ curve given by (3.11), and when the concentration decreases from the maximal φ^∗ (here = 0.15) the pressure vanishes; for a subsequent loading cycle the pressure remains zero until φ reaches φ^∗, and then follows the “virgin” curve again.

(32)

3.2.3 Inter-Phase Momentum Transfer: Permeability

The phase interaction term, m, takes into account the relative motion between the particles and the ﬂuid and can be modeled as follows

m = 1− φ

β(φ)(v− u) + ρfφC(φ)D

Dt(v− u), (3.12) This model, valid only for low relative velocities v− u, has also been used in [2]

and [3].

The first term represents a drag force with a drag coefficient β(φ), related to the permeability of the material.The second term is a virtual mass force proportional to the mass of fluid accelerated by a particle in motion. C(φ) is called the virtual mass coefficient.The virtual mass term is an inertia force and will not be further discussed here since inertial effects are neglected based on the order of magnitude estimates in section 3.3.

In [17] the authors split m into a dynamic and a hydrostatic part, m = md+ mb.The dynamic part, md is similar to equation (3.12). The hydrostatic part is chosen such that the momentum equation for the ﬂuid yields∇pf = ρfg in equilibrium.This is already included in our deﬁnition of the momentum equation, (3.4).

Permeability

Let k(φ) denote the intrinsic permeability with the dimension [m²] so √ k is a length scale representative of some pore diameter.Sludge experimental data is obtained for K(φ) = k(φ)/µ. K(φ) is called the Darcy function or simply permeability function and µ is the viscosity of water.The data were ﬁtted to a power law function of φ,

K(φ) = Cφⁿ

C = 4.6× 10⁻¹⁸ m²/P as n =−7.41

(3.13)

Note the extreme variation for low values of φ, Figure 3.4.

The model has to be modiﬁed in the limit φ → 0 to give ﬁnite velocities.

The modification should be done in such a way that the sharp limit between the clear fluid and the suspension observed in practice is kept.The modification is defined and discussed in Section 4.2.2.

In a two dimensional process, the velocity field can diverge, and the solids fraction tends to decrease, hence changing the permeability.How the permeability changes, depends on details of floc sizes and structure which cannot be captured only by the volume fraction φ.As an example, the two structures in Figure 3.5 have the same volume fraction but different permeabilities.

A heuristic model for how permeability depends on micro structure as de- scribed by φ and φ^∗may be devised as follows.It has one parameter α related to

(33)

0 0.1 0.2 0.3 0.4 0.5 0.6 10⁻²⁰

10⁻¹⁵ 10⁻¹⁰ 10⁻⁵ 10⁰ 10⁵

φ [−]

K(φ) [m2/Pas]

Figure 3.4. K(φ) =_µ^k as a function ofφ for the sludge experimental data

ﬂoc sizes of a network breaking up in expansion.We have no direct experimental evidence for the value of α.However, the model produces qualitatively correct eﬀects.

Much more detailed theories have been put forward for yield stress models for weakly aggregated dispersions as in [54], using relations for fractal properties of ﬂocs, see [41].Inspired by such theories, considerations of “conceptual models”

for force balances, see [20], have been shown to produce results in correspondence with experiments.Other studies, such as [21], report only limited success in relating the observed power laws to fractal exponents.

Figure 3.5. Same φ, diﬀerent permeabilities

(34)

h

D

D+h h*

D+h*

1

Figure 3.6. Conceptual model of structure

Consider the medium in Figure 3.6. The n ﬂocs per unit volume are imper- meable of size S:

φ∼ nS^e, 1≤ e ≤ 3.

The exponent e may be related to the fractal properties of the ﬂocs, see [41].In compression, all inter ﬂoc channels have the same width h^∗, and the permeability is the one measured in 1D experiments.In subsequent expansion, we assume that only a fraction α of the channels widen:

1 = (1− α)n(S + h^∗)^e+ αn(S + h)^e

This expresses the relation between particle sizes and interparticle distances.

The resulting structure has two scales, S and S/α¹³ where it is assumed that the larger aggregates are compact so the relevant fractal exponent is 3.

Assuming a f -power law dependence on channel widths, with p liquid pres- sure, and D the φ-independent factor in the Darcy coeﬃcient gives a volumetric ﬂow rate Q (m/s)

Q

D∇p = (1− α)nh^∗f+ αnh^f.

Again f may be related to fractal exponents of the channel shapes.It follows, that with δ = 1/e, the ﬂow rate Q^∗in initial compression when h = h^∗ is

Q^∗

S^f^−eD∇p = φ^∗(1−fδ)(1− φ^∗δ)^f,

for low concentrations a power law which for many suspensions can be well ﬁtted to data.The resulting expression for the ﬂow rate in expansion becomes

(35)

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

−15

−10

−5 0 5 10

φ [−]

log k(φ,φ*)

φ^*=φ φ^*=0.1 φ^*=0.2 φ^*=0.3 φ^*=0.4

Figure 3.7. Permeability as function of volume fraction solids. The dashed- dotted line is the experimentally ﬁtted power law function including the modiﬁ- cation for lowφ for the permeability as a function of volume fraction only.

Q

S^f^−eD∇p = φ

(1− α)(1/φ^∗δ− 1)^f+

+α^1−fδ((1/φ− (1 − α)/φ^∗)^δ− 1)^f .

(3.14)

Figure 3.7 shows φ, φ^∗ dependence with f = 5, α = 0.005 and the exponent e = 1. f and D were chosen to make the permeability function ﬁt the values obtained from measurements for φ > φ^∗.Note the 10log-scale on the ordinate.

3.3 Dimensional Analysis

A dimensional analysis is performed to give the orders of magnitude of the diﬀer- ent terms in equations (3.1)-(3.4). A dimensionless formulation of the equations provides indicative parameters, insight about the relative size of diﬀerent terms and their relative importance.

Non-dimensional quantities for a ﬂow in a H× H box under g gravity are introduced as

x^∗= x/H y^∗= y/H t^∗= t/T u^∗= u/U v^∗= v/U p^∗= p/P ρ^∗= ρ/∆ρ ∆ρ = ρs− ρf e^∗_g= g/G.

(36)

A characteristic time scale is given by T = H/U where U is the free settling velocity of an individual ﬂoc at concentration φ₀. U is related to the permeability of the material as U = ∆ρGK₀ where K₀= K(φ₀) and is used as a scale factor for K(φ).Since the consolidation process is dominated by gravity, the pressure scale is chosen to balance the gravitational forces, P = ∆ρGH.

Equation (3.3) together with (3.6) and (3.12) can now be written in dimensionless form (the superscript * is dropped for simplicity)

FrρsφDu

Dt + φ∇p − 1

Pe_s∇ · τs=− 1

Pe∇ps(φ, φ^∗) + ρsφeg+

β(φ)(v− u) + FrC0φC(φ)d

dt(v− u)

(3.15)

where C₀= C(φ₀).We assume that C₀ isO(1), [63].

In the same way we can write equation (3.4) in dimensionless form together with equation (3.12) and under the assumption that τf can be neglected,

FrρfDv

Dt + ∇pf = ρfeg−

β(φ)(v− u)

− FrC₀φC(φ)d

dt(v− u).

(3.16)

The dimensionless P´eclet numbers, Pe, Pe_s, and the Froude number, Fr, are deﬁned as

Pe = G∆ρH

Ps0 Pe_s= G∆ρH

τ⁰ Fr = U²

GH (3.17)

with τ₀ = uwallη₀/H.The stress scale factors are η₀ = η(φ₀, uwall/H) and Ps0= ps(φ₀, φ₀) where uwall/H is a measure of the shear rate.Pe and Pe_s are P´eclet numbers which relate compression and shear strength of the material to the gravitational force.They can be viewed as a relation between the parameters related to the process and the material properties, yield pressure and shear strength.The Froude number, Fr, is the ratio of the inertia force to the gravity force.

3.3.1 Indicative Froude and P´ eclet numbers

Since the suspension is dense, the permeability will become extremely small in areas where the concentration of the particle phase is high.In these regions the sedimentation velocity will be low and Fr will be small, typically of order 10⁻⁹, compared to 1/Pe and 1/Pe_swhich are of order 10⁻⁵ and 10⁻³.For full details about indicative Froude and Peclet numbers, see Table 3.1.

To compute indicative numbers of the dimensionless variables φ = φ₀= 0.1 is chosen as a reference state.This is the initial volume fraction of particles in the suspension. Using (3.11) and (3.13) we obtain

K₀= K(φ₀)≈ 10⁻¹⁰ m²/P as (3.18)

(37)

and

Ps0= ps(φ₀, φ₀)≈ 10²P a. (3.19) Note that K₀and Ps0 are related only to the properties of the material and not to the speciﬁc process the suspension is subjected to.

From equation (3.9) the reference value of the viscosity, η₀ is computed as η₀= η(φ₀, uwall/H).

This value will depend on the bottom wall speed and is given in Figure 3.2. By using the relation between the sedimentation velocity, U and the permeability:

U = ∆ρGK₀, the dimensionless numbers, deﬁned in equation (3.17), can be written in the following form

Pe = ∆ρH

Ps0 G Fr = (∆ρK₀)²

H G Pe_s=∆ρH² η₀

G uwall

Approximative values of Pe, Pe_s and Fr for diﬀerent values of G and uwall, are given in Table 3.1. In all computations ∆ρ = 10³ kg/m³ and H = 0.1 m.

The values of K₀ and Ps0 are given by (3.18) and (3.19).

uwall[m/s] G [m/s²] 1/Pe_s 1/Pe Fr

10³ 10⁻² 10⁻³ 10⁻¹⁰

0.01 10⁴ 10⁻³ 10⁻⁴ 10⁻⁹

10⁵ 10⁻⁴ 10⁻⁵ 10⁻⁸

10³ 10⁻² 10⁻³ 10⁻¹⁰

0.1 10⁴ 10⁻³ 10⁻⁴ 10⁻⁹

10⁵ 10⁻⁴ 10⁻⁵ 10⁻⁸

Table 3.1. Peclet and Froude numbers for diﬀerent values of the wall speed and the G-number.

As long as uwalland G are in the range given in Table 3.1 the inertial terms can be neglected since Fr is smaller than both 1/Pe and 1/Pe_s.However, if the G-number or uwall is increased further the terms will eventually be of the same order and the inertial eﬀects have to be considered in the equations.

In the application of sludge dewatering some realistic values on G and uwall

are: G≈ 15000 m/s^sand uwall≈ 0.02 m/s.Most of the numerical computations are performed with G = 10000 m/s² and uwallranging from 0 to 0.01 m/s.

(38)

3.4 Reduced Model

Assuming that all dimensionless variables are of order O(1), then since Fr is much smaller than both 1/Pe and 1/Pe_s, all terms of O(Fr) can be discarded from the equations (3.15) and (3.16) yielding

φ∇pf − 1

Pe_s∇ · τs=− 1

Pe∇ps(φ, φ^∗) + ρsφeg+

β(φ)(v− u) (3.20)

∇pf = ρfeg−

β(φ)(v− u) (3.21)

From equation (3.21) we obtain a relation between u and v as

v = u − β(φ)(∇pf − ρfeg). (3.22) If we choose β(φ) = K(φ)/(1− φ) and introduce the reduced pressure

p = pf − p₀− ρfeg· x (3.23) where p₀ is a constant, this equation is recognized as the Darcy equation, see [57],

(1− φ)(v − u) = −K(φ)∇p.

If equations (3.1) and (3.2) are added, we obtain one equation for the mixture from which the ﬂuid velocity ﬁeld can be eliminated by using the relation (3.22).

∇ · (u − β(φ)(1 − φ)(∇p) = 0. (3.24) The ﬂuid velocity can also be eliminated from equations (3.20) and (3.21) by adding them

∇pf− 1

Pe_s∇ · τs=− 1

Pe∇ps(φ, φ^∗) + φ∆ρeg+ ρfeg. (3.25) From equations (3.1), (3.24) and (3.25), a coupled hyperbolic-elliptic system for the unknowns φ, p and u is obtained.With the constitutive relation for the stress tensor, τs, given by equation (3.7), it can be written as (here in dimensional form)

φt+∇ · (φu) = 0 (3.26)

∇ · (u − K(φ)∇p) = 0 (3.27)

∇p − ∇ ·

η(φ,| ˙γ|)(∇u + ∇u^T −2

3(∇ · u)I)

=−∇ps(φ, φ^∗) + φ∆ρGeg.

(3.28)

The permeability K is given by equation (3.13), the solid phase pressure, ps

is given by (3.10) and (3.11) and the viscosity by equation (3.9).

Mathematical and Numerical Modeling of 1-D and 2-D Consolidation

and 2-D Consolidation

Katarina Gustavsson

Doctoral Dissertation Royal Institute of Technology

Department of Numerical Analysis and Computer Science

Contents

Acknowledgments

Introduction

1.1 One Dimensional Gravity and Buoyancy Sep- aration

1.1.1 Numerical Results

1.2 Gravity and Shear Induced Separation

1.2.1 Numerical Results

1.3 Outline

Physical Model

2.1 Consolidation of Flocculated Suspensions

2.2 Characterization of the Suspension

2.3 Consolidation Under Gravity

2.4 Box Model - Consolidation by Shear and Gra- vitational Forces

2.5 Related Work

Mathematical Models

3.1 Eulerian Two-Fluid model

3.2 Constitutive Relations

3.2.1 Irreversibility and Micro Structure

3.2.2 Stresses: Yield Pressure and Non-Newtonian Vis- cosity

3.2.3 Inter-Phase Momentum Transfer: Permeability

3.3 Dimensional Analysis

3.3.1 Indicative Froude and P´ eclet numbers

3.4 Reduced Model