Keywords: Slurry -Rheology -Visco-plastic flows -Finite Volume Method -Computational Fluid Dynamics, iii Suspensions -Rhéologie -Écoulements visco-plastiques -Méthode des volumes finis -Dynamique numérique des fluides, 36

Phosphate ore slurry is a suspension of insoluble particles of phosphate rock, the primary raw material for fertilizer and phosphoric acid, in a continuous phase of water. This suspension has a non-Newtonian flow behavior and exhibits yield stress as the shear rate tends toward zero, which can be described as viscoplastic behavior. First, we proposed an innovative methodology to handle the numerical simulation of viscoplastic flows, based on a multigrid initialization algorithm in conjunction with the SIMPLE procedure. Using the proposed solver in combination with the regularization scheme of Papanastasiou, we chose the square lid-driven cavity flow and pipe flow as test cases for validation and discussion. In doing so, we study the influence of the Bingham number and the Reynolds number on the development of rigid areas and the features of the vortices within the flow domain. Regarding viscoplastic flows, our experiments demonstrate that our approach based on using the multigrid method as an initialization procedure makes a significant contribution by outperforming the classic single grid method. Afterwards, various grades and phosphate ore concentrations were chosen for this rheological investigation. We created some experimental protocols to determine the main characteristics of these complex fluids and established relevant rheological models with a view to simulate the numerical flow in a cylindrical pipeline. Rheograms of these slurries were obtained using a rotational rheometer and were accurately modeled with commonly used yield-pseudoplastic models. The results show that the concentration of solids in a solid-liquid mixture could be increased while maintaining a desired apparent viscosity. Then, a new rheological model for the flow of phosphate-water suspensions is proposed and it is found to be in good agreement with experimental data. The proposed model is then incorporated into the OpenFoam numerical code. The results demonstrate that the model is capable of reproducing the rheological behavior of phosphate suspensions at both low and high concentrations by comparing it with suitable models for modeling the rheological behavior of phosphate suspensions. The proposed model can provide a pathway for accurately simulating and monitoring phosphate slurry flows. The overall contributions of this work address a significant gap in the current understanding and simulation of phosphate ore slurry flow.

First and foremost, I would like to express my deepest appreciation to my esteemed supervisors, Pr. Fayssal Benkhaldoun, Pr. Driss Ouazar, and Pr. Abderrahim Maazouz. Their unwavering support, invaluable guidance, and innovative insights have been instrumental throughout my research journey. I am particularly indebted to them for introducing me to the intricate world of mathematical modeling and simulation. They have been much more than just mentors in research; their wisdom and kindness have enriched various facets of my life, and words fall short to express my profound gratitude.

Further, I wish to convey my gratitude to Pr Jean-Michel Ghidaglia, lead of the Modeling Simulation and Data Analysis (MSDA) program. His support and positivity have been both encouraging and uplifting. Similarly, I am indebted to Mr. Nabil Atlassi, the head of Flux, Planification et Performance at OCP Khouribga site, and his dedicated team. Favour of them, I have earned many soft skills., which will undoubtedly benefit me in future endeavors.

On a personal note, my journey would not have been the same without the support of my colleagues and friends. I am grateful to all those who shared in my daily life and activities, especially Issam, Ilias, Ayoub, Nidal, Ismail, Hamza, Raqi, Imad, Tarek, Moussa, Badr, Asmaa, and many others. Their encouragement, insights, and shared moments have made this voyage all the more memorable. I would also like to extend my gratitude to the esteemed jury members for their rigorous review and insightful feedback on my work. Your expertise and constructive comments have been invaluable in refining this thesis.

To each and every individual mentioned above, your support and belief in me have been the wind beneath my wings, and I remain forever grateful.

List of Figures

Introduction

Phosphate slurry pipelines have been shown to be an appropriate mode of transport for phosphate ore in terms of energy use and the generated emissions [START_REF] Rusconi | The Design and Engineering of the 187 Km Khouribga to Jorf Lasfar Phosphate Slurry Pipeline[END_REF]. Since phosphate rock is known to be insoluble, it cannot be used directly as a fertilizer. Consequently, it is normally processed in a wet-process phosphoric acid (WPA) plant [START_REF] Hore-Lacy | Production of Byproduct Uranium and Uranium from Unconventional Resources[END_REF]. Transporting wet phosphate from the washing plants over the pipeline is a more rational way to use water and energy resources, as no drying at the mine is required and the slurry can simply be driven down the pipeline by gravity.

Phosphate ore is found in sub-horizontal beds as sedimentary deposits with high phosphorus (P) concentrations of up to 30% alongside limestone, marls, and clays. In addition to potassium and nitrogen, phosphorus plays a crucial role in biological processes and is an essential constituent of all living matter, in both its animal and plant forms. Phosphate minerals are, therefore, raw materials for agricultural fertilizers [START_REF] Ilyin | Rare-Earth Geochemistry of 'old' Phosphorites and Probability of Syngenetic Precipitation and Accumulation of Phosphate1In Memory of Richard P. Sheldon1[END_REF]. Mining is accomplished by removing the overburden, followed by the extraction of phosphate ore for processing [START_REF] Dar | Sedimentary: Phosphates[END_REF]. Right after excavation, the phosphate is separated from the unwanted sand and clay minerals in a washing plant [5]. Then, the ore is mixed with water, turned into pulp, and the resulting slurry is transported through pipes to a chemical processing plant where it is mixed with sulfuric acid to produce phosphoric acid [START_REF] Becker | Phosphates and Phosphoric Acid: Raw Materials[END_REF].

The motivation for this rheological study was triggered by evidence that indicates that the transport characteristics and processing of phosphate ore slurry is heavily dependent on its rheological properties. Deep knowledge of the rheological behavior is therefore essential, particularly when transporting a large amount of slurry. It has been shown that when designing a pipeline, prediction of the flow regime (laminar or turbulent) and the head loss cannot be done properly without a rigorous understanding of the rheology [START_REF] Eshtiaghi | The Laminar/Turbulent Transition in a Sludge Pipeline[END_REF].

The rheology of a material is usually described using a constitutive equation that relates the applied shear stress τ to the resulting shear rate γ. The main objective of such a study is to predict the flow that would appear in a given piece of equipment under the action of applied forces. A highly viscous liquid requires more power to be pumped than one with a low viscosity. Knowledge of a material's rheological behavior is therefore useful when controlling pumping and piping systems [START_REF] Coussot | Yield Stress Fluid Flows: A Review of Experimental Data[END_REF][START_REF] Hart | A Review of Technologies for Transporting Heavy Crude Oil and Bitumen via Pipelines[END_REF].

Obtaining suspensions of particles in a fluid matrix is challenging and complex, and numerous investigations of particle rheology can be found in the literature. As reported in prior studies, concentrated suspensions may exhibit shear-thickening or shear-thinning behaviors (with or without a yield point), irreversibility under oscillating shear, and many other complex performances [START_REF] Stickel | Fluid Mechanics And Rheology Of Dense Suspensions[END_REF][START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF]. At high solids loadings (>25% by volume), observations of the rapid growth in the apparent viscosity, η = τ / . γ , and an increase in the non-Newtonian behavior (particularly shear-thinning), with Newtonian limiting behavior at both low and high shear rates, are reported in several articles [START_REF] Stickel | Fluid Mechanics And Rheology Of Dense Suspensions[END_REF][START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF][START_REF] Kelessidis | Optimal Determination of Rheological Parameters for Herschel-Bulkley Drilling Fluids and Impact on Pressure Drop, Velocity Profiles and Penetration Rates during Drilling[END_REF][START_REF] Eshtiaghi | Rheological Characterisation of Municipal Sludge: A Review[END_REF][START_REF] Núñez-Santiago | Rheological Evaluation of Non-Newtonian Mexican Nixtamalised Maize and Dry Processed Masa Flours[END_REF]. In piping systems, a shear-thinning slurry is often desired since the viscosity of the mixture can be reduced by increasing the shear rate, making the pumping more efficient. In most suspensions, the important physical properties that affect the rheology include a solid volume fraction, the particle size and shape, and the mechanical and physical interactions between particles and the fluid matrix.

Understanding the flow behavior of materials with yield remains a challenging problem at both experimental and computational levels. Nevertheless, numerical resolution of visco-plastic flows has significantly improved over the last few decades. Indeed, Computational Fluid Dynamics has firmly established itself across a wide range of applications and research areas. Subsequently, it has spread throughout engineering industry. This has gone hand in hand with major improvements in high-performance computing hardware and treatment of more sophisticated mathematical models. Notably, numerical simulations of visco-plastic flows have assisted in better understanding the behavior and fundamental properties of fluids exhibiting a yield stress in different applications relevant of natural and engineering sciences [START_REF] Apostolidis | Modeling of the blood rheology in steady-state shear flows[END_REF][START_REF] Bittleston | Mud removal: Research improves traditional cementing guidelines[END_REF][START_REF] Frigaard | Bingham's model in the oil and gas industry[END_REF].

Visco-plastic fluid flows account for an important branch of non-Newtonian fluid mechanics, as a considerable proportion of fluids flows occurring both naturally and industrially are known to exhibit yield stress. Suspensions of particles in a fluid matrix, such as slurries, gels, and nanocomposites, typically display viscoplasticy and the study of their behavior has led to a large literature. These fluids behave as solids as the shear rate tends to zero, and as liquids beyond a certain critical shear stress level, i.e. the yield stress. The fluid flow field is therefore divided into unyielded and yielded fluid regions. A comprehensive review of Visco-plastic fluids has been carried out by Barnes [18]. The most widely known model for inelastic visco-plastic fluids was first proposed by Bingham [START_REF] Bingham | An investigation of the laws of plastic flow[END_REF], which describes a linear interrelation between shear stress and shear rate during flow. The constant of proportionality of this linear relationship is often called the plastic viscosity of the yielded fluid, µ. Thus, the Bingham model can be expressed, in tensorial form, as follows:.

   γ = 0 if τ ≤ τ y τ = µ + τ y γ γ if τ > τ y (1.1)
where τ is the stress tensor and γ is the shear rate tensor, γ = ∇u + (∇u) T , ∇u being the velocity gradient tensor. The magnitudes τ and γ are evaluated based on the second invariants of the stress and the rate of strain tensors respectively, as :

τ = 1 2 τ : τ 1/2 (1.2) γ = 1 2 γ : γ 1/2 (1.3)
Later, more general forms of Bingham model were adopted in an attempt to capture the non-linear flow behavior beyond the critical shear stress. These slighlty more complex empirical models include, among others, Herschel and Bulkley, Casson and Robertson-Stiff models. The resulting constitutive equations combine the behavior of solids in the so-called 'rigid' regions and of non-Newtonian liquids in the 'flow' regions through discontinuous equations. Thus, simulating visco-plastic flows is likely to encounter many difficulties due to the non-differentiability of the constitutive equations and the indeterminacy of the stress tensor below the yield stress. In most cases, it is necessary to determine the location and the shape of the interface between the yielded and unyielded regions at which the flow must switches from one branch of the constitutive law to the other. In Eq. 5.1, the yield surface is represented by the locus of points where τ = τ y .

Since performing numerical simulations of a yield stress material flow is not a straightforward task, two main families of solution approaches were suggested over the past decades, the regularization method and augmented Lagrangian algorithm. The former approach includes methods which approximate Eq. 5.1 by one regularised and smooth constitutive equation which is well determined regardless of the shear rate magnitude. The regularised equation treats the whole material domain as a fluid of variable vis-cosity and locally assigns a large but finite value of viscosity to the unyielded regions. The most popular regularisation methods in the literature are proposed by Bercovier and Engelman [START_REF] Bercovier | A finite-element method for incompressible non-Newtonian flows[END_REF], and that by Papanastasiou [START_REF] Papanastasiou | Flows of materials with yield[END_REF]. Many authors have extensively studied viscosity regularisation methods for a range of different problems including the following list [START_REF] Tang | Bingham fluid simulation with the incompressible lattice Boltzmann model[END_REF][START_REF] Syrakos | Performance of the finite volume method in solving regularised Bingham flows: Inertia effects in the lid-driven cavity flow[END_REF][START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF][START_REF] Mitsoulis | Flow of Bingham plastics in a liddriven square cavity[END_REF][START_REF] Sverdrup | Highly parallelisable simulations of time-dependent viscoplastic fluid flow with structured adaptive mesh refinement[END_REF]. The other solution method in solving visco-plastic flows consists of using variational inequalities whose solutions are equivalent to the solutions of original problems. The underlying problem is solved as the minimization of shear rate or maximization of shear stress, which form the basis of the Augmented Lagrangian Method. The method also involves using Lagrange multipliers and the corresponding solution algorithm is of the Uzawa type. The development of this approach is due to Duvaut and Lions [START_REF] Duvaut | Rigid visco-plastic Bingham fluid[END_REF], together with Glowinski [START_REF] Glowinski | SUR L'ECOULEMENT D'UN FLU-IDE DE BINGHAM DANS UNE CONDUITE CYLINDRIQUE[END_REF]. Numerically solving visco-plastic flows using regularization methods is typically faster and easier to implement than the augmented Lagrangian algorithm, and therefore it has been common to use the first strategy [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF].

Several numerical methods were proposed for flow analysis of such materials in various bounded domains. The governing equation set of incompressible flows provides, along with momentum equations, a coupled continuity equation that should be satisfied. While every velocity component appears in each momentum equation and in the continuity equation, there is evidently no equation for the pressure which makes the formulation of an equation for pressure based on the governing flow equations non-trivial. The pressure-velocity coupling gives rise to a constraint in the solution: the applied pressure field in the momentum equations should verify continuity. Both issues associated with the pressure-velocity linkage and the further appearance of non-linearities in the set of equations have lead to a variety of iterative solution algorithms [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][START_REF] Ferziger | Computational methods for fluid dynamics[END_REF].

Historically, pressure and velocity were coupled using a staggered variable arrangement, in early works by Harlow and Welch [START_REF] Harlow | Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface[END_REF] and Chorin [START_REF] Chorin | Numerical Solution of the Navier-Stokes Equations[END_REF], with velocity typically defined at the cell faces, while other scalar variables including pressure were evaluated at the cell centres. However, iterative solution strategies were widely adopted for incompressible flows in the 1970s, notably through the development of the well-known segregated SIMPLE algorithm (semi-implicit method for pressure linked equations) [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF]. Since its formulation, the SIMPLE algorithm has provided the basis for the development of a number of SIMPLE-like algorithms devoted to algorithm improvements, a review of which is reported in [START_REF] Acharya | Pressure-based finite-volume methods in computational fluid dynamics[END_REF] whereas detailed descriptions are available in [START_REF] Versteeg | An introduction to computational fluid dynamics: the finite volume method[END_REF]. besides, the work of Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] expanded the application area of the SIMPLE-like algorithms by enabling the use of a collocated variable arrangement and setting the ground for a geometric flexibility similar to that of the finite element method (FEM) [START_REF] Schneider | Finite element analysis of incompressible fluid flow incorporating equal order pressure and velocity interpolation[END_REF]. The main issue arising from solving Navier-Stokes equations in primitive variables on collocated grid is the well known checkerboard pressure effect. The oscillation is caused by decoupling of pressure and velocity in discretised equations [START_REF] Miller | Use of a pressure-weighted interpolation method for the solution of the incompressible Navier-Stokes equations on a nonstaggered grid system[END_REF], and it is necessary to use special face velocity interpolation techniques for the calculation of mass flux through the control volume faces. In the early 1980s, Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] adopted a modified momentum interpolation to eliminate the irrational pressure field in co-located grid. Since then, the momentum interpolation of Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] has been widely used to approximate face velocity for SIMPLE based algorithms. Although SIMPLE is an old algorithm, it is still very popular and it has been used successfully in numerous studies to solve complex flow phenomena, [START_REF] Sverdrup | Highly parallelisable simulations of time-dependent viscoplastic fluid flow with structured adaptive mesh refinement[END_REF] being a recent example involving non-Newtonian flow. However, it has been reported that some limitations of this method may arise when dealing with visco-plastic flows. Reported results of simple iterative techniques to solve flows of fluids with yield show low convergence rates and significant computational costs for the solutions in some cases. SIMPLE is a slowly converging algorithm, but its performance can be greatly enhanced by using it in a multigrid (MG) context [START_REF] Syrakos | Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method[END_REF][START_REF] Roy | A parallel multigrid finitevolume solver on a collocated grid for incompressible Navier-Stokes equations[END_REF].

Multigrid techniques are implemented in conjunction with a linear solver using a sequence of coarsening and refining grids. The method is motivated by its capability to offer fast convergence by allowing all wavelengths of the algebraic error to decay uniformly fast [START_REF] Tong | A review of current progress in multiscale simulations for fluid flow and heat transfer problems: The frameworks, coupling techniques and future perspectives[END_REF]. The multigrid method was first adopted in the context of the finite volume methods by Sivaloganathan and Shaw [START_REF] Sivaloganathan | A multigrid method for recirculating flows[END_REF] and Peric et al. [START_REF] Peric | A finite volume multigrid method for calculating turbulent flows[END_REF] to solve recirculating Newtonian flows on staggered and co-located grids, respectively. These studies have shown that the efficiency of calculating incompressible flows using implicit methods based on SIMPLE algorithm can be substantially improved using the multigrid method for outer iterations. However, implementing a multigrid algorithm is a non-trivial task that requires considerable thought and care. One of the challenges involves inter-grid transfers, which may be dependent on a dynamic criterion based on the residuals. However, an initial solution can be obtained with relative ease and lower computational cost. While some methods are well defined and classified in a mathematical sense, they can still be a source of new insight and of unexpected new practical applications.

This thesis expects to contribute in a number of ways to the literature of slurry-related rheology and numerical simulations and especially their understanding. The major contributions are:

1. First contribution relates to an innovative methodology for handling the numerical simulation of viscoplastic flows using a multigrid initialization algorithm in conjunction with the SIMPLE procedure.

2. Second contribution concerns a study of rheological properties of phosphate suspensions. Based on experimental data collected from a series of rheometric tests, we investigated the effect of particles concentration and size, and shear rate on the viscosity and flow behavior of these suspensions.

3. Third contribution relates to a new rheological model that accurately reproduces the rheological behavior of phosphate suspensions at low and high concentrations, making it suitable for industrial simulations.

The remaining part of this Ph.D. dissertation is structured as follows:

Chapter 2: General Aspects of Rheology and Fluid Dynamics

In this introductory chapter we provide simple physical explanations of the phenomena under study, including fluid, flow, viscosity, and forces. We then present the various tools that may be useful for a mathematical treatment of viscous flows under the continuum hypothesis and in the absence of thermal changes, including forces and stresses, mass and flow equations, and the constitutive equation. While our goal is not to provide a comprehensive presentation of fluid mechanics and rheology, we review the basic knowledge required to understand the developments presented in the following chapters.

Chapter 3: Numerical Study Of Viscoplastic Flows

In this chapter, we propose an innovative methodology for handling the numerical simulation of viscoplastic flows using a multigrid initialization algorithm in conjunction with the SIMPLE procedure. We solve the governing equations for incompressible flow on a collocated grid using finite volume discretization and Rhie and Chow interpolation for pressurevelocity coupling. Using the square lid-driven cavity flow and pipe flow as test cases, we validate our approach and demonstrate that it outperforms the classic single grid method, achieving a computation speed-up ratio of 6.45 for the finest grid size (320 × 320). Our experiments also show that the representation of the yield surface and the plug zone is influenced by the chosen stress growth parameter values, which have implications for the development of rigid areas and the features of the vortices within the flow domain.

Chapter 4: Rheological Characterization of Phosphate Slurry

In this chapter, we investigated the rheological properties of phosphate ore slurries, which are commonly used in fertilizer and phosphoric acid production. We developed experimental protocols to determine the main characteristics of these complex fluids and established relevant rheological models to simulate the numerical flow in a cylindrical pipeline. Our findings suggest that increasing the concentration of solids in these slurries can be achieved without affecting their desired viscosity, which could have practical implications for their processing and transport.

Chapter 5: A New Rheological Model for Phosphate Slurry Flows

In this chapter, we propose a new rheological model for the flow of phosphate-water suspensions and evaluate its ability to replicate their rheological characteristics using rheometric tests. We present a comprehensive methodology for obtaining the model parameters and incorporate it into the OpenFoam numerical code. Our results show that the proposed model can accurately reproduce the rheological behavior of phosphate suspensions at both low and high concentrations, which makes it suitable for simulating and monitoring phosphate slurry flows in industrial applications.

Chapter 6: Concluding remarks and perspectives

The conclusion highlights the main findings of this Ph.D. thesis. Moreover, several future research directions are discussed.

Chapter 2

General Aspects of Rheology and Fluid Dynamics

In this chapter, we shall present the basic tools for a mathematical approach to phosphate slurry rheology and dynamics. Within a certain range of solid concentrations, phosphate slurry or phosphate suspensions flow materials undergo very large, approximately continuous deformations. As a result, the most immediate way of describing these phenomena is to consider them within the frame of fluid mechanics. However, the materials under consideration generally offer resistance (viscosity) to shear, and thus to flow, which is much larger than water flow resistance and is a priori unknown. Consequently, we must start by studying the mechanical behaviour of this kind of material in greater detail. Indeed, this is the aim of rheology. Relationships between stresses and deformations are the fundamental concepts of continuum mechanics, which are discussed in this chapter.

In the following pages, we shall start by providing some simple physical explanations of the phenomena under study: fluid, flow, viscosity, forces, etc. Then we shall present the various tools that may be useful for a mathematical treatment of viscous flows under the continuum hypothesis and in the absence of thermal changes: forces and stresses, mass and flow equations, and the constitutive equation. It is not our aim to provide a complete presentation of fluid mechanics and rheology. We shall restrict ourselves to reviewing the basic knowledge required to understand the developments presented in the following chapters. Readers wishing to explore these fields in greater depth should consult various references such as the one written by [START_REF] Coussot | Mudflow Rheology and Dynamics[END_REF], [START_REF] Stickel | Fluid Mechanics And Rheology Of Dense Suspensions[END_REF] and [START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF] for rheology, and [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF], [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF] for both aspects in conjunction with an overview of the flow behaviour and solution of all types of materials.

Generalities

Stress Tensor

It is possible to classify the applied forces on a fluid element of volume dV in two categories according to their ranges:

-Volume forces that are long range forces (gravity, electrical or magnetic forces). These forces when applied on a volume element are proportional to its volume.

-Surface forces that are short-range forces. These short-range interactions will only involve a thin outer layer on a given volume element.

The overall force applied by these interactions is proportional to the area of the surface limiting the fluid element and is independent of its volume.

Surface forces act on the surface of the fluid element and arise by virtue of contact. In general, the orientation of the contact area may be such that the outward normal to the contact area does not coincide with the direction of the exerted force. Since two vectors (the outward normal vector and the force vector) are involved, a complete description of the surface force exerted on a fluid element requires the definition of the stress tensor (or the stress matrix). Stress is of course, force divided by the area. For an arbitrary macroscopic volume element, the forces acting on its surface element dS, of normal unit n, are due to pressure and viscous stresses which can be expressed in term of the total stress tensor σ, as follow:

dF = σ • ndS (2.1)
The Stress is the quantity that represents applied forces on the surfaces of a volume, which will be the cause of any resulting deformation in that volume. The force vector has three components and the area of contact can be projected onto three planes each normal to a coordinate direction. In the 3D cartesian coordinate system (x, y, z), the Stress σ is a 2 nd rank tensor (i.e. 3 × 3 Matrix) which has 9 components σ ij as you can see in Figure 2.2. The coefficient σ ij of the tensor represents the component in the i-direction of the stress exerted on a surface with normal oriented in the j-direction. Accordingly:

σ yx is the y-component of the force exerted on a unit area with normal pointing in the x-direction. This is a tangential or shear stress. The stress tensor σ (or σ ij ) is symmetric, thus σ ij = σ ji . Hence in cartesian coordinates, the stress tensor at a point can be written as:

σ =   σ xx σ xy σ xz σ yx σ yy σ yz σ zx σ zy σ zz   (2.2)
In a fluid at rest, the force exerted by the fluid element, located on one side of this element, on the fluid located on the other side is normal to the surface elements and its magnitude is independent of the orientation of these elements. Since this stress is isotropic, a single value is sufficient to define it at each point; this is the hydrostatic pressure. Thus, the stress tensor is :

σ ij = -pδ ij (2.3)
where p is the pressure and the negative sign associated with it indicates that a fluid at rest is usually under compression, and thus the stress is acting opposite to the outward normal n. δ ij is the Kronecker delta defined as: -the orientation of the surface dS in space, defined by means of the unit vector n normal to the surface;

δ ij = 1 if i = j 0 if i ̸ = j
-the values of the three components of the force per unit area in the x-, y-, and z-directions for the three orientations of unit surfaces normal to these respective axes.

It is generally very convenient to separate out of the stress tensor σ the part which corresponds to the pressure stresses, which are the only ones acting in the absence of velocity gradients, for a fluid at rest or in uniform translational motion. In practice the stress tensor is split into two terms such that:

σ ij = -pδ ij + τ ij (2.4)
Where τ is deviatoric component known as the deviatoric stress tensor, which tends to distort it. It is also called the shear stress or viscous stress tensor, which is a function of the fluid velocity. The crucial question is then: How is this tensor is related to the deformation (and its history) for a given fluid? The model describing this relationship is termed the stress constitutive equation. Obtaining this relationship is one of the main objectives in rheology.

Strain Rate Tensor

The strain can be seen as the quantity that represents the amount of transformation (deformation) of a body from a reference configuration to a current configuration after a time interval ∆t. The shape, location and orientation of the fluid element at time t + ∆t are determined by the outcome of the following:

-pure translation -pure rotation -shear -volume dilatation.

For the sake of simplicity, and without any loss of generality, we will consider the deformation of the face marked ABCD, see Fig. It can be seen from Fig. 2.3 that both translation and rotation of the fluid element does not result in a deformation of the fluid element. On the other hand, the fluid element undergoes a deformation under a shear or dilatation. As seen in Fig. 2.3, a pure shearing of the fluid element causes a distortion of the fluid element. It is important to note that this distortion occurs without any net rotation, as the diagonals AC and A'C' coincide. However, the edges A'D' and A'B', which were initially perpendicular, are no longer so. Such a shearing strain would result from the action of the σ xy and σ yx components of the stress tensor (see Fig. 2.2) on the fluid element.

Dilatational strain of the fluid element is illustrated in Fig. 2.3 as a an elongation in the x-direction and a compression in the y-direction, resulting in a variation of the volume. It can easily be seen that such a strain would arise from the action of either the σ xx or σ yy (or both) components of the stress tensor. The off-diagonal components of the stress tensor, in general, are the cause of shear strains (which are also the off-diagonal components of the strain rate tensor) and the diagonal components are the cause of dilatational strains (which again are the diagonal components of the strain rate tensor). The quantitative evaluation of these strains is carried out below.

If we consider the face ABCD of the fluid element illustrated again in Fig. 2.4. At time instant t, the velocity components at vertex A are u and v along the x and y coordinate directions respectively. Thus, if we ignore terms of order ∆x 2 and ∆y 2 , the velocity components at the vertices B', C' and D' can be written as: After a time interval of ∆t, vertex A is displaced to a location A' by an amount u∆t and v∆t along the x-and y-directions respectively. Other vertices locations relative to A' are shown in Fig. 2.4. Rotation of the fluid element can be evaluated by the angular velocity with respect to each coordinate axis. With regard to Fig. 2.4, the angular velocity of the fluid element directed along the z-axis can be quantified as the arithmetic mean of the rates of rotation of edges AD and AB. Thus,

Ω z = 1 2 dα dt - dβ dt (2.5)
where the sign convention for angular displacement in the anti-clockwise direction is positive. For small angular motions, we can write:

Ω z = 1 2 ∂v ∂x - ∂u ∂y = ω z 2
The given expression demonstrates the correlation between a fluid element's vorticity and its angular velocity. Analogous expressions can be established for the Ω x and Ω y components.

Strain rate tensor is the variation of the strain tensor with respect to time denoted as γ or, equivalently, as the symmetric part of the gradient (derivative with respect to space) of the flow velocity vector.The shear strain rate of the fluid element can be defined as the average rate of decrease of the angle between two initially perpendicular lines, such as lines AB and AD. Hence:

γxy = 1 2 dα dt - dβ dt (2.6)
Replacing dα/dt and dβ/dt with their expressions in Eq. Strains, like stresses, have normal strain and shear strain components (i.e. perpendicular to or along the face of an element respectively). Normal strain components are γxx , γyy , and γzz in the direction of the axes x, y, and z respectively. The dilatational strain of a fluid element is determined by calculating the rate of change in length along each coordinate direction and dividing it by the original length. Taking edge AB in Figure 2.4 as an example, the dilatational strain in the x-direction is represented as follows:

γxx = 1 ∆t ∆x + ∂u ∂x ∆x∆t -∆x ∆x = ∂u ∂x
If we consider edge A'D', the dilatational strain along the y-direction, γyy , is written as ∂v/∂y. Demonstrating that γzz = ∂ω/∂z is a straightforward extension.

In the 3D cartesian coordinate system (x, y, z) the strain rate γ is also a 2 nd rank tensor which has 9 components and can be finally written as : The rate of strain tensor can be summarized as:

γij = 1 2
∂v i ∂x j + ∂v j ∂x i (2.11)
Where v is the velocity vector, x is coordinate, and the two indices i and j can range over the three coordinates (x, y, z) in three dimensional space. We note in passing that the strain rate tensor is also symmetric, as otherwise, it cannot be related to the stress tensor, which is symmetric.

Relation between the deviatoric stress tensor and the strain rate tensor

Stokes postulates the following with regards to the relationship between the components of the deviatoric stress tensor and the strain rate tensor:

-A linear relationship exists between the components of the deviatoric stress tensor and the components of the strain rate tensor, similar to the Hooke's law in solid mechanics.

-The relationship is unchanged under rotations or exchanges in the coordinate system.

-The deviatoric stress tensor and the strain rate tensor have congruent principal axes.

The last two items in the above list ensure isotropy of fluid deformation in all directions. These postulates lead to the result:.

τ = 2η γ + λ(∇ • u)I (2.12)
In this expression, η represents the well-known dynamic viscosity coefficient. I is the identity matrix (or, the Kronecker delta tensor). It is important to note that the first term in Eq. 2.12 signifies deformation without a change in volume, whereas the second term represents isotropic dilation. This latter term is equal to zero for an incompressible fluid. The second coefficient λ which appears in Eq. 2.12 is commonly referred to as the second viscosity or bulk viscosity. The associated stresses, represented by the diagonal components of the tensor τ in the form λ(∇ • u), are associated with variations in the fluid volume as a result of compression effects. This term becomes insignificant in the analysis of incompressible fluids, as the divergence of velocity ∇ • u becomes equal to zero in that case.

For a Newtonian fluid, we make the assumption that the stress is a linear function of the velocity gradient. In general, the linear relation can be written as:

τ = 2η γ (2.13)
where η is the coefficient of viscosity of the fluid and γ is the symmetric traceless part of the rate of deformation tensor. However, the stress could be a non-linear function of the strain rate in complex fluids such as phosphate slurry.

τ = 2η( γ) γ (2.14)

Concept of viscosity

Viscosity is a measure of the fluid friction or the resistance to flow of a liquid or a gas. This can be illustrated by considering an experiment where the gap between two infinite parallel plates and separated by a distance a in

v x = V 0 y a (2.15)
The resulting flow is commonly referred to as simple shear flow or plane Couette flow. It is of interest to determine the proportional relationship between the frictional force, F , acting on an area S of the plate and opposing the relative motion of the plates, and the velocity gradient between them. The direction of the frictional force F is along the negative x-axis.

F x S = η V 0 L = -η ∂v x ∂y (2.16)
The ratio of F x /S is referred to as the shear stress, which has the dimensions of pressure. The constant η is reffered to as the dynamic viscosity of the fluid, reflecting its association with a force, or simply as viscosity. Its dimensional form is, according to the following:

[η] = [M ][L][T ] -2 [L] -2 [L][T ] -1 [L] -1 = [M ][L] -1 [T ] -1
(2.17)

In the SI system, the unit of dynamic viscosity is the Pascal-second (P a.s) (1P a.s = 1kg/(m.s)).

In many applications it is convenient to employ the combination viscosity/density, denoted ν, which depends on the properties of the fluid, and is known as the kinematic viscosity; it obeys the relationship: The shear rate dependence of fluids is an important consideration, since many fluids have complex shear rate dependence. It is first important to make a distinction between simple 'Newtonian' fluids and 'non -Newtonian' fluids. Newtonian fluids are characterised by the following behaviour:

ν = η ρ (2.
-Viscosity is independent of shear rate.

-Viscosity is independent of time of shear at a constant shear rate.

-The normal stress differences are zero.

-Viscosities measured by different types of deformations, such as uniaxial and biaxial extension, are proportional to each other.

Any fluid showing a deviation from this behaviour is a non-Newtonian fluid. All gases, as well as simple liquids like water, are Newtonian liquids. Non-Newtonian liquids include all materials which have a complex microstructure such as gels, suspensions, phosphate slurry, etc.

Non-Newtonian fluids

Thus far, we have limited our discussion to the case of Newtonian fluids, characterized by a straightforward proportionality between stress and strain rates. However, it is now imperative to discuss the case of non-Newtonian fluids, in which this correlation is no longer linear and may vary based on the fluid's flow history. These unique properties are often a result of the presence of macroscopic objects within the fluid, which are considerable in size compared to the atomic scale but still comparatively small in relation to the overall flow scale. Examples of such objects include macro-molecules in polymer solutions, particles in suspensions, droplets or vesicles in emulsions, and biological fluids. Furthermore, these entities may aggregate to form larger structures that have a significant impact on fluid flow, such as aggregates of platelets in clays, clumps of particles, or entangled macromolecules. Non-Newtonian fluids can be commonly encountered in both natural environments (such as snow, mud, blood, and cream) and in daily life (such as paints, shaving cream, mayonnaise, yoghurt, cosmetics, etc.) as well as in industrial contexts (e.g. phosphate slurry). In order to understand the flow behavior of non-Newtonian fluids, it is crucial to understand their response to an imposed stress. This is the aim of the field of rheology, which was established in the 1920s and named by E.C. Bingham, who is considered, along with M. Reiner, as a founding figure in the discipline. We might recall the statement attributed to Heraclitus: "panta rhei" (everything flows).

As previously stated, a non-Newtonian fluid is characterized by its nonlinear flow curve, which represents the relationship between shear stress and shear rate. This relationship deviates from the linear relationship described by Newtonian fluids and is dependent on a variety of factors, including flow conditions such as flow geometry and shear rate, as well as the kinematic history of the fluid element being considered. Non-Newtonian fluids can be conveniently grouped into three general categories:

-Time-independent or purely viscous fluids, also known as inelastic or generalized Newtonian fluids, for which the rate of shear at any point is determined only by the value of the shear stress at that point at a given instant.

-Time-dependent fluids, for which the relationship between shear stress and shear rate depends, in addition, upon the duration of shearing and their kinematic history.

-Visco-plastic fluids, which exhibit characteristics of both ideal fluids and elastic solids and show partial elastic recovery after deformation.

Among these three general classes of Non-Newtonian fluids, the timeindependent group is the most popular and easiest to analyze. Hereafter, we will exclusively focus on this specific classification of Non-Newtonian fluids.

The rheological properties of a material are typically described by the manner in which the shear stress τ varies with the shear rate γ ( γ = ∂v x /∂y for a simple shear flow). Depending upon the form of this relationship, these fluids may be further subdivided into three types: shear-thinning or pseudoplastic, shear-thickening or dilatant and visco-plastic. Figure 2.6 presents a graphical representation of this relationship for various non-Newtonian fluids on a linear scale. It should be noted that the illustration assumes that the dependence of τ on γ is constant over time or timeindependent. Reproduced from [START_REF] Coussot | Mudflow Rheology and Dynamics[END_REF].

Shear thinning fluids

Fluids that exhibit a reduction in apparent viscosity η app = τ / γ with an increase in shear rate are referred to as shear-thinning fluids or pseudoplastic fluids (the subscript app is omitted in the following). This behavior is commonly observed in dilute solutions of high molecular mass polymers and is attributed to the progressive separation and alignment of entangled macro-molecules as a result of the shear. A considerable number of dilute suspensions of solid particles also exhibit shear thinning behavior, which is the result of the disruption of structures established by interparticle attraction due to flow. Other instances of such suspensions include shampoos and fruit juice concentrates. Similarly, phosphate slurries, which consist of solid particles suspended in liquid matrix, demonstrate comparable characteristics.

Many mathematical expressions of differing complexity have been proposed in the literature to model shear-thinning behaviour. Some of these expressions are simplistic attempts at curve fitting, providing empirical re-lationships between the apparent viscosity (or shear stress) and shear rate, while others are based on theoretical foundations, as the extension of kinetic theory to the liquid state or the theory of rate processes, as informed by statistical mechanics. A selection of the more widely used rheological models is given in the next chapter.

For small values of the shear rate γ, the viscosity of shear thinning fluids, particularly slurries, demonstrates a tendency towards a constant value, referred to as the "Newtonian plateau," denoted as η 0 . As the concentration of suspended particles increases, η 0 increases accordingly and can be several thousand times greater than the viscosity of the suspending liquid, such as water. The relationship between shear stress and shear rate for this type of fluid can be mathematically expressed as follows:

τ = K γn (2.19)
So the apparent viscosity for the so-called power-law fluid is thus given by:

η = K γn-1 (2.20)
In these equations, k and n serve as empirical curve-fitting parameters and are referred to as the fluid consistency coefficient and the flow behavior index, respectively. It is worth noting that for a fluid exhibiting shearthinning behavior, the flow behavior index may range from 0 to 1. The degree of shear-thinning increases as the value of n decreases. Conversely, if the fluid displays shear-thickening behavior, the value of the flow behavior index n will be greater than 1. It is noteworthy that when n = 1, Eq. 2.19 transforms into the constitutive equation of a Newtonian fluid. While the power-law model provides a basic representation of shear-thinning behavior, it is important to acknowledge its limitations. Specifically, it is typically valid only within a restricted range of shear rates, which results in the fitted values of K and n being dependent on the chosen range of shear rates. Additionally, the power-law model fails to predict the zero and infinite shear viscosities.

Shear thickening fluids

The shear thickening fluids don't have a flow threshold and are referred to as dilatant fluids since their viscosity increases with applied stress. This type of fluid behavior was first observed in concentrated suspensions, and a possible explanation for their dilatant behavior is as follows: when at rest, the void space is minimal, and the available liquid is sufficient to fill it. At low shear rates, the liquid functions as a lubricant, allowing each particle to move past others with minimal stress. At high shear rates, however, the material expands or dilates slightly (as seen in the transport of sand dunes), causing the available liquid to be insufficient to fill the newlycreated voids. This leads to direct solid-to-solid contact, which results in increased friction and higher shear stresses. This mechanism causes the apparent viscosity to rise quickly with increasing shear rate. This behavior can also be observed in certain polymer solutions where the initial entanglement of macro-molecules can be unwound into longer chains by the stresses induced by fluid flow, leading to a higher effective viscosity. Additionally, some fluids exhibit shear thinning followed by shear thickening regimes as the shear rate is increased. This can be explained as polymer chains unwinding and aligning, thereby reducing the apparent viscosity η, until the chains begin to interact and cause an increase in η.

Although the changes in viscosity due to the observed trends can be incorporated into a model for apparent viscosity (η), the generalized constitutive law (Eq. 2.20) is strictly a phenomenological model that attempts to replicate the impact of applied flow on a fluid's internal structure by adjusting the viscosity in response to shear. When a material displays both shear-thinning and shear-thickening, it suggests that multiple mechanisms for stress creation and relaxation are competing with each other. As a result, a simple model like Eq. 2.20 would likely fail to accurately describe even basic flows of these fluids. Consequently, generalized Newtonian models should be reserved for use in complex fluid flows, where it is reasonably believed that the primary dynamics are caused by the changes in viscosity induced by shear, and only in the simplest of flows.

Viscoplastic fluid behaviour

This type of fluid behavior is defined by the presence of a flow threshold (τ 0 ), often called the yield stress, that must be exceeded before the fluid will start to deform or flow. In contrast, the material will elastically deform when the externally applied stress is lower than the yield stress. Once the magnitude of the external stress has exceeded the value of the yield stress, the flow curve may follow a linear or non-linear pattern, but it will not pass through the origin (see Figure 2.6). As a result, in the absence of surface tension effects, such a material will not settle under gravity to form a completely flat free surface. This type of fluid behavior can be explained by the theory that the material at rest is composed of three-dimensional structures that are rigid enough to resist any external stress below a certain threshold value τ 0 . However, when the applied stress exceeds this threshold value, the structure breaks down and the material behaves like a viscous fluid. In some cases, the formation and breakdown of the structure are reversible, meaning that the material can regain its initial yield stress value.

A Bingham plastic fluid is a type of fluid that has a linear flow curve when the applied stress |τ yx | > |τ 0 | and is characterized by having a constant plastic viscosity (the slope of the shear stress versus shear rate curve) and a yield stress. On the other hand, a substance that exhibits a yield stress and a non-linear flow curve on linear coordinates (for |τ yx | > |τ 0 |) is referred to as a yield pseudoplastic material. A notable observation is that a viscoplastic material exhibits a decreasing apparent viscosity as the shear rate increases. When the shear rate is very low, the apparent viscosity is virtually infinite just before the substance yields and starts to flow. Therefore, these materials can be considered to have a specific type of shear-thinning behavior. Figure 2.7 illustrates, in semi-log coordinates, the relationship between shear rate and stress for a mixture of phosphate and water, which represents a phosphate particle suspension. As the stress gradually increases from zero, the material does not flow until the stress reaches a threshold value around 10 Pa. Some concentrated suspensions of solids in a liquid, as well as some polymer solutions, exhibit a flow threshold. If an increasing pressure head is applied to such a fluid in a cylindrical tube, immediately above the yield stress, we observe a solid-like flow of almost all the fluid with a velocity independent of the distance from the walls, resulting in localized velocity gradients near the walls, where the stress required for shear flow is achieved. This is called plug flow. As the stress is further increased, the velocity gradient becomes progressively non-zero throughout the volume. In Section 2.3.3, we provide a quantitative calculation of such profiles.

Over the years, several empirical equations have been suggested through simple curve-fitting methods, but a theory-based model has not yet been established. We discuss the frequently used models for viscoplastic fluids in Section 4.1.2. Figure 4.10 displays the relation between the stress and the shear rate in semi-log coordinates for a mixture of kaolin and water which corresponds to a suspension of clay particles: as the stress is gradually increased from zero, we observe no flow until a threshold value of the order of 10 Pa is reached.

Rheology of suspensions

Rheological behavior of solid-liquid suspensions are dependent on the interactions between the dispersed particles and the fluid matrix which cause the diversion of the shear stress versus shear rate relation from the linearity observed in Newtonian fluids. The type and extent of this diversion will again depend on the strength of the attractive and repulsive forces, and the ratio of the time for reorientation of the components to their original configuration to the timescale of measurement.

The mechanics of dilute and semi-dilute suspensions have been extensively studied and are well understood. However, constitutive equations that describe the relationship between stress and the rate of strain for concentrated suspensions are generally unknown. Consequently, their rheology remains a topic of much research, despite a significant amount of work conducted over the past century. The rheological behavior of fluid suspensions, including shear-thinning and thixotropy, shear-thickening, and yield stresses, has been identified through experimentation [START_REF] Stickel | FLUID MECHANICS AND RHE-OLOGY OF DENSE SUSPENSIONS[END_REF]. This type of material constitutes a substantial portion of non-Newtonian fluids, and will be further discussed in this section.

Dimensional analysis

The dimensional analysis presented here is based on the approach of Krieger (1963Krieger ( , 1972) ) and, more recently, Jomha et al. (1991). This analysis takes into account the suspension's viscosity, which is considered a function of various system parameters:

η = f (a, ρ p , n, η 0 , ρ 0 , kT, γ or τ, t) (2.21)
In this equation, we have a set of particle properties including radius a, density ρ p , and number concentration n. Additionally, there are several suspending medium properties such as viscosity η 0 and density ρ 0 , as well as thermal energy kT . The shear variable can be expressed as either shear rate γ or shear stress τ , and time t is also considered. In suspensions with multiple particle sizes, the average radius should be used for a, and one or more additional terms are required to account for the particle-size distribution. All of the terms in Eq. 2.21 can be expressed in units of mass, length, and time. By creating dimensionless groups, this equation can be simplified to 6 variables, which is obtained by subtracting 3 from the total number of terms [START_REF] Hart | A Review of Technologies for Transporting Heavy Crude Oil and Bitumen via Pipelines[END_REF] as:

η r = f (ϕ, ρ r , P e γ , Re γ , kT, t r ) (2.22)
where

η r = η η 0 , ϕ = 4π 3 na 3 , (2.23
)

ρ r = ρ p ρ 0 , P e γ = 6πη 0 a 3 γ kT , (2.24
)

Re γ = ρ 0 a 2 γ η 0 , t r = tkT η 0 a 3 (2.25)
We selected the shear rate γ instead of the shear stress τ as the shear variable. Equation 2.22 can be further simplified in various significant scenarios. In steady-state, neutrally buoyant systems, the variables ρ r and t r can be disregarded:

η r = f (ϕ, P e γ , Re γ ) (2.26)
Krieger made the assumption that the Reynolds number was negligible (Re γ → 0) for his systems with particle radii of around 1µm, and he proposed a semi-empirical equation for η r = f (ϕ, P e γ ). An alternative approach is to consider non-Brownian systems in which the Peclet number is very large (P e γ → ∞) and use the equation η r = f (ϕ, Re γ ). Some authors, such as [START_REF] Chang | Hydrodynamic transport properties of concentrated suspensions[END_REF], have assumed that both P e γ and Re γ can be neglected and assumed that:

η r = f (ϕ) (2.27)
This implies that the viscosity is a unique value at every concentration and hence the suspensions are Newtonian. Both P e γ and Re γ can be neglected, i.e., Re γ ≲ 10 3 and P e γ ≳ 10 3 , for only a relatively narrow window of shear rates, given values of a, η 0 , and ρ 0 . The size of this "window" scales according to the Schmidt number, Sc = P e γ /Re γ . A suspension may be expected to behave as a Newtonian fluid for greater ranges of shear rate as particle size and fluid viscosity increase, such that Sc ≫ 1 (Figure 2.8). This indicates that the viscosity is a unique constant value for every concentration, and therefore the suspensions exhibit Newtonian behavior. Neglecting both P e γ and Re γ , specifically when Re γ ≲ 10 3 and P e γ ≳ 10 3 , is applicable only within a relatively narrow range of shear rates, provided that a, η 0 , and ρ 0 are given. The width of this range, referred to as the "window," scales based on the Schmidt number:

Sc = P e γ Re γ .
(2.28)

As the particle size and fluid viscosity increase, a suspension is expected to behave like a Newtonian fluid over a greater range of shear rate, such that Sc ≫ 1 (refer to Figure 2.8).

Suspension viscosity

There is a large body of data available in the scientific literature regarding the relationship between viscosity and either shear rate or shear stress in fluid-particle suspensions. When the solids concentration exceeds 0.4 by volume, it is common to observe non-Newtonian behavior. Previous studies have suggested that concentrated suspensions exhibit shear-thinning behavior, with Newtonian limiting behavior at both low and high shear rates, as a convenient assumption. However, it should be noted that this assumption may not hold true in all cases.

Yield stress is also a common rheological property exhibited by many types of suspensions, including colloidal suspensions and emulsions. The determination of a yield stress and its empirical quantification has been a topic of considerable controversy and discussion, as documented in [START_REF] Barnes | The yield stress -a review or ' ' -everything flows?[END_REF] [START_REF] Barnes | The yield stress -a review or ' ' -everything flows?[END_REF] and Coussot (2013) [START_REF] Coussot | Yield stress fluid flows: A review of experimental data[END_REF]. Nevertheless, a visible yield stress has been clearly observed in suspensions, which is signified by a non-zero shear stress without any accompanying deformation over extended experimental time frames or, alternatively, by a viscosity that increases towards infinity as shear rates approach zero. Yield stresses are predominantly evident at extremely high concentrations (ϕ > 0.5). Yield stress behavior in suspensions has been linked to a solid-liquid phase transition by certain authors, as reported in [START_REF] Heymann | On the solid-liquid transition of concentrated suspensions in transient shear flow[END_REF] [START_REF] Heymann | On the solid-liquid transition of concentrated suspensions in transient shear flow[END_REF].

In order to obtain a deeper understanding of the non-Newtonian rheology of suspensions, examining the microstructure of the suspension can be a helpful approach. The microstructure of suspensions containing spherical particles is characterized by the relative positions of the suspended particles with respect to each other. Several types of configurations are possible, including strings, sheets, clumps, and semicrystalline groups. It is widely acknowledged that alterations in the microstructure due to shear forces are responsible for the non-Newtonian behavior of dense suspensions. Based on the analysis of dimensions outlined in the previous section (Section 2.2.1), it can be reasonably inferred that suspensions exhibit shear-thinning behavior when the Peclet number for shear rate (P e γ ) is significant, shear-thickening behavior when the Reynolds number for shear rate (Re γ ) dominates, and behave as Newtonian fluids when both can be neglected (as illustrated in Figure 2.8). This conclusion is reinforced by the findings of microstructural analyses, which suggest that suspensions demonstrate an ordered structure at low shear rates and a disordered structure at high shear rates. However, it is important to note that neither the Peclet number nor the Reynolds number can fully account for suspension rheology. The influence of particle polydispersity, particle roughness, electrostatic forces, and van der Waals forces must be considered in addition to the Peclet and Reynolds numbers when examining the structure of suspensions. While laboratory experiments may attempt to minimize the impact of these non-hydrodynamic factors, even minor effects can significantly alter microstructure (Brady Morris, 1997 [50]).

Several researchers have attempted to establish a correlation between relative suspension viscosity and solids volume concentration. However, early studies did not adequately consider the shear rates and particle polydispersity, leading to a significant amount of variation in these correlations [START_REF] Rutgers | Relative viscosity of suspensions of rigid spheres in Newtonian liquids[END_REF]. Some researchers, such as Storms et al. 1990 [START_REF] Storms | Low shear rate viscosity of bimodally dispersed suspensions[END_REF], exclusively conducted experiments in what they considered the low shear rate Newtonian limit, while others, such as Chang and Powell 2002 [START_REF] Chang | Hydrodynamic transport properties of concentrated suspensions[END_REF], focused on the high shear rate Newtonian limit. Based on intuition, one may observe that the relative viscosity of a suspension will tend towards an infinite value when the volume fraction (ϕ) approaches a maximum value (ϕ m ):

lim ϕ→ϕm η r = ∞ (2.29)
Physically, the maximum value of the volume fraction (ϕ m ) corresponds to the highest possible packing fraction for a given suspension composition and arrangement. Nevertheless, in viscosity models of the form, ϕ m is frequently used as a variable parameter as:

η r = f ϕ ϕ m (2.30)
In this context, the maximum packing fraction serves as a scalar measure of the microstructure of a suspension. As it is dependent on the various parameters affecting suspension microstructure, it is not considered an independent system variable and was therefore not included in the dimensional analysis presented in the previous section. By plotting the relative viscosity (η r ) against the ratio of volume fraction to maximum packing fraction (ϕ/ϕ m ) and using appropriate values for ϕ m , viscosity-concentration plots collapse onto a single curve [START_REF] Chang | Hydrodynamic transport properties of concentrated suspensions[END_REF].

To model the correlation between viscosity and concentration for suspensions, it is crucial to take into account non-Newtonian characteristics. To achieve this, an approach that can be adopted is the inclusion of a term that is dependent on the shear rate in the model. In the case of a Newtonian fluid with a dilute suspension of monodisperse particles (ϕ → 0), Einstein (1956) [START_REF] Einstein | Eine neue Bestimmung der Moleküldimensionen[END_REF] demonstrated:

η r = 1 + [η]ϕ (2.31)
where the constant [η] is variously referred to as the 'Einstein coefficient' or the 'intrinsic viscosity', and takes the value [η] = 2.5 in Eq. (2.31) for spherical particles. [START_REF] Batchelor | The Effect of Brownian Motion on the Bulk Stress in a Suspension of Spherical Particles[END_REF] [START_REF] Batchelor | The Effect of Brownian Motion on the Bulk Stress in a Suspension of Spherical Particles[END_REF] expanded this relationship to the second order, given by,

η r = 1 + [η]ϕ + Bϕ 2 (2.32)
In Eq 2.32, B = 6.2 for Brownian suspensions in any flow, and B = 7.6 for non-Brownian suspensions. Their formulas are adapted for monodisperse particles for maximum solid fractions of 0.15 to 0.25 and link the increase of the bulk viscosity of the suspension η to the viscosity of the liquid η 0 , and the solid fraction in the suspension [START_REF] Thomas | Transport characteristics of suspension: VIII. A note on the viscosity of Newtonian suspensions of uniform spherical particles[END_REF]. However, linearity is restricted to low solid fractions ϕ. For dense or highly concentrated suspensions, the apparent relative viscosity (η r = η/η 0 ) increases rapidly and without bound with ϕ. At large solid fraction, small changes in concentration will lead to large changes in viscosity. If non-hydrodynamic inter-particle forces are absent, suspensions with high solid fraction can be described with the satisfactory model of Kriege (1959) [START_REF] Krieger | A Mechanism for Non-Newtonian Flow in Suspensions of Rigid Spheres[END_REF] , valid for low and high shear rates,

η r = 1 + ϕ ϕ m -[η]ϕm (2.33)
One additional input parameter, defined as the maximum packing fraction of particles ϕ m and related to particles properties, is required in this model. Many other functions have been proposed and provide excellent fits to experimental data [START_REF] Mooney | The viscosity of a concentrated suspension of spherical particles[END_REF][START_REF] Frankel | On the viscosity of a concentrated suspension of solid spheres[END_REF][START_REF] Chong | Rheology of concentrated suspensions[END_REF]. As the solid volume fraction ϕ increases, the apparent viscosity η also increases. However, this increase is not uniform across different maximum packing fraction ϕ m values. The rate of viscosity increase is dependent on the ϕ/ϕ m ratio. The maximum packing fraction ϕ m is significantly influenced by particle shape and size distribution since smaller particles can fill the gaps between larger particles. The maximum packing fraction of polydisperse particle systems is higher than that of monodisperse particle systems [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF][START_REF] Ouchiyama | Porosity of a Mass of Solid Particles Having a Range of Sizes[END_REF]. Many experimental studies have been conducted to establish a relationship between ϕ m and suspension viscosity, particularly for bidisperse and polydisperse suspensions [START_REF] Chong | Rheology of concentrated suspensions[END_REF][START_REF] Stickel | FLUID MECHANICS AND RHE-OLOGY OF DENSE SUSPENSIONS[END_REF]. Generally, these studies demonstrate that the apparent viscosity η decreases with an increasing maximum packing fraction ϕ m at a constant volume fraction ϕ.

Another form for the relationship between suspension viscosity and shear rate is due to [START_REF] Krieger | A Mechanism for Non-Newtonian Flow in Suspensions of Rigid Spheres[END_REF] and [START_REF] Cross | Kinetic interpretation of non-newtonian flow[END_REF]:

η r = η 0 + η 0 -η ∞ 1 + aD n (2.34)
where D = τ for [START_REF] Krieger | A Mechanism for Non-Newtonian Flow in Suspensions of Rigid Spheres[END_REF], and D = γ for [START_REF] Cross | Kinetic interpretation of non-newtonian flow[END_REF]. This model can generally be used to fit a shear-thinning viscosity.

Maximum Packing Fraction

Recent research has shown that the maximum packing fraction (ϕ m ) of suspend particles is a crucial parameter that influences the rheological behavior of suspensions (e.g., [START_REF] Stickel | FLUID MECHANICS AND RHE-OLOGY OF DENSE SUSPENSIONS[END_REF]). When the volume fraction is below ϕ m , the deformation of the suspension primarily involves the flow of liquid between particles. However, as the volume fraction (ϕ) increases, the average inter-particle distance decreases, leading to greater gradients in the velocity of the interstitial liquid and a higher viscosity of the bulk suspension. Despite considerable research, the a priori prediction of the maximum packing fraction (ϕ m ) for a system of particles remains an open question. By considering various geometric arrangements of monomodal spheres, a range of theoretical values for ϕ m can be calculated, ranging from the simple cubic value of 0.524 to the hexagonal close-packed value of 0.740 [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF]. However, in practice, a well-mixed suspension does not self-assemble into one of these theoretical arrangements but instead forms a so-called random close-packed (RCP) arrangement. The value of ϕ m for RCP was determined to be 0.64 using settling experiments [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF].

For multimodal systems, it is more difficult to arrive at a theoretical value of ϕ m . Qualitatively, small spheres may fit into the spaces between packed large spheres. Experiments of bimodal systems bear this out: ϕ m increases with size ratio up to about 10:1, at which point the small spheres can completely fit into the empty spaces of the packed large spheres [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF]. The volume fraction ratio between large and small spheres is also important, with maximum packing obtained at about 60%-75% large particles [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF][START_REF] Shapiro | Random packings of spheres and fluidity limits of monodisperse and bidisperse suspensions[END_REF]. Trimodal, multimodal, and polydisperse systems can obtain even higher packing fractions. Rheological experiments with bimodal, multimodal, and polydisperse systems emphasize the relation between ϕ m and suspension viscosity. As expected, suspensions with multimodal distributions exhibit lower viscosities than their monodisperse counterparts at equal total particle volume fractions [START_REF] Chang | Effect of particle size distributions on the rheology of concentrated bimodal suspensions[END_REF][START_REF] Chong | Rheology of concentrated suspensions[END_REF]. [START_REF] Wildemuth | Viscosity of suspensions modeled with a shear-dependent maximum packing fraction[END_REF] [START_REF] Wildemuth | Viscosity of suspensions modeled with a shear-dependent maximum packing fraction[END_REF] put forward a formula that explicitly links the maximum packing fraction, ϕ m , to the applied shear stress. They derived an equation in which ϕ m varies from an initial value of ϕ m 0 at zero shear stress to a maximum value of ϕ m∞ as the shear stress tends to infinity. The proposed relationship between ϕ m and shear stress is an explicit flow-dependent expression and is given by:

ϕ m = 1 ϕ m 0 - 1 ϕ m 0 - 1 ϕ m∞ 1 1 + τ -m -1 (2.35)
Once substituted into an appropriate η r versus ϕ ϕm correlation, this results in a model that expresses the viscosity as a function of particle concentration and shear rate. Eq. 2.35 for ϕ m (τ ) predicts shear-thinning behavior and a yield stress. If the concentration of a suspension exceeds ϕ m 0 but is lower than ϕ m∞ , such as ϕ 0 , it will not undergo deformation unless subjected to a shear stress greater than that needed for ϕ m (τ ) to reach ϕ 0 . This critical stress is known as the yield stress and is clearly dependent on the concentration of the suspension.

The Equations of fluid motion

In this section, the equations that govern the motion of an incompressible fluid, also known as Navier-Stokes equations, are derived. The governing equations of incompressible fluid flow represent mathematical statements of the conservation laws of physics:

-The mass of a fluid is conserved -The rate of change of momentum equals the sum of the forces on a fluid particle (Newton's second law)

Conservation of mass

The principale of conserving mass states that if there are no sources or sinks of mass present, a given area will maintain its mass at a local level . When we express the overall balance of fluid mass inside a fixed volume, we derive a local equation for the conservation of mass in that particular region. This equation is commonly referred to as the continuity equation.

dm dt = d dt V (ρdV ) = V ∂ρ ∂t dV = - S ρu • ndS (2.36)
The vector unit n, normal to the enclosed surface S, is directed outward from the contained volume V . As V is constant, the order of integration and differentiation with respect to time can be exchanged. Furthermore, through the application of Gauss' divergence theorem to the final term of Equation 2.36, we arrive at the following result:

V ∂ρ ∂t + ∇ • (ρu) dV = 0 (2.37)
Given that this equation is applicable to any arbitrary volume V , it follows that the integral must be identically zero, resulting in the equation of continuity:

∂ρ ∂t + ∇ • (ρu) = 0 (2.38)
One can write the ∇ • (ρu) term of Eq. 2.38 as:

∂ρ ∂t + u∇ρ + ρ∇ • u = 0 (2.39)
It is important to acknowledge that the quantity enclosed in the parentheses denotes the variation (dρ/dt) of the fluid density over time for an element of fluid that is tracked along with the flow (the convective derivative that corresponds to the Lagrangian description). Consequently, it is possible to reformulate Eq. 2.38 in the following form:

dρ dt + ρ∇ • u = 0 (2.40)
When there are no considerable variations in absolute pressure or temperature, it is reasonable to assume that the flow of the fluid is incompressible. This means that changes in pressure do not significantly affect density. While this is generally true for liquids, it can also be a good approximation for gases traveling at speeds significantly lower than the speed of sound. It is worth noting that sound waves are compressible. As a result, the application of the mass conservation (continuity) equation in fluid dynamics cannot be used to determine the density in such scenarios.

The condition of incompressibility implies that the density (ρ) of the fluid does not vary with the flow, and this can be represented mathematically as dρ/dt = 0. In accordance with the mass conservation equation provided in Eq. 2.41, it can be inferred that the continuity equation for incompressible flow can be expressed as follows:

∇ • u = 0 (2.41)
It is noteworthy that although dρ/dt = 0 does not necessarily signify that the density ρ of the fluid is uniform throughout (although it may be the case in several hydraulic applications), it does imply that the density does not vary along a streamline. To be more precise, the incompressibility assumption dictates that each fluid element maintains its initial density as it moves.

Conservation of momentum

The principle of the conservation of momentum asserts that when an external force is absent, an body maintains its total momentum, defined as the product of its mass and velocity vector. Because momentum is a vector quantity, the conservation of its components holds in all directions.

The dynamics of a fluid volume, represented by V , can be analyzed using Newton's second law of motion. This involves equating the time derivative of momentum with the net forces exerted on the volume, including those from its surface and volume. The matter within V moves with the fluid and remains confined within it as:

d dt V ρudV = V ρf dV + S σ • ndS (2.42)
σ is the tensor of all the forces, namely pressure and viscosity, applied on the differential surface element dS. Conversely, the volume force per unit mass of fluid f may comprise a range of physical effects, including gravitational or electrostatic forces on charged fluids, Coriolis forces in rotating reference systems, or magnetic forces exerted upon fluids containing suspended magnetic particles, known as ferrofluids.

It is worth noting that the derivative (d/dt) refers to the Lagrangian derivative, calculated within a reference system that moves with the fluid. Within this reference frame, the product ρdV , denoting the mass of a small fluid element, remains constant. This occurs because any fluid element, by definition, has the same molecules as they traverse through the local velocity field of the flow. Therefore, it becomes possible to derive the factor u in the first term of Eq. 2.42 with respect to time. The resulting expression can be represented as follows:

d dt V ρudV = V ρ du dt dV (2.43)
The last term in Eq. 2.42 can be transformed into a volume integral by means of Gauss' divergence theorem. Eq. 2.42 can therefore be written as:

V ρ du dt dV = V ρf dV + V ∇ • σdV (2.44)
In this context, the term ∇ • σ denotes a vector consisting of components ∂σ ij /∂x j . It is also worth noting that in Eq. 2.44, the integrals are computed across a volume V that moves concurrently with the fluid. By taking the limit as this volume approaches zero and dividing by the value of the volume element, we obtain the local equation of motion for a fluid particle:

ρ du dt = ρf + ∇ • σ (2.45)
rmal stresses. We then have:

We can now proceed to partition the tensor σ into two components: one corresponding to pressure and the other to viscous forces. This separation mirrors the procedure outlined in Eq. 2.4, where we stated: σ ij = -pδ ij + τ ij ; this expression is a reminder that pressure only contributes to normal stresses. Thus, Eq. 2.45 becomes:

ρ du dt = ρf -∇p + ∇ • τ (2.46)
This equation is applicable to any fluid, as we have made no assumption about the form of the stress tensor τ . Most often, the acceleration of the fluid particle is be expressed as :

du dt = ∂u ∂t + (u • ∇)u (2.47)
Thus, Eq. 2.48 becomes:

∂u ∂t + (u • ∇)u = ρf -∇p + ∇ • τ (2.48)

One-dimensional flow of viscous fluids

This subsection reviews mathematical models of the flow of time-independent, non-Newtonian, viscous fluids through a straight pipe in laminar flow regime. The fluid will be regarded as homogeneous. The equations of particular relevance are given in borders, and the other equations provide useful analytical solutions to which numerical approximations can be compared later.

Consider the flow of fluid through a straight pipe as illustrated by the figure below.

Steady-incompressible Navier-Stokes Equations

We begin with a statement of the governing equations, the steady, incompressible Navier-Stokes equations in polar coordinates, in the absence of body-force terms. We list these here as :

∂u ∂z + 1 r ∂ ∂r (rv) = 0, (2.49a)
ρ(u ∂u ∂z + v ∂u ∂r ) = - ∂p ∂z + [ ∂ ∂z (µ ∂u ∂z ) + 1 r ∂ ∂r r(µ ∂u ∂r )] (2.49b) ρ(u ∂v ∂z + v ∂v ∂r ) = - ∂p ∂r + [ ∂ ∂z (µ ∂v ∂z ) + 1 r ∂ ∂r r(µ ∂v ∂r )] (2.49c)
Note that there are no θ-direction derivative terms due to the axisymmetric assumption. If we invoke the fully-developed flow assumption so that u z = 0 and v = 0, these equations can be readily reduced to:

∂u ∂z = 0, (2.50a) 1 r ∂ ∂r r(µ ∂u ∂r )] = ∂p ∂z (2.50b) ∂p ∂r = 0 (2.50c)
we observe from Eq.(2.50c) that pressure does not depend on the radial coordinate, or more formally,

p = C(z)
In particular, pressure can depend only on the z direction. Differentiation of the above with respect to z gives

∂p ∂z = dC dz
indicating that the z component of the pressure gradient can depend only on z. Now from the fact that the flow is fully developed it follows that u can be a function only of r, implying that the left-hand side of Eq. (2.50b) can depend only on r. But we have just seen that the right-hand side, ∂p/∂z, depends only on z. Thus, this pressure gradient must be a constant, and we set

∂p ∂z = ∆p L (2.51)
Substitution of Eq 2.51 into Eq. (2.50b) yields

1 r ∂ ∂r r(µ ∂u ∂r ) = 1 r ∂ ∂r rτ = ∆p L (2.52)
Using the initial condition τ = 0 at r = 0, the expression above integrates to give

τ = r 2 ∆p L (2.53)
At r = R, this gives the following relationship between wall stress and pressure gradient:

τ w = R 2 ∆p L (2.54)
The above equation pair can be combined to give a linear relationship between shear stress and radial distance, which is

τ rz τ w = r R (2.55)
Since Equations 2.50 are in terms of shear stress, the equations so far make no assumptions about the rheology of the fluid being used. These equations are valid for any fluid, and will often be referred to in the following discussion of laminar, critical and turbulent flow.

Dimensionless Forms

Flow equations can be expressed as relationships between dimensionless groups that are useful for defining graphical or empirical relationships of a given flow situation. Dimensionless groups are more than just a convenient way of expressing flow relations as they can, for instance, be used to identify the critical point that divides the laminar flow regime from the transitional/turbulent flow regime.

The first step in this process is identification of independent and dependent variables, and parameters, that fully describe the system. Once this has been done, we introduce "typical values" of independent and dependent variables in such a way as to render the system dimensionless. Then, after some rearrangement of the equations, the dimensionless parameters that characterize solutions will be evident, and it is these that must be matched between flows about two geometrically similar objects to guarantee dynamic similarity. We will demonstrate this scaling procedure using the 2D incompressible continuity and Navier-Stockes equations:

∂u ∂x + ∂v ∂y = 0, (2.56a)
∂u ∂t + u ∂u ∂x + v ∂u ∂y = - 1 ρ ∂p ∂x + µ ρ ( ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ) (2.56b) ∂v ∂t + u ∂v ∂x + v ∂v ∂y = - 1 ρ ∂p ∂y + µ ρ ( ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 ) -g (2.56c)
The independent variables of this system are x, y and t; the dependent variables are u, v and p, and the parameters are g, the gravitational acceleration in the y direction (taken as constant), density ρ and viscosity µ, also both assumed to be constant. We next introduce "typical" values of independent and dependent variables needed to make the equations dimensionless. These values must be chosen by the analyst, and experience is often important in arriving at a good scaling of the equations. Here we will demonstrate the approach with a simple flow in a pipe as depicted in Fig. 2.12. In this case we have indicated a typical length scale to be the height H of the duct, and we have taken the velocity scale to be the centerline speed U c (which is the maximum for these types of flows).

It is often the case that the correct time scale can be obtained the length and velocity scales:

u L T =⇒ t s = H Uc
where t s is the time scale for this flow situation. We introduce the notation P s for the pressure scale. We can now formally scale all independent and dependent variables:

x * = x/H, y * = y/H, t * = t/t s ,
and

u * = u/U c , v * = v/U c , p * = p/P s ,
Notice that the " " quantities are all dimensionless.

Substitution of the dependent variables into Eq. (2.56a) and noting that U c is constant yields

∂u * ∂x * + ∂v * ∂y * = 0 (2.57)
We can now apply the same procedure to the x-momentum equation. Analogous to what we have just done with some rearrangements, we can write,

∂u * ∂t * + u * ∂u * ∂x * + v * ∂u * ∂y * = - P s ρU 2 c ∂p * ∂x * + µ ρU c H ( ∂ 2 u * ∂x * 2 + ∂ 2 u * ∂y * 2 ) (2.58)
We set P s = ρU 2 c the coefficient on the pressure gradient term in Eq. (2.58). The quantity ρU 2 c occurs widely in fluid dynamics; it is two times what is termed the dynamic pressure, denoted p d : that is,

p d = 1 2 ρU 2 c
the reciprocal of the term µ/ρU c H is called the Reynolds number after Osbourne Reynolds who identified it as a key parameter in his studies of transition to turbulence. In general we express the Reynolds number for a Newtonian fluid as:

Re = ρU H µ (2.59)
and it is a measure of the ratio of inertial forces, ρU 2 /D, to viscous forces, µU/D 2 . In 1883,Reynolds conducted some classical experiments of Newtonian fluid flow through pipes, and showed that fluid flow may either exhibit laminar or turbulent motion. For low Reynolds numbers of less than about 2 100, the viscous forces dominate the inertial forces, and laminar flow prevails. However, for high Reynolds numbers greater than about 4000, the inertial forces dominate the viscous forces and turbulent flow prevails. For intermediate Reynolds numbers, fluid flows in a transitional state between laminar and turbulent flow.

If we now suppress the " " notation we can express Eqs. (2.58) in dimensionless form as

∂u ∂x + ∂v ∂y = 0, (2.60a
)

∂u ∂t + u ∂u ∂x + v ∂u ∂y = - 1 ρ ∂p ∂x + 1 Re ( ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ) (2.60b) ∂v ∂t + u ∂v ∂x + v ∂v ∂y = - 1 ρ ∂p ∂y + 1 Re ( ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 ) - 1 F r 2 (2.60c)
where r is the Froude number, defined as

F r = U √ gH (2.61)
We shall now introduce the Fanning friction factor is a dimensionless group defined as the following ratio between frictional forces to inertial forces:

f = τ w ρU 2 2 (2.62)
Notice that this ratio is indicative of the relative importance of the wall stress with respect to the kinetic energy per unit volume of fluid.

Pressure Loss through Straight Circular Pipe

For steady slurry flows in a piping system with straight conduit of uniform cross section, the total pressure change, ∆p, in general comprises,

∆p = p f + p f t + ρg∆z (2.63)
where p f is the frictional pressure loss in the pipe, p f t is the frictional pressure loss resulting from flow through fittings. The last term, ρg∆z, is the pressure loss or gain resulting from elevation changes and is negative for a reduction in static head and positive for an increase in static head z is the net change in pipeline elevation, assuming that the pipeline is full of slurry.

The frictional pressure loss over a straight pipe with a length L and an internal diameter D is given by

p f = f 4L D ρ U 2 2 (2.64)
where f is the friction factor of pipe flow. This equation shows that the evaluation of f is the major task in the determination of the frictional pressure loss. Note that the friction factor can also be defined in terms of the shear stress at the pipe wall, τ w , as shown in Eq. (2.62).

Flow rate

Consider a one-directional flow of fluid through a circular tube with radius R (or D/2), Figure 3.5. The volumetric flow rate through an annular element of area perpendicular to the flow and of width δr is given by :

δQ = 2πrδru (2.65)
and, consequently, the flow rate through the whole tube is

Q = 2π R 0 rudr (2.66)
Integrating by parts yields

Q = 2π r 2 2 u R 0 + R 0 r 2 2 - du dr dr (2.67)
If there is no slip at the pipe wall (u(R)=0), the first term in Eq. (2.67) vanishes, then can be written as

Q = π R 0 r 2 ( γ)dr (2.68)
Changing variables in Eq. (2.67), using Eq. (2.55), and dropping the subscripts rz, Eq. (2.67) can be written as

Q = π τw 0 τ 2 R 2 τ 2 w γ R τ w dτ = π R 3 τ 3 w τw 0 τ 2 γdτ
Multiplying both sides of the above equation with 4 and rearranging, we get

4Q πR 3 = ϕ = 4 τ 3 w τw 0 τ 2 γdτ (2.69)
The shear rate term, , can be eliminated from this equation with a relevant rheological equation. This gives a relation, known as the pseudoshear flow function, between the pipe wall stress and the pseudo-shear rate. The pseudo-shear flow function is a pipe flow equivalent of the true shear flow function. For Newtonian fluids, using the equation (τ = µ γ) to eliminating γ from Eq. (2.69) gives,

τ w = µϕ
It therefore follows that ϕ = γw for Newtonian fluids, though this is not generally true for other fluids. Using Eq. (2.70) and (2.59), a well-known relationship between f the friction factor and Re is revealed by the equation above, and is given by

f = 16 Re (2.70)
For any pseudoplastic fluid, the shear rate γ can be obtained from Eq. (2.69) to give the Mooney-Rabinowitsch equation by differentiating with respect to τ w . It is convenient first to multiply Eq. (2.69) by 3τ w , then differentiating with respect to τ w , we get: where

3τ 2 w ϕ + τ 3 w dϕ dτ w = 4τ
n ′ = dlnτ w dlnϕ
This equation gives a relation between the wall shear rate and the pseudo-shear rate where n ′ can be estimated by plotting ln(ϕ) against ln(τ w ) and measuring the gradient. Thus from a practical viewpoint, values of n' can be estimated from a log-log plot of the pipe flow data, and used to estimate corresponding values of γw .The measurements and the calculation procedure are as follows:

• Measure Q at various values of ∆P/L.

• Calculate τ w from the pressure drop measurements and the corresponding values of the flow characteristic ϕ = 4Q πR 3 from the flow rate measurements.

• Plot ln(ϕ) against ln(τ w ) and measure the gradient at different points on the curve.

• Calculate the pseudo-shear rate at the wall from Eq. (2.71) with n'.

Calculation of flow rate for laminar flow using shear stress -shear rate data

Depending on the flow information available, flow rate-pressure drop calculations for laminar non-Newtonian flow in a pipe ca be performed in different ways. When the flow data are in the form of flow rate and pressure gradient measured in a tubular viscometer or in a pilot scale pipeline, direct scale-up can be done. When the data are in the form of shear stressshear rate values, the flow rate can be calculated directly using Eq. (2.69), where D is the diameter of the pipe to be used and τ w is the wall shear stress corresponding to the specified pressure gradient. Whether obtained with a rotational instrument or with a tubular viscometer, the data provide the relationship between τ and ġamma. Numerical evaluation of the integral in Eq. (2.69) may be done using selected pairs of values of τ and γ ranging from 0 to τ w .

When shear stress -shear rate data relationship can be accurately represented by a simple algebric expression, such as the Casson model or the rheological model that we propose, over the required range then this may be used to substitute for γ, in Eq. (2.69), allowing the integral to be evaluated analytically.

At a given pressure drop, the wall shear stress is given by Eq. (2.54):

τ w = D 4 ∆p L (2.73)
The equation above gives as a value of τ w that is necessary to evaluate the integral in Eq. (2.69) analytically,

4Q πR 3 = 4 τ 3 w τw 0 τ 2 γdτ
For Newtonian slurries,

Q = πD 3 τ w 32µ = πD 4 ∆p 128µL (2.74)
For Bingham fluids,

Q = πD 3 τ w 32η B 1 - 4 3
τ yB τ w + 1 3
τ yB τ w 4 (2.75)
For Casson fluids,

Q = πD 3 τ w 32η c 1 + 4 3 τ yc τ w - 16 7
τ yc τ w 1/2 - 1 21
τ yc τ w 4 (2.76)
For Herschel-Bulkley fluids,

Q = πD 3 8τ 3 w 1 K 1/n n(τ w -τ yH ) (1+3n)/n 1 + 3n + 2nτ yH (τ w -τ yH ) (1+2n)/n 1 + 2n + nτ 2 yH (τ w -τ yH ) (1+n)/n 1 + n (2.77) For Robertson-Stiff fluids, Q = - π 3 γ0 R 3 + π ∆p 2KL 1/n n 1 + 3n R (1+3n)/n (2.

78)

Velocity distribution

In what follows, the velocity distribution u(r) -the distribution of axial direction velocity against radial distance-will be given. The shear rate is written as γ = -du/dr. Assuming the no-slip condition u(R) = 0, this integrates to give

u(r) = R r γdr (2.79)
For the Newtonian case, substituting (τ = µ γ) and Eq. (2.55) into Eq. (2.79) gives the velocity distribution,

u(r) = 2U 1 - r R 2 (2.80)
Where U is the mean cross-sectional velocity.

Considering the linear relation between τ and r defined by Eq. (2.55) and the rheological models, we obtain the following expressions :

For Bingham fluids,        u(r) = 1 η B R τ w 2 -τ 0 -r τ w r 2R -τ 0 if r > r p u = R η B τ 2 0B 2τ w + τ w 2 -τ 0 if r ≤ r p
Where r p = Rτ y /τ w is the radius of plug flow region, where the velocity is constant. For Casson fluids,

                   u(r) = Rτ w 4η c 16 3
τ 0 τ w 1/2 r R 3/2 -1 + 4 τ 0 τ w 1 - r R + 2 1 -r R 2 if r > r p u = Rτ w 4η c 2 - 2 3
τ 0 τ w 2 - 16 3
τ 0 τ w 1/2 + 4 τ 0 τ w if r ≤ r p (2.81)
For Herschel-Bulkley fluids,

           u(r) = 1 K 1/n R τ w n (n + 1) (τ w -τ 0 ) (1+n)/n - τ w R r -τ 0B (1+n)/n if r > r p u = 1 K 1/n R τ w n (n + 1) (τ w -τ 0 ) (1+n)/n if r ≤ r p For Robertson-Stiff fluids,            u(r) = 1 2K ∆p L 1/b ( n n + 1 ) R (n+1)/n -r (n+1)/n -γ0 (R -r) if r > r p u = 1 2K ∆p L 1/b ( n n + 1 ) R (n+1)/n -r (n+1)/n p -γ0 (R -r p ) if r ≤ r p (2.82)
Where the plug region radius is r p = 2(K γn 0 ) (∆p/L) .

Chapter 3

Numerical Study Of Viscoplastic Flows

Introduction

In light of the advantages of the finite volume method for solving fluid flows, it has been employed within many Computational Fluid Dynamics (CFD) solvers. The most attractive features of the method lie in its potential to produce highly accurate solutions with low computational cost and the possibility to be applied to a variety of physical problems with relative ease. The aim of this work is the development of numerical simulations based on the finite volume method and multigrid initialization procedure, with a view to solving viscoplastic flows. As the main novelty, we present an easy to implement multigrid method that significantly improves the convergence rate in terms of CPU time and number of iterations. To the authors's knowledge, this approach has not been addressed in previous research on non-Newtonian flows. The pressure-velocity coupling lie on a momentum weighted interpolation and a structured collocated grid. Multigrid is implemented with a cell-by-cell relaxation procedure. For solving the flow and pressure field, the popular SIMPLE-like algorithm was chosen as the non-linear solver. First, pressure-velocity coupling and multigrid algorithms for steady-state simulations are described together with a summary of the iterative procedure followed in each of them. In what follows, we examine the accuracy of the numerical solution for several different situations and verify the capability of the regularization approach to describe the flow properties of a liquid with a yield stress. Our test examples can be considered in two general categories. We first perform an accuracy test for our present code using different parameters to establish the qualitative credibility of the obtained results. Secondly, we examine the effect of the growth stress parameter in determining the shape of the yielded surface in a steady flow. Finally, the performance of the multigrid initialization algorithm in different situations is discussed.

Numerical methodology

Governing equations

We consider the steady-state, two-dimensional flow. Let ρ and η = η( γ) denote the density and the viscosity of any generalized-Newtonian fluid.

Incompressible steady Navier Stokes equations are given by a pair of partial differential equations:

∇.u = 0 (3.1) ρu.∇u = -∇p + ∇.τ + f (3.2)
Eq. (3.1) represents conservation of mass whereas Eq.(3.2) represents conservation of momentum. The velocity field is defined as u(x, t), with components u and v. The Cauchy stress tensor σ is introduced as the sum of isotropic and deviatoric parts, σ = -pI + τ . Here, the pressure p(x, t) is multiplied by the identity tensor, while the deviatoric part of the stress tensor is denoted as τ (x, t). We have introduced f to describe external body forces such as gravity acting on the fluid. Note that the mass conservation equation is simplified to (3.1) due to incompressibility, i.e., constant density within a fluid volume.

Newtonian flow is characterized by a dynamic coefficient of viscosity µ > 0 which is independent of the strain, yielding a linear relationship between rate-of-strain and stress in the rheological equation in a simple shear flow,

τ = µ γ (3.3)
For fluids where the dependency of the stress on the shear rate tensor is nonlinear, the apparent viscosity is a useful concept when considering rheological responses to strain. It is a generalization of the constant viscosity for Newtonian flow, where we allow the viscosity to be a function of the magnitude of the shear rate tensor. Denoting the apparent viscosity by η, we thus have:

η = τ γ (3.4)
When considering viscoplastic fluids, Eq. (5.1) does not present any problems for yielded regions τ > τ y , but the corresponding apparent viscosity η = µ + τ y / γ has a singularity when the stress fall below the yield (and τ → τ y ). To help address this issue, we employ the popular Papanastasiou regularization [START_REF] Papanastasiou | Flows of materials with yield[END_REF], which lies on an exponential relaxation according to:

1 γ → 1 -e -m γ γ (3.5)
When (m γ) >> 1, it is a suitable approximation that provides a continuous dependence between the shear and shear rate. While near the small shear rate magnitude limit, we have:

lim γ→0 1 -e -m γ γ = lim γ→0 1 - ∞ n=1 (-m γ) n! = m (3.6)
Substitution of Eq. (3.5) into the apparent viscosity of Bingham model yields the following regularized viscosity: The effect of m = 1/ϵ is related to the curves shown in Fig. 3.1. For twodimensional flow, the explicit form of the second invariant of rate-of-strain tensor, γ is given by:

η = µ p + τ y γ (1 -e -m γ ) (3.7)
γ = 2 ∂u ∂x 2 + 2 ∂v ∂y 2 + ∂u ∂y + ∂v ∂x 2 1/2 (3.8)
Under these considerations, the conservation equation for a general flow dependent variable ϕ in the steady laminar case described above, can be written as:

∂ ∂x (ρuϕ) + ∂ ∂y (ρvϕ) = ∂ ∂x Γ ∂ϕ ∂x + ∂ ∂y Γ ∂ϕ ∂y + S ϕ (3.9)
where u and v are the x and y components of the velocity field and Γ is the diffusion coefficient. The terms on the left represent the net convection flow. The two terms on the right represent the net diffusion and the last term the source generation. The mass conservation equation is obtained by setting:

ϕ = 1, Γ = 0, S ϕ = 0
Similarly, u-momentum and v-momentum equations can be obtained from Eq. (3.9), respectively, by setting:

ϕ = u, Γ = η, S ϕ = - ∂p ∂x + f x and ϕ = v, Γ = η, S ϕ = - ∂p ∂y + f y 3.2.

Discretization of the equations

The domain is divided into a number of control volumes (CVs) using a Cartesian grid of equally spaced horizontal and vertical grid lines as shown in Fig. 3.2. We denote the unit vectors in the x and y directions by i and j, respectively. In finite volume method, a discrete approximation of the continuity and momentum equations is obtained by integrating each equation over every cell. The resulting algebraic expressions involve the values of the unknowns u and p at the centroid of each cell and at the centers of neighboring CVs. The surface flux integrals are evaluated separately on each face using various interpolation schemes. Fig. 3.2 shows a control volume P and its neighbors, S, E, N and W . The letters P, S, E, N and W also denote the position vectors of the centers of the respective CVs. Integrating Eq. (3.9) over the volume ∆V of the computational cell, in the absence of any source term, yields the following form: By integrating the governing equation over a control volume and applying Green-Gauss theorem, the semi-discretized form of the equation, at its nodal point P, can be written as follows :

(ρuAϕ) e -(ρuAϕ) w + (ρvAϕ) n -(ρvAϕ) s = ΓA dϕ dx e -ΓA dϕ dx w + ΓA dϕ dy n -ΓA dϕ dy s (3.11)
where A is the cross-sectional area of the control volume face. With the aim of obtaining an useful form of the discretized equation, the diffusion coefficients Γ, the mass fluxes ρu and ρv and the gradients dϕ/dx and dϕ/dy at the cell faces e, w, n and s are required. At the control volume faces an approximate distribution of properties between nodal points is used. Following well-established practice, linear approximations seem to be the obvious and simplest way of calculating interface values and the gradients. This practice is called central differencing. In a uniform grid linearly interpolated values for mass fluxes at cell face ′ e ′ is given by:

(ρu) e = (ρu) P + (ρu) E 2 (3.12)
Similarly, the coefficient Γ is obtained by:

Γ e = Γ P + Γ E 2 (3.13)
And the diffusive flux term is evaluated as:

ΓA dϕ dx e = Γ e A e ϕ E -ϕ P δx e (3.14)
It is convenient to define two variables F and D to represent the convective mass flux per unit area and diffusion conductance at cell faces. For x direction it is written as:

F = ρu (3.15a) D = Γ δx (3.15b)
Assuming that A e = A w = ∆x and A n = A s = ∆x and considering Eqs. (3.12), (3.14) and (3.15), the integrated general conservation equation (3.10) can now be written as :

∆y(F e ϕ e -F w ϕ w )+∆x(F n ϕ n -F s ϕ s ) = ∆y (D e (ϕ E -ϕ P ) -D w (ϕ P -ϕ W )) + ∆x (D n (ϕ N -ϕ P ) -D s (ϕ P -ϕ S )) (3.16)
We assume that the velocity field is known, which takes care of the mass fluxes values F . To solve Eq. (3.16), we need to calculate the transported property ϕ at each cell face using appropriate numerical schemes. In the discretization of divergence terms, many schemes are available in literature (first-order upwind, second-order upwind, central differencing scheme, QUICK, etc). Second and higher-order schemes have been widely used in applications involving orthogonal and uniform meshes. However, the stability of higher-order schemes for applications involving complex flows in complex geometries is not guaranteed and convergence may be difficult to achieve. A stable and fast converging variant of QUICK schemes can be summerized as follows:

ϕ w = ϕ W + 1 8 (3ϕ P -2ϕ W -ϕ W W ) if F w > 0 ϕ e = ϕ P + 1 8 (3ϕ E -2ϕ P -ϕ W ) if F e > 0 ϕ w = ϕ P + 1 8 (3ϕ W -2ϕ P -ϕ E ) if F w < 0 ϕ e = ϕ E + 1 8 (3ϕ P -2ϕ E -ϕ EE ) if F e < 0 (3.17)
First or second-order methods are best suited for such applications. The first order upwind differencing scheme takes into account the flow direction when determining the value at a cell face: the convected value of ϕ at a cell face is taken to be equal to the value at the upstream node.

For the sake of brevity, we make use of the upwind scheme and consider an one dimensional general transport equation. We also assume that ∆x = ∆y. When the flow is in the positive direction, F w > 0 and F e > 0, the upwind scheme sets

ϕ e = ϕ P , ϕ w = ϕ W (3.18)
and the discretized Eq. (3.16) becomes

F e ϕ P -F w ϕ W = D e (ϕ E -ϕ P ) -D w (ϕ P -ϕ W )
which can be rearranged as

(D w + D e + F e )ϕ P = (D w + F w )ϕ W + D e ϕ E
to give

[(D w + F w ) + D e + (F e -F w )]ϕ P = (D w + F w )ϕ W + D e ϕ E (3.19)
When the flow is in the negative direction, F w < 0 and F e < 0, the scheme yields

[D w + (D e -F e ) + (F e -F w )]ϕ P = D w ϕ W + (D e -F e )ϕ E (3.20)
Identifying the coefficients of ϕ W and ϕ E as a W and a E , Eq. (3. [START_REF] Bingham | An investigation of the laws of plastic flow[END_REF]) and (3.20) can be written in the general form:

a P ϕ P = a W ϕ W + a E ϕ E (3.21)
with central coefficient

a P = a W + a E + (F e -F w ) (3.22)
and neighbor coefficients

a W = D w + max(F w , 0) a E = D e + max(-F e , 0)
Following the same methodology described above, the two-dimensional discretized equation can be written in the following final discretized form:

a P ϕ P = a W ϕ W + a E ϕ E + a S ϕ S + a N ϕ N + b ϕ (3.23)
Where b ϕ is the source term and :

a P = a W + a E + a S + a N + (F e -F w ) + (F n -F s ) a W = D w + max(F w , 0) a E = D e + max(-F e , 0) a S = D s + max(F s , 0) a N = D n + max(-F n , 0)
In order to reach convergence, under-relaxation factors are typically applied to independent variables to limit the change in consecutive iterations. From Eq. (3.24) we can write the final expression of ϕ P as:

ϕ P = α ϕ a P (a W ϕ W + a E ϕ E + a S ϕ S + a N ϕ N + b ϕ ) + (1 -α ϕ )ϕ 0 (3.24)
where superscript 0 is used for the quantities calculated at the previous iteration and α ϕ is the under-relaxation factor for variable ϕ.

Pressure-velocity coupling

Rhie and chow's interpolation

Since momentum equations are integrated over co-located control volumes, checkerboard oscillations may appear if central difference is used to approximate face velocities in discretized continuity equation and pressure gradient in momentum equations [START_REF] Kawaguchi | Checkerboard pressure predictions due to the underrelaxation factor and time step size for a nonstaggered grid with momentum interpolation method[END_REF]. Indeed, in this case, the interpolated expression for face velocity will not include a pressure gradient across the faces of the control volume and the resulting solution can be oscillatory. Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] interpolation is used to prevent spurious oscillations in the solution by ensuring strong pressure-velocity coupling. Their procedure has been extremely successful in co-located grids because of its inherent advantages, such as convenience in implementation on an unstructured grids and economical storage of velocity and pressure data.

By taking out the pressure gradient from the source term in Eq. (3.24), the final discretization equation for x-component of the velocity field u P takes the following form, into which under-relaxation is incorporated:

u P = nb a nb u nb + b u a P -∆y (p e -p w ) a P (3.25)
The two terms p w and p e are linearly interpolated from the neighboring nodes as:

(p w ) P = f + w p P + (1 -f + w )p W (p e ) P = f + e p E + (1 -f + e )p P
where f + e and f + w are weighting factors that can be determined, respectively, as f + e = ∆x/2δx e and f + w = ∆x/2δx w for the structured mesh shown in Fig. 3.2. A similar procedure can be applied for velocity component v p and the related terms p s and p n , and the corresponding interpolation coefficients f + s and f + n . v p is written as:

v P = nb a nb v nb + b v a P -∆x (p s -p s ) a P (3.26)
Eq. (3.25) can also be expressed in general form:

u P = H P -D P (∇P ) P (3.27)
where

Hp = nb a nb u nb + b p a P
and,

D p = V a p
For the sake of brevity, interface velocity interpolation is shown only for face e. From Eq. (3.25), we could write the following equations for two adjacent cells E and P:

u P =   nb a nb u nb + b u a P   P -∆y
(p e -p w ) P (a P ) P (3.28)

u E =   nb a nb u nb + b u a P   E -∆y (p e -p w ) E (a P ) E (3.29)
Following the formulation for u E and u P , Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] proposed similar equation for the interface velocity u e :

u e =   nb a nb u nb + b u a P   e -∆y (p E -p P ) (a P ) e (3.30)
As in the staggered grid, Eq. (3.30) includes the difference of pressure (p E -p P ). As it is, if a checkerboard pressure distribution occurs, pressure will not be properly represented in the discretized momentum equations and the continuity equation will reject it as an acceptable solution. Hence, special treatment is needed. Eq. (3.30) can also be expressed as: yielding the classical expression of Rhie-Chow momentum interpolation, which does not include the under-relaxation factor of velocity. However, in general under-relaxation is required for the momentum equations. When considering under-relaxation, Eq. (3.37) becomes:

u e =
u e = α u H e - 1 a P e ∆y(p E -p P ) + (1 -α u )[f e + u 0 E + (1 -f e + )u 0 P ] (3.37)
In this way, Rhie and Chow's special pressure-velocity coupling involves the addition of a higher order pressure gradient term. The restored linkage between the pressure differences across the control volume faces and the face velocities provides damping of the spurious oscillations in the solution due to the co-located arrangement.

Incompressible flow solver : SIMPLE algorithm

Due to the ease of their implementation and the lower peak memory requirements, segregated pressure-velocity coupling algorithms are commonly used for solving incompressible flow equations written in terms of the primitive variables velocity and pressure. The basic idea of SIMPLE-like algorithms is to develop an iterative procedure which constructs and solves a number of linear systems within each iteration. Both velocity and pressure fields are updated during each iteration so that the continuity equation is always satisfied, and velocity equations converge progressively to their final solution. The common characteristic of these algorithms is that a pressure-correction equation is constructed by combining discretized continuity equation and approximate forms of the momentum equations. The approximate pressure correction linear system is used afterwards to improve the current pressure estimate and the intermediate velocity solved from the momentum equation so as to force the modified velocity to satisfy the continuity condition for each control volume at each iteration level. The method is illustrated by considering the two-dimensional laminar steady flow equations in Cartesian co-ordinates.

In the first step of the SIMPLE calculation process, the so-called predictor step, a pressure field p * is guessed. Then, intermediate solutions u * and v * are obtained by solving discretized momentum equations of the x-momentum and y-momentum equations, respectively, using the guessed pressure field. The improved velocities can be expressed as:

u * P =   nb a nb u * nb + b u a P   P -∆y (p * e -p * w ) P (a P ) P (3.38) v * P =   nb a nb v * nb + b v a P   P -∆x (p * n -p * s ) P (a P ) P (3.39)
Following the predictor step, the calculated velocities, denoted as u * P and v * P , are based on the guessed pressure p * or from previous iterations. Now, let us introduce the correction p ′ as the difference between correct pressure field p and the guessed pressure field p * , so that

p = p * + p ′ (3.40)
Substitution of the correct pressure field p into the momentum equations yields the correct velocity field u. Since u * P and v * P do not satisfy continuity yet, similar procedure is applied to velocity components by introducing the corrections u ′ and v ′ as:

u = u * + u ′ (3.41) v = v * + v ′ (3.42)
The equations for the velocity correction terms, u ′ and v ′ , can be derived by subtraction of equations (3.38)

u ′ P =   nb a nb u ′ nb a P   P -∆y (p ′ e -p ′ w ) P (a P ) P (3.43) v ′ P =   nb a nb v ′ nb a P   P -∆x (p ′ n -p ′ s ) P (a P ) P (3.44)
The main approximation of the SIMPLE method lies in neglecting the terms

v ′ P = -∆x (p ′ n -p ′ s ) P (a P ) P (3.46)
Since a collocated grid is employed, a special interpolation scheme needs to be applied in order to avoid oscillating distribution of pressure. By adopting Rhie and chow interface interpolation, the cell face velocity corrections for faces e and n take the following form:

u ′ e = -∆y (p ′ E -p ′ P ) (a P ) e (3.47) v ′ n = -∆x (p ′ N -p ′ P ) (a P ) n (3.48)
Then the corrections are applied to face velocities to yield the following expressions:

u e = u * e -d e (p ′ E -p ′ P ) (3.49)
v n = v * n -d n (p ′ N -p ′ P ) (3.50)
where

d
a p p ′ P = nb a nb p ′ nb + b p (3.52)
where

a P = a W + a E + a S + a N b p = (ρu * A) w -(ρu * A) e + (ρv * A) s -(ρv * A) n
By applying Eq. (3.52) on all control volumes of the domain, the pressure corrections can be obtained at all nodes. These are used afterwards to correct the velocity field, which will now satisfy the continuity condition. However, since velocity components may not satisfy the momentum equations, another iteration must be performed using the solutions from the previous iteration as initial guesses and so forth with the iteration procedure until convergence is achieved.

In the solution procedure, the apparent viscosity η can be obtained from the guessed or intermediate velocity field and then be used to perform the current iteration. In this way, the discretization method of convective and diffusive terms described above can refer to both methods for Newtonian and non-Newtonian flows.

The sequence of operations in a solution procedure which employs the SIMPLE algorithm on co-located grid is summarized below:

1. An initial velocity field u 0 and v 0 is assumed 2. Apparent viscosity η is evaluated using regularization approach with Eq. (3. It is to be noted that, typically, all changes are under-relaxed to obtain numerical stability. Improved pressures p are obtained with:

p = p * + α p p ′ (3.53)
where α p is the pressure under-relaxation factor. The iteratively improved velocity components are also under-relaxed and would be computed as:

u = α u u (m) + (1 -α u )u (m-1) (3.54) and , v = α v v (m) + (1 -α v )v (m-1) (3.55)
where α u and α v represent velocity under-relaxation factors, u (m) and v (m) are corrected velocity components in the current iteration, and u (m-1) and v (m-1) represent their values determined by the previous iteration. Adequate values of under-relaxation factors α are essential to stabilize the iterative procedure. Unfortunately, the optimum values of under-relaxation factors are flow dependent and should be determined on a case-by-case basis.

Multigrid procedure

Two-grid algorithm

Iterative techniques are attractive because of their low storage overheads, specially for the solution of large systems of equations resulting from refined grids. In a CFD simulation, the efficiency of a solution technique is measured in particular by the computational cost put into achieving the desired accuracy. It is common knowledge that the discretization error is drastically reduced with highly refined meshes. However, the convergence rate towards the exact numerical solution of iterative methods, such as the Jacobi and Gauss-Seidel, rapidly gets worse as the mesh spacing is reduced. The application of the multigrid idea results in an approximately linear increase of computing time with grid refinement, allowing much finer grids to be used and therefore more accurate solutions to be obtained.

In what follows, we provide an overview of FAS Multigrid concept applied to a two-grid system. A two grid algorithm consists of performing a few smoothing iterations on the fine grid, approximation of the required correction on the coarse grid, prolongation of the coarse grid correction to the fine grid and again smoothing on the fine grid. Consider a linear system of equations, whose exact solution for any variable on grid h, ϕ h , satisfies the following equation :

A h ϕ h = S h (3.56)
Here, ϕ is the exact solution at all CV centers of the grid with spacing h, A is the coefficient matrix and S is the source term. After few SIMPLE iterations, the unconverged equation, having residual R h , can be expressed as:

A * h ϕ * h = S * h + R h (3.57)
where A * h and S * h are evaluated using the approximate solution ϕ * h . Subtracting (3.56) from (3.57) yields

A h ϕ h = S h + A * h ϕ * h -S * h -R h (3.58)
The obtained equation is used as the basis for multigrid coupling. The algebraic system (3.58) is then restricted on a coarser grid of spacing 2h, as follows:

A 2h ϕ 2h = S 2h + A 2h [I 2h h ϕ h ] -S 2h -[I 2h h R h ] (3.59)
The restriction of variable values from fine h to coarse 2h grid has to be performed by interpolation denoted by the operator I 2h h as ϕ 2h = I 2h h ϕ h . In this way, the coarse grid equations are derived. Variables and operators on grid 2h based on the restricted approximate solution on the finer grid h are denoted by (∼), while those being modified in the course of iterations on the coarser grid 2h by (∧). The underlined terms are kept unchanged during the coarse grid iteration and can be considered as an extra source term in Eq. (3.59). If the residual is zero ( R = 0), the above equation is satisfied, i.e. ( ϕ = ϕ). For consistency reasons, A 2h and S 2h are evaluated on the coarse grid in the same way as A 2h and S 2h using the restricted values [I 2h h ϕ h ]. It is noted that pressure operator is linear and therefore, no restriction of pressure from fine to coarse grid is needed. At each coarse grid iteration, the pressure is initialized as zero and the corrected pressure is then transferred to the fine grid to adjust the fine grid solutions. In contrast, the values ϕ 2h will in the course of coarse grid iterations depart from their initial values

[I 2h h ϕ h ].
After certain number of SIMPLE iterations, the system (3.59) can be solved to obtain ϕ 2h . Then, the coarse grid corrections for the full approximation velocity variables can be computed as:

δϕ 2h = ϕ 2h -I 2h h ϕ h (3.60)
The corrections obtained on the coarsest grid are transferred back to the fine grid through prolongation process denoted by I h 2h δϕ 2h , yielding a better estimate of the fine grid solution:

ϕ new h = ϕ old h + α M G [I h 2h δϕ 2h ] = ϕ old h + α M G [I h 2h ( ϕ 2h -I 2h h ϕ h )] (3.61)
where α M G represents the under-relaxation factor for coarse grid corrections. Then, few smoothing iterations are performed on the fine grid to eliminate any high frequency error components introduced by the prolongation process. A two-grid correction scheme is carried out according to the following sequence of instructions : Given (3.59), we assemble the equations for u-momentum and v-momentum on the coarse grid, whose solutions are u , v , as follows:

a p u p = nb a nb u nb -∆y(p ′ e -p ′ w ) + S u + a p u p - nb a nb u p -S u -R u (3.62)
a p v p = nb a nb v nb -∆x(p ′ n -p ′ s ) + S v + a p v p - nb a nb v p -S v -R v (3.63)
The SIMPLE procedure applied to the coarse grid equations involves the following steps:

1. The primitive variables are initialized as:

u p = u p , v p = v p , p ′ = 0
v n =    nb a nb v nb + S v + f v a P    n -∆x (p ′ N -p ′ P ) (a P ) n (3.65)
where f u and f v denote the underlying terms in (3.62) and (3.63), respectively. Analogous expressions for the other face velocities u w and v s can be derived in a similar manner 4. The face velocities are substituted into the mass continuity equation to derive the coarse-grid pressure correction equation. The pressure correction equation is p ′′ , given by :

a p p ′′ P = nb a nb p ′′ nb + b p (3.66)
where In practice, it is necessary to perform the above two-grid cycle continuously until desired accuracy is achieved. The coarse grid problem (3.59) can be solved using an even coarser grid 4h and so on. This whole process of going down from grid h to some coarsest grid and then back again up to grid h constitutes a multigrid cycle. In general, a number of multigrid cycles will be required to reduce the residuals by a given amount. Different kinds of cycles have been suggested and used.

b p = ρ[( u w -u w ) -( u e -u e )]∆y + ρ[( v s -v s ) -( v n -v n )]∆x
The MG initialization method used in the present study and summarized in Algorithm 1 consists of constructing the desirable number of geometric grid levels using the procedure outlined above. To begin the process, the initial solution is restricted all the way down to the coarsest level. The FAS multigrid cycle is then applied until a given order of residual reduction is obtained or the maximum number of cycles is reached. The obtained solution is then prolonged to the finer grid where the solution of the problem is required, see Algorithm 2. The coarsest grid used in the cycles was the 40 × 40 grid. remarkable power would be missed. The lowest frequency errors need to be reduced on a very coarse grid (8 h spacing or higher ) to decay quickly. The two-grid v-cycle can be extended in a natural way to more grids, going down to coarser grids (2 h, 4 h, 8 h, . . . ). As the coarse grid iterations are much faster than fine grid iterations, efficiency may be improved by the decision to switch from one grid to another on the rate of convergence through the combination of V-W cycles. Therefore, the W-cycle, which stays on coarse grids longer, is generally superior to a V-cycle and is used in the present work. Let it be noted that the optimum choice of parameters remains problem-dependent, but their effect on performance is not as dramatic as for the single-grid method.

Convergence criteria

Once the solver iteration is completed and the primitive variables corrected, convergence may be monitored by either the relative change of a variable or the residual for each of the discretization equations. The details on convergence criteria can be found in [START_REF] Kim | A study on convergence criteria for a simple-based finite-volume algorithm[END_REF]. The discretized equation for general flow variable ϕ at control volume i can be expressed as follows:

(a P ϕ P ) i = nb a nb ϕ nb i + (b ϕ ) i (3.67)
The relative changes of the primitive variables over two consecutive iterations can be computed to check the following convergence criteria:

max(|(ϕ P ) (k) i -(ϕ P ) (k-1) i |) < ϵ (3.68)
where ϵ is a specified small number. The final solution will not satisfy Eq. (3.67) exactly at all cells in the mesh, but after k iterations there will be a difference between the left and right hand sides. The absolute value of this difference at mesh cell i is termed the local residual (R ϕ )

(k) i : (R ϕ ) (k) i = (a P ϕ P ) (k) i - nb a nb ϕ nb (k) i -(b ϕ ) (k) i (3.69)
To get an indication of the convergence behavior across the whole flow field, we define the global residual R (k) ϕ , which is the sum of the local residuals over all N control volumes within the computational domain.

After k iterations we have:

R (k) ϕ = N i=1 (R ϕ ) (k) i = N i=1 (a P ϕ P ) (k) i - nb a nb ϕ nb (k) i -(b ϕ ) (k) i (3.70)
Inspection of equation (3.70) shows that the magnitude of the global residual R ϕ decreases as we get closer to the final solution, since the size of the local residuals should decrease in a converging sequence. Thus, it would seem that R ϕ might be a satisfactory single number indicator of convergence. However, the global residual will be larger in simulations where the flow variable ϕ has a larger magnitude, so we would need to specify different truncation values for R ϕ . This can be resolved if we use a global residual that is scaled to take out the magnitude of ϕ. Thus, we define the normalized global residual R ϕ for flow variable ϕ after k iterations as follows:

R ϕ (k) = R (k) ϕ /F ϕ R (3.71)
The normalisation factor F ϕ R is a reference level of the residuals for flow variable ϕ. Two common normalization methods are given below:

F ϕ R = R (k 0 ) ϕ (3.72) F ϕ R = N i=1 |(a P ϕ P ) (k) i | (3.73)

Test cases

The lid-driven cavity flow is considered by many researchers when testing computational fluid dynamics codes due to the wide range of fluid flow phenomena that can be observed in this simple flow configuration. Thereupon, the effectiveness of the solution procedure developed in this work are first demonstrated for flow in a square cavity of side L, where the top boundary (lid) moves horizontally towards the right with a uniform velocity U , while the remaining sides are fixed. The moving lid in conjunction with the shear properties of the fluid lead to a recirculation of the flow in the cavity, see Fig. 3.4.

Next, we assess the method described above on a steady pipe flow test as shown in the schematic diagram of the domain in Figure 3.5. For this

Numerical Results

All the numerical steady-state solutions presented in this section are obtained on Cartesian grids consisting of square control volumes. The Papanastasiou model is employed for the evaluation of viscosity and the stress growth parameter (m) is fixed to 200, a sufficiently high value that allows the ideal Bingham behavior to be approximated with satisfactory accuracy. The approximation of the convection terms in the momentum equations is carried out using the Quadratic Upwind Interpolation for the Convective Kinematics (QUICK) scheme.

In the present study, Reynolds number (Re) and Bingham number (Bn) are taken as independent parameters. These two numbers are defined as follows:

Re = ρU L µ p (3.74) Bn = τ y L µ p U (3.75)
with zero yield stress (Bn = 0) in case of Newtonian fluid. The Bingham number standing for the ratio of yield stress to viscous force is a crucial dimensionless number to describe viscoplastic materials.

Verification of the Numerical Method: Lid Driven Cavity

We now turn our attention to validating the numerical results obtained using the Simple/regularization procedure of our code. To this end, we consider the test case of lid-cavity flow, which is the prototypical recirculation flow and has long been used as a standard problem for testing and evaluating Navier-Stokes solvers. Many authors provided high quality benchmark results for this particular problem, most famous of which are Ghia et al. [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF] who used a coupled strongly implicit multigrid (CSI-MG) in a uniform mesh of 257 × 257. They presented a second-order accurate results that have served as "The" result for Newtonian flows to compare against ever since.

All the results presented in this section are obtained on Cartesian grids consisting of 320 × 320 square control volumes. To validate the present algorithm, first, a Newtonian flow (Bn = 0) in the Re = 10 3 case is calculated. Figure 3.7aa,b show the evolution of the horizontal and vertical velocities along the mid-planes x = 0.5 and y = 0.5, with the corresponding numerical result of Ghia et al. [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF]. The present results are in good agreement with those numerical data. Now we consider the steady viscoplastic flow at various Bingham numbers (Bn) in the case of Re = 10 3 . For this particular flow, many of the previous studies (Mitsoulis et al. [START_REF] Mitsoulis | Flow of Bingham plastics in a liddriven square cavity[END_REF]; Vola et al. [START_REF] Vola | Laminar unsteady flows of Bingham fluids: a numerical strategy and some benchmark results[END_REF]; Prashant et al. [START_REF] Prashant | Direct simulations of spherical particle motion in Bingham liquids[END_REF]; Syrakos et al. [START_REF] Syrakos | Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method[END_REF]) formulated the main quantitative results in terms of the features of the flow vortices. To this end, we need to calculate the stream function on the whole domain through integration of the velocity field (u = ∂ψ/∂y, v = -∂ψ/∂x). The location of the primary vortex formed in lid-driven cavity flow is significantly influenced by both Bingham and Reynolds numbers. This shifting of location is easily noticeable, which allows us to further assess the qualitative aspect of the simulation results. The evolution of the center position of the principal vortex is compared to The evolution of the principal vortex location as a function of the yield stress drawn in Figure 3.7 shows that the primary vortex goes up toward the moving lid as Bn increases. Increasing the yield stress leads to restricting the flow region closer to the upper moving boundary due to the resistance to lid motion, which gets higher. In viscoplastic flows, rigid zones are located at the bottom of the cavity where the shear stresses are very low as the fluid in contact is motionless due to the no-slip boundary condition. These unyielded regions expand significantly with the increase of the Bingham number, leading the main vortex to move toward the upper boundary as there is less space for flow to happen. The streamline contour evolution shown in Figures , respectively. When Re = 10, the streamline contours are symmetric for all Bingham numbers, 2, 10, and 100. However, the unyielded zone is proportional to the yield stress magnitude. At Re = 10 3 , the primary vortex shifts slightly to the right for the case of Bn = 100, and moves further in the same direction for Bn = 10, then towards the center of the cavity for Bn = 2. We see that increasing the Reynolds number moves the vortex first toward the right, and then downwards and left to the center depending on the value of the Bingham number. This trend can be observed for higher Bingham numbers if the Reynolds number is large enough as greater stresses are required in order to make the material flow. In the works of [START_REF] Syrakos | Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method[END_REF][START_REF] Prashant | Direct simulations of spherical particle motion in Bingham liquids[END_REF], similar conclusions were drawn, the flow field at high Bn numbers follow the same trend of any lower Bn number if the Re number is also sufficiently decreased. It is clear that with an increasing exponent m, we can achieve quick stress growth at very low shear rates, which is consistent with the behavior of the material in its unyielded region. In fact, regularized models yield results that do not contain truly unyielded or solid zones and hence the solution may, in some cases, be sensitive to the chosen stress growth exponent. Previous studies, including [START_REF] Syrakos | Performance of the finite volume method in solving regularised Bingham flows: Inertia effects in the lid-driven cavity flow[END_REF], found that different values of 100, 200, and 400 for m do not influence the location of the yield surface significantly, for the particular case of cavity flow. However, very large values of the exponent parameter may be needed if the region where the stress close to the yield stress is relatively large.

Regarding the viscoplastic flows though pipes, two distinct unyielded regions can possibly occur. First, in the layer adjacent to no-slip boundary where the velocity of the fluid is set to zero, it is called the stagnant zone. Then, a second unyielded region with non-zero fluid velocity is referred to as the plug zone, see Figure 3.11. One result of importance in a range of applications is the location of the plug zone, where γ tends to zero, as previously stated. In Figure 3.12, we present a comparison of the velocity profiles for different values of m: 100, 200, 400, and 800. As m is increased, the location of the unyielded region would be slightly affected.

The difference especially vanished between m = 400 and m = 800 for both cases. Figure 3.12 shows velocity profiles, which deviate slightly from the reference result considered here in red dotted lines (m = 800) and plug zones which get flatter as m is increased. Indeed, beyond a certain value of stress growth exponent, yield lines (τ = τ y ) are less sensitive to its variation. This verifies the assumption that the yield areas calculated using a regularization approach will converge to the true yield areas as m → ∞. It should be pointed out that, unlike the velocity field, there is no theoretical proof that yield surfaces converge to the exact solution [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF]. In Figure 3.13, we compare the simulated results of the apparent viscosity profiles along the circular cross-section of the pipe calculated with m = 100 and m = 800. It is shown that using a high value of the exponent m leads to fast growing apparent viscosity as the shear rates are decreasing, meaning that the larger the value of m, the larger the region of the solid-like behavior of the Bingham model is reproduced. Figure 3.13 also shows that, with regularized models, the unyielded material (plug zone) is no longer a rigid solid but a highly viscous fluid that approximates the ideal viscoplastic behavior. The degree of approximation depends on the adequate choice of the stress growth parameter in order to satisfy the von Mises criterion. An important issue when using the Papanastasiou regularization is that the degree of non linearity increass sharply with higher values of m that ensure better approximations. Therefore, as a matter of practicality, the choice of exponent m has to be reasonably reduced depending on the problem under study. In this work, we limited our study to the comparisons showing the influence of m on the quality of the results. The reader may refer to the interesting assessments provided by Syrakos et al. [START_REF] Syrakos | Performance of the finite volume method in solving regularised Bingham flows: Inertia effects in the lid-driven cavity flow[END_REF] and Frigaard et al. [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF] on the issues of feasibility and the computational effort required as a function of the stress growth exponent. They suggested that, in some cases, to achieve a fully-converged solution with a high value of m, one can start with a very low value of m as an initial guess, and progressively increase it every certain number of SIMPLE iterations.

Algebraic Convergence of the SIMPLE/Regularization Procedure

In viscoplastic flows, derivatives of the flow variables take significantly higher values across the yield surfaces in the limit between the unyielded and yielded regions. This disparity is due to the fact that these derivatives are zero inside the rigid zones and non-zero outside. In addition, increasing the Bingham number and growth stress parameter would further widen this gap. As a result, solution errors are found to be much larger there than the rest of the domain, causing inefficiencies and deterioration in convergence rates of iterative solvers. In such cases, using fine grids can considerably reduce the numerical errors and address more efficiently the high values of the derivatives. However, it is well known that the number of iterations required for convergence increases almost linearly with the number of control volumes for the single grid set [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF][START_REF] Peric | A finite volume multigrid method for calculating turbulent flows[END_REF]. To address this limitation, one may notice that using an appropriate initial solution is very important for fine grids and contributes to the enhancement of the convergence of iterative solvers. To this end, a steady lid driven cavity flow simulation is carried out for different Bingham numbers 0, 1, and 10 at Re = 10 2 using single-grid (SG) and multi-grid (MG) initialization methods.

Figure 3.14 demonstrates the convergence of the L ∞ norm of u and v residuals as a function of the number of SIMPLE iterations on the finer grids. It shows that the performance of SIMPLE, as a single-grid solver, deteriorates, and more iterations are needed as the number of control volumes increases for all cases. In contrast, computations with MG initialization exhibit near constant number for both 160 × 160 and 320 × 320 grids. In all cases, the number of iterations has drastically been reduced compared to the single grid computation. One may also observe that convergence becomes faster in the case of Bn = 10, but it also becomes more unstable and oscillatory. It should be noted that for higher Bingham numbers, convergence difficulties may be encountered. In order to overcome this difficulty, low under-relaxation factors for pressure and velocity field, denoted by α p , α u , and α v can be considered within the numerical solver [START_REF] Syrakos | Performance of the finite volume method in solving regularised Bingham flows: Inertia effects in the lid-driven cavity flow[END_REF]. This technique was found to be efficient in making the solver more robust, but we observed that it leads to a slow down in convergence rate. Another remedy to this problem could be the suggestion of Ferziger and Peric [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF] in a multigrid context. They recommend updating the viscosity on the finer grid only and holding its value constant on coarse grids within the multigrid cycle. The performance of multi-grid (MG) initialization is based on the speedup characteristics, and a comparison between the number of iterations and computation times, as summarized in Tables 3.1 and 3.2. In all cases, CPU time is recorded until the convergence of residuals is achieved. In the case of MG Initialization, the total simulation time includes computation time on both multilevel grids during initialization, and the finer grid, which is the true representation of the simulation time. The simulations are performed on a Core i7-6700HQ CPU @ 2.60GHz × 8, 8 GB RAM-based computer. According to our numerical experiments, the speed-up ratio between single grid and multigrid initialization methods is around 6 for 320 × 320 grids and up to 3 for 160 × 160 grids. This confirms the potential of multigrid initialization technique for improving the convergence over single-grid method by significantly reducing CPU time and number of iterations. Since MG initialization does most of the computations on coarse grids, this initial-ization procedure is computationally inexpensive and, for very dense grids, a good initial solution can be obtained in a fraction of the time spent to converge on a final solution. The computational effort can be further reduced when parallel algorithms and high performance computing machines are used [START_REF] Roy | A parallel multigrid finitevolume solver on a collocated grid for incompressible Navier-Stokes equations[END_REF].

Conclusions

This chapter introduced the idea of multigrid initialization to overcome the drawbacks of solvers using single grids for viscoplastic flows. The proposed solver is easy to implement and significantly improves the convergence rate by reducing CPU time and number of iterations. The computational procedure employed in this work is based on the finite volume / SIMPLE algorithm in conjunction with the regularization scheme of Papanastasiou and is applied for the solution of the lid-driven cavity and pipe flow problems. The convection terms are discretized using the QUICK scheme in order to avoid any artificial diffusion that may be caused by interpolation schemes of low order of accuracy. The capability of the algorithm for solving viscoplastic flows was verified against a number of benchmark solutions for a range of Bingham numbers in the case of lid-driven cavity flow. A qualitative assessment was also conducted by monitoring the positioning of the principal vortex with regard to Bingham and Reynolds numbers. Results from the current study were accurate and compare favorably with the published results. In order to analyze the influence of the growth stress parameter m of the regularization model on fluid motion, pipe flow tests for different values of m and Bn were considered. Numerical simulations confirmed that a large value of m reproduced a large threshold-like region. The result also showed that as m increases, the yield surface converges to the true yield region, with no significant influence beyond a value of 400 in this particular case. Furthermore, the results clearly reveal that introducing the multigrid method as an initialization procedure yields better convergence for Navier-Stokes equations for viscoplastic flows compared to the single-grid method. A multi-grid initialization speed-up ratio as high as 6.45 was achieved for the finest grid size (320 × 320).

The present study can be extended to the numerical solution of non-Newtonian problems with a different definition of effective viscosity and accounting for complex geometries encountered in many industrial processes. This part is the subject of an ongoing investigation that is to be published in our future works.

Chapter 4

Rheological Characterization of Phosphate Slurry

4.1 The assessment of the rheological characteristics of Phosphate Slurry

Introduction

Particles in a liquid matrix can be considered as obstacles, hindering the liquid's flow and, therefore, increasing the flow resistance, i.e., the viscosity. The mechanics of dilute and semi-dilute suspensions representing very low particle concentrations (volume fraction ϕ < 25%) are well understood. Simple models for describing the increased viscosity at a given shear rate (or apparent viscosity) have been proposed by Einstein [START_REF] Einstein | Eine Neue Bestimmung Der Moleküldimensionen[END_REF] and Batchelor [START_REF] Batchelor | The Effect of Brownian Motion on the Bulk Stress in a Suspension of Spherical Particles[END_REF]. Their formulas were adapted to monodisperse particles for maximum solid volume fractions of 0.15 to 0.25 and link the increase in bulk viscosity of the suspension η to the viscosity of the liquid η 0 , as well as the solid fraction in the suspension. However, linearity is restricted to low solid fractions ϕ. For dense or highly concentrated suspensions, the apparent relative viscosity (η r = η/η 0 ) increases rapidly and without bound with ϕ. At a large solid fraction, small changes in concentration lead to large variations in viscosity. If nonhydrodynamic interparticle forces are absent, suspensions with a high solid fraction can be described with the satisfactory model of Krieger and Dougherty [START_REF] Krieger | A Mechanism for Non-Newtonian Flow in Suspensions of Rigid Spheres[END_REF], which is valid for low and high shear rates. This model requires an additional input parameter, defined as the maximum packing fraction of particles ϕ m , which is related to the particles' properties. Physically, ϕ m refers to the maximum solid fraction possible in suspension for a given particle size distribution and packing arrangement [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF][START_REF] Ouchiyama | Porosity of a Mass of Solid Particles Having a Range of Sizes[END_REF]. Many other functions have been proposed and provide excellent fits to experimental data [START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF][START_REF] Mooney | The Viscosity of a Concentrated Suspension of Spherical Particles[END_REF][START_REF] Frankel | On the Viscosity of a Concentrated Suspension of Solid Spheres[END_REF].

Dense particulate suspensions often display a variety of non-Newtonian flow properties, including shear thinning, shear thickening, and reversible, shear-induced solidification called shear jamming. Recently, it has been shown that interparticle friction drastically influences this behavior. Shear jamming (SJ) is a phenomenon where a suspension is fluid-like at low stresses and jams into a solid-like state at high shear. This shear-induced solidification occurs at particle concentrations below the traditional, frictionless jamming packing fraction [START_REF] Barnes | Dilatancy") in Suspensions of Nonaggregating Solid Particles Dispersed in Newtonian Liquids[END_REF][START_REF] Bender | Reversible Shear Thickening in Monodisperse and Bidisperse Colloidal Dispersions[END_REF][START_REF] Maranzano | The Effects of Interparticle Interactions and Particle Size on Reversible Shear Thickening: Hard-Sphere Colloidal Dispersions[END_REF][START_REF] Brown | The Role of Dilation and Confining Stresses in Shear Thickening of Dense Suspensions[END_REF][START_REF] Wagner | Shear Thickening in Colloidal Dispersions[END_REF].

With an increased solid volume fraction ϕ, the apparent viscosity η increases but this increase is not the same for different maximum packing fractions ϕ m . In fact, the rate of the viscosity increase depends on the ratio ϕ/ϕ m . The most important factors affecting the maximum packing fraction ϕ m are the shape and size distribution of the particles. Since smaller particles can fit into the voids created by larger ones, the maximum packing fraction of polydispersed particles is significantly higher than in monodispersed particle systems [START_REF] Mcgeary | Mechanical Packing of Spherical Particles[END_REF][START_REF] Ouchiyama | Porosity of a Mass of Solid Particles Having a Range of Sizes[END_REF][START_REF] Zou | Packing of Multi-Sized Mixtures of Wet Coarse Spheres[END_REF]. Many experiments with bidisperse and polydisperse suspensions have been conducted to establish a relationship between ϕ m and the suspension viscosity [START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF][START_REF] Kelessidis | Optimal Determination of Rheological Parameters for Herschel-Bulkley Drilling Fluids and Impact on Pressure Drop, Velocity Profiles and Penetration Rates during Drilling[END_REF][START_REF] Farris | Prediction of the Viscosity of Multimodal Suspensions from Unimodal Viscosity Data[END_REF][START_REF] Chang | Effect of Particle Size Distributions on the Rheology of Concentrated Bimodal Suspensions[END_REF]. In general, these studies showed that the apparent viscosity η tends to decrease when increasing the maximum packing fraction ϕ m at a constant volume fraction ϕ. This indicates that the concentration of solids in a solid-liquid mixture can be raised markedly while maintaining a desired apparent viscosity.

At a constant solid fraction, a reduction in particle size leads to an increase in apparent viscosity. With smaller particle sizes, the total surface area becomes larger, which binds to water molecules and results in a higher number of particles in a given volume. When the volume fraction increases, interactions between particles become significant, particularly for very small particles, thus substantially raising the applied force required to shear the sample. For larger particles that are fewer in number, the effect of particleparticle interactions is considerably reduced [START_REF] Austin | Process Engineering of Size Reduction: Ball Milling[END_REF][START_REF]Thixotropy. In Colloidal Suspension Rheology[END_REF]. However, for most real suspensions used in industry and mineral processing, the particle size is not uniform and is more conveniently described as a particle size distribution. Under small shear stress, Brownian motion dominates the behavior of the finest particles, while for coarse particles, frictional contacts or hydrodynamic forces tend to dominate [START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF][START_REF] Denn | Rheology of Non-Brownian Suspensions[END_REF]. In this study, we seek to investigate the rheological behavior of phosphate suspension with the aim of optimizing their dispersity and stability throughout future research works. Regarding this concern, a few studies on phosphate composites have been conducted for diverse applications. An early study by Bujake [START_REF] Bujake | Rheology of Concentrated Dicalcium Phosphate Suspensions[END_REF] investigated the effect of additives on the rheological and sedimentation behavior of a phosphate slurry. He suggested that particle-particle interactions play an important role in the flow of these systems. Benretem et al. [START_REF] Benretem | Factors Influencing Slurry Rheology[END_REF] stated that the rheological characteristics of water-phosphate slurries strongly depend on the solids concentration and increase markedly with decreasing particle size. In a comparable study, Belbsir et al. [START_REF] Hamza | Study of the Rheological Behavior of Phosphate Slurry and Its Derivatives Products[END_REF] drew similar conclusions and showed that viscosity tends to decrease with increasing temperatures and can also be reduced using appropriate additives. In a more recent study, Shaikh [START_REF] Shaikh | Rheological, and Mechanical Properties of Polypropylene/Phosphate Ore Composites[END_REF] reported the fabrication of a geopolymer of phosphate ore and polypropylene. He found that the composites of smaller phosphate ore particles displayed a highly complex viscosity. Since the transported slurry has a continuous particle size distribution, it is important to discriminate the role of this distribution in the process. Thus, in this work, we put our attention toward the influence of the D-value ratio on apparent viscosity, as well as the solid concentrations. The following sections present some basic information about the experimental work and the methodology carried out to provide a description of the non-Newtonian behavior that is typically exhibited by concentrated phosphate slurries. The main objective of this study was to evaluate the predictive or fitting capabilities of commonly used yield-pseudoplastic models for the rheological data description of concentrated phosphate slurries. The rheological parameters for the Bingham, Casson, Herschel-Bulkley, and Robertson-Stiff models were computed for six samples. A variety of samples with phosphate ore concentrations ranging from 51 to 56 wt.% in water, i.e., 29 to 33% v /v, were employed for this rheological study. These slurries with relatively coarse particles are a representation of what is typically found in phosphate rock processing. The rheological models were used to obtain analytical expressions relating the mean flow velocity to the pressure drop and to determine velocity profiles in the laminar regime.

Materials and Methods

Concentrated Phosphate Slurry Samples

Phosphate rock is mined at different locations and then processed to be stored in agitated tanks at the main slurry pipeline head station. The solid particles in the present study were provided by OCP (Khouribga, Morocco). The samples prepared at our laboratory were a mixture containing rigid randomly shaped solid particles with an equivalent spherical diameter ranging from 1 to 500 µm. The particle density was 2.4t/m 3 when the pores were filled with air and 2.6t/m 3 when the pores were filled with water and stabilized by a specific soluble polymer, such as PEG. Various grades and phosphate ore concentrations of 51, 54, and 56 wt.% in water were used for this rheological study. These samples, labeled S1 through S6, differed mainly in their solids concentrations and particle size distributions.

The main ingredients of the concentrated phosphate slurry (CPS) sam-ples were CaO and P 2 O 5 , even though they also contained other elements, such as SiO 2 , CO 2 , Fe, MgO, SO 3 , and Al 2 O 3 . The phosphate content of phosphorite P 2 O 5 is often expressed as a percentage of bone phosphate of lime (BPL) [START_REF] Fertilizer | United Nations Industrial Development Organization (UNIDO) and International Fertilizer Development Center (IFDC)[END_REF]. Table . 4.1 summarizes the elemental concentrations of each sample. The samples have a similar composition that can be considered a representation of what is specifically mined in Morocco. Phosphorite deposit compositions may differ from one location to another. The analysis of the particle size distribution provides essential data for forecasting the slurry behavior in the transport process or its physical properties. Measurements were carried out for each sample using a Mastersizer 3000 granulometer from Malvern (Worcestershire, UK). The analysis of the particles' size distribution allowed us to define the percentage of the total volume in each sample containing particles with diameters less than the following values: 500 µm, 400 µm, 315 µm, 250 µm, 150 µm, 74 µm, 44 µm, 40 µm, and 37 µm. The most commonly used metrics when describing particle size distributions are the D-values: D 20 , D 50 , and D 80 , which are the intercepts for, respectively, 20%, 50%, and 80% of the cumulative mass. Although a comparison of D 50 and D 80 values indicated little difference between the samples, the D 20 values were dissimilar. The particle sizes in the CPS samples ranged from around 1 µm to 500 µm, where the cumulative particle size distribution of all CPS samples is shown in Fig. 4.1.

Suspension Characteristics

The suspended particles in the present study were assumed to be noncolloidal, which means that Brownian motion in the suspensions could be neglected. The interaction of the forces describing the rheophysical properties of the suspensions can be defined using nondimensional numbers that were first defined by Krieger (1959). In particular, the impact of Brownian motion is governed by the Peclet number, Pe, which is defined as the ratio of the hydrodynamic force and the Brownian force [START_REF] Stickel | Fluid Mechanics And Rheology Of Dense Suspensions[END_REF][START_REF] Frankel | On the Viscosity of a Concentrated Suspension of Solid Spheres[END_REF]:

P e = 6πη 0 a 3 . γ kT , (1)
where η 0 is the suspending fluid viscosity, a is the average radius of particles,

.

γ is the shear rate, k = 1.38 × 10 -23 J K -1 is the Boltzmann constant, and T is the absolute temperature. For Pe > 10 3 , hydrodynamic forces dominate [START_REF] Stickel | Fluid Mechanics And Rheology Of Dense Suspensions[END_REF][START_REF] Frankel | On the Viscosity of a Concentrated Suspension of Solid Spheres[END_REF]. The relative importance of inertial forces at the particle scale is governed by the particles' Reynolds number Re p , which is given by:

Re p = ρ 0 a 2 . γ µ 0 , (2)
where ρ 0 is the density of the suspending liquid. For Re p > 10 -3 , inertial forces cannot be neglected [START_REF] Stickel | Fluid Mechanics And Rheology Of Dense Suspensions[END_REF][START_REF] Mueller | The Rheology of Suspensions of Solid Particles[END_REF]. Thus, for low shear rates in this work, the suspended particles were strongly influenced by inertial forces.

Experimental Data

The shear rheology of the slurry was measured using a RheolabQC rotational rheometer from Anton Paar (Ostfildern, Germany) with a concentric cylinder geometry, where the inner and outer cylinders represented, respectively, the test bob and cup (see Fig . 4.2). This device uses the Searle principle with a rotating inner cylinder and stationary outer cylinder and can measure the dynamic viscosity of low-viscosity samples up to semi-solid specimens [START_REF]Plastics-Polymers/Resins in the Liquid State or as Emulsions or Dispersions-Determination of Viscosity Using a Rotational Viscometer with Defined Shear Rate[END_REF][START_REF] Macosko | Rheology: Principles, Measurements, and Applications[END_REF]. The Searle system is the most commonly used system for mineral pulps. Coaxial (or concentric) cylinder systems are absolute measuring systems that conform to "Deutsches Institut für Normung" (DIN) and the international standards organization (ISO) standards. Measuring the rotational speed Ω of the electronically communicated (EC) motor and produced torque M allowed us to calculate the shear stress and the shear rate according to the following relationships (standard: ISO 3219 [START_REF] Van Wazer | Viscosity and Flow Measurement. A Laboratory Handbook of Rheology[END_REF] (δ ≤ 1.2)):

.

γ = 1 + δ 2 δ 2 -1 Ω, (3)
τ = 1 + δ 2 2δ 2 M 2πLR 2 i C L , (4)
where δ = R e /R i is the ratio of the outer cylinder and inner cylinder radii (R e and R i , respectively), L is the bob length, and C L = 1.28 is an end effect correction factor. The rheological behavior of the CPS samples was determined by measuring the shear stress at various applied shear rates. The shear stresses were measured at each point while decreasing the shear rate from 1000 s -1 to 112 s -1 using a coaxial geometry with a rotating cylinder and a cylindrical cup, with diameters of 38.713 mm and 44 mm, respectively. The estimated shear rate range occurring during the transport process was 200 to 600 s -1 , which is well within the shear rate range tested. The fixed distance between the two measuring cylinders was 2.644 mm. The experiments were conducted at controlled shear rates in the same range as those expected in pipe flows. The experimental data of the CPS samples were obtained at room temperature. For each measurement, the suspensions were vigorously mixed in a container beforehand. A pre-shear is required to erase material memory and to have similar initial conditions for all samples. Then, after prolonged mixing (>30 s), the necessary volume was quickly poured into the outer cylinder. The inner cylinder was subsequently immersed into the slurry. The stress τ exerted on this cylinder was recorded as a function of the shear rate .

γ. The rotational tests were carried out in a preset shear rate ramp mode, descending in steps, while recording the shear stress for each set shear rate. In order to allow the sample to adapt itself to each shear rate step, a duration of three seconds was maintained for each measuring point. It is important to note that the chosen Couette geometry (Fig. 4.2) with the adapted gap allowed us to assume that we would not observe common sudden changes in the rheogram plots at high shear rates and in the ranges of the yield-stress calculations. We intended for the applied pre-shear and the preset shear rate ramp mode to allow for homogenized suspensions and, hence, we could assume their homogeneity. Furthermore, we had other work in progress using a helical geometry and our calculated yield stress seemed to corroborate with those obtained using our adapted geometry and procedure. The density of the solid particles was higher than that of water, and the settling experiments performed on a mixture of water with 10 vol.% of particles with sizes between 100 and 500 m revealed small settling velocities. As a result, the effects of the particle settling on the sample concentration were negligible during the performed experiments (less than one minute).

Rheological Models

Phosphate ore slurry is a suspension of insoluble particles in a continuous phase of water. This suspension displays non-Newtonian flow behavior and exhibits yield stress as the shear rate tends toward zero. These materials only flow due to a certain constraint called the yield stress τ y . Below this stress, they behave like solids, while above this stress, they may flow like a Newtonian fluid or exhibit shear-thinning. Knowledge about the yield stress is essential for transporting suspensions, especially for resuspending particles after they have settled in a pipeline or channel. This property is found in certain polymer solutions and very concentrated suspensions, and in everyday life with toothpaste, yogurt, ketchup, clays (laponite and bentonite, for example), and drilling muds. This behavior can be interpreted as a presence of aggregates in the material that remains up to the yield stress. However, some authors find the idea of yield stress controversial [START_REF] Barnes | The "yield Stress Myth?" Paper-21 Years On[END_REF] and doubt whether true yield stress actually exists. The main reason for this debate and discussion in the literature is that no equipment has so far allowed researchers to measure the shear stress at very low shear rates without being affected by wall-slip or end effects. For further reading, we refer to the work by Coussot et al. [START_REF] Coussot | Macroscopic vs Local Rheology of Yield Stress Fluids[END_REF].

For Newtonian fluids, most pure liquids, and gases, shear stress is a linear relation with respect to the shear rate:

τ = µ . γ, (5)
where the constant of proportionality µ is known as the viscosity of the fluid. Non-Newtonian materials, suspensions of solids, polymer solutions, and polymer melts exhibit a different response, for which an apparent viscosity η is often assigned:

τ = η( . γ) . γ. ( 6
)
The simplest representation of such complex fluids is the Bingham model, which gives the following relationship between the shear stress and shear rate:

τ = τ yB + η B . γ, (7)
where τ yB is the yield stress and η B is the constant plastic viscosity. This model has the advantage of giving a linear relationship after an initial yield and defining the minimum stress value at which the material starts to flow, which is also referred to as the Bingham yield point τ B . However, recent studies show that this model fails to conform to complex fluid behavior at low shear rates [START_REF] Thota Radhakrishnan | Rheological Characterisation of Concentrated Domestic Slurry[END_REF] and may result in unrealistically high values of τ yB [START_REF] Khataniar | Evaluation of Rheological Models and Application to Flow Regime Determination[END_REF].

A general form of the Bingham model, proposed by Herschel and Bulkley [START_REF] Herschel | Consistency of Rubber Benzene Solutions[END_REF], is used to describe the flow of pseudoplastic fluids, which require yield stress to initiate the flow and then exhibit a nonlinear flow curve. The Herschel-Bulkley model can describe the rheological behavior of mineral slurries [START_REF] Turian | Characterization, Settling, and Rheology of Concentrated Fine Particulate Mineral Slurries[END_REF][START_REF] He | Slurry Rheology in Wet Ultrafine Grinding of Industrial Minerals: A Review[END_REF]. Other investigations have shown that this model can also describe sludge behavior and is most commonly used to characterize concentrated sludge [START_REF] Einstein | Eine Neue Bestimmung Der Moleküldimensionen[END_REF][START_REF] Baudez | Rheology of Aging, Concentrated, Polymeric Suspensions: Application to Pasty Sewage Sludges[END_REF], as well as most drilling fluids or cement slurries [START_REF] Núñez-Santiago | Rheological Evaluation of Non-Newtonian Mexican Nixtamalised Maize and Dry Processed Masa Flours[END_REF][START_REF] Bailey | Investigation of Methods for Direct Rheological Model Parameter Estimation[END_REF]. The Herschel-Bulkley model can be used with the constitutive equation:

τ = τ yH + K . γ n , (8)
where the yield stress τ yH in this equation represents the intercept of the line fitting the yield stress region, K is the consistency index, and n is the flow index. With the use of the third parameter, this model often provides a better fit to experimental data.

Another commonly used rheological model for yield-pseudoplastic fluids is the so-called Casson model, which was initially generated to describe ink and paint behavior. However, it has demonstrated a good approximation for non-Newtonian mineral suspensions [START_REF] Turian | Characterization, Settling, and Rheology of Concentrated Fine Particulate Mineral Slurries[END_REF][START_REF] He | Slurry Rheology in Wet Ultrafine Grinding of Industrial Minerals: A Review[END_REF] and complex fluids in the petroleum industry [START_REF] Núñez-Santiago | Rheological Evaluation of Non-Newtonian Mexican Nixtamalised Maize and Dry Processed Masa Flours[END_REF][START_REF] Bailey | Investigation of Methods for Direct Rheological Model Parameter Estimation[END_REF]. The Casson model is given by the following equation:

√ τ = √ τ yc + k c . γ. (9)
The Robertson-Stiff model is a three-parameter model that was developed in 1976 by Robertson and Stiff to describe the rheology of drilling fluids in the petroleum industry and cement slurries, though it can also be used for fluids as diverse as maize flour pastes [START_REF] Batchelor | The Effect of Brownian Motion on the Bulk Stress in a Suspension of Spherical Particles[END_REF] and gels. It is expressed as:

τ = K( . γ + . γ 0 ) n . (10)
The parameters K and n can be considered similar to the parameters (K and n) of the Herschel-Bulkley model. However, the parameter γ 0 appears as a correction to the shear rate rather than the shear stress and the yield stress is defined as:

τ yRS = K . γ n 0 . (11)
For Equations ( 7)-( 9) and ( 11), the following condition can be applied:

. γ = 0, |τ | < |τ y | . (12)

Determination of Model Parameters

When comparing the measuring data, it is not useful to compare all values of one test with those of another mostly due to the large number of individual measuring points. Mathematical model functions for curve fitting are therefore used to characterize complete flow or viscosity curves, resulting in only a small number of curve parameters. This simplifies the comparison of measuring curves since there are only a few model parameters left to compare. Fitting is also called approximation and the corresponding functions are often referred to as regression models. Not every model function can be used for each kind of flow behavior. If the correlation value (e.g., in %) indicates insufficient agreement between the measuring data and model function, it is useful to try another model function. It is also important to keep in mind that both model-specific coefficients and exponents are purely mathematical variables and do not represent real measuring data in principle.

Here, with a controlled shear rate (CSR), the rotational speeds (or shear rates) are preset in the form of steps. However, when using this kind of testing, a yield point cannot be determined directly. It is instead calculated using a fitting function, which is adapted to the available measuring points of the flow curve. Curve fitting is carried out using one of the various model functions. For all approximation models, the dynamic yield stress value τ y is determined via extrapolation of the flow curve toward the shear rate value . γ = 0 or at the point of intersection of the fitting function and the τ -axis. The obtained value of the τ y method strongly depends on the chosen rheological model and the accuracy of the measurements.

To determine the rheological coefficients of each model, curve fitting is performed. Curved relationships between variables are not as straightforward to fit and interpret as linear relationships, but there are standard ways for dealing with this problem [START_REF] Núñez-Santiago | Rheological Evaluation of Non-Newtonian Mexican Nixtamalised Maize and Dry Processed Masa Flours[END_REF][START_REF] Ohen | Golden Section Search Method for Determining Parameters in Robertson-Stiff Non-Newtonian Fluid Model[END_REF]. In total, six sets of rheometer data for the shear stress τ i and shear rate . γ were analyzed. Three statistical indicators were computed and reported: the correlation coefficient R 2 , the sum of square errors SSE, and the root mean square error RMSE. The SSE was calculated according to:

SSE = m i=0 (τ i -τi ) 2 , ( 13
)
where τi is the predicted or modeled shear stress value. The root-meansquare error is given as:

RMSE = m i=0 (τ i -τi ) 2 N m , (14)
where N m is the number of measurements.

Results and Discussions

Rheological and Statistical Evaluation

For each test, the Anton Paar RheolabQC rotational rheometer provided the variation of the shear stress with the shear rate. Figs. 4.4, 4.5 and 4.6 show the rheograms obtained for the six samples at solids concentrations ranging from 51 to 58 wt.%. In all cases, the flow curves appeared to exhibit non-Newtonian behavior with a nonzero intercept on the stress axis and could be described as yield pseudoplastic. It was also observed that the shear-thinning was accentuated and the dynamic yield stress increased as the solids concentrations were increased for all cases. In this study, a statistical model fitting analysis was performed on all five models previously discussed to determine the model that was best able to predict the rheological behavior of the CPS samples. 4.2 lists the measured shear stresses of the fluids at different rotational speeds. The raw data obtained were analyzed using several models to determine the various rheological parameters, such as the yield value, consistency, and shear-thinning index. For discussion purposes and in order to limit the size of this chapter, only six samples of the total number of experimental runs are presented. The plots of experimental data and shear stress predicted by each model versus the shear rate for all samples are shown in Fig. 4.4, 4.5 and 4.6. 4.3 and 4.4 summarize, respectively, the rheological parameters and statistical indicators obtained after the fitting processes. It can be seen that all the models fit the general trend well but the Bingham model failed to match the experimental stress as closely as the other models. All the models had high R 2 values ranging from 0.90 to 0.99 and the highest value was from the Robertson-Stiff fit. The same procedure was applied to the data from all the tests and similar results were obtained; the Robertson-Stiff fit model appeared to fit the data better than the other models tested. The plots show that increasing viscosity or decreasing size ratios resulted in shear-thinning (curved line) for intermediate shear rates. The data fit to the Bingham and Casson equations resulted in better correlation coefficients for samples with weak shear-thinning: S1, S2, S3, and S4. Thus, for the diluted phosphate slurry, the Bingham and Casson models gave a fair description of the rheological behavior.

However, the Robertson-Stiff and Herschel-Bulkley models provided a better fit for all samples, with the correlation coefficients R 2 ranging from 0.97 to 0.99, thus reflecting the accuracy of these models. The sum of the square errors was also lower in comparison with the other models and ranged from 0.48 to 93.80. These indicators were significantly lower than for the Bingham and Casson models, perhaps indicating that three-parameter models were more adequate to describe the rheological behavior of phosphate slurries. Furthermore, the shear-thinning behavior was satisfactorily captured by these models.

Effect of Particles Concentration

The rheology of a particle suspension is a complex function of its physical properties and of processes that occur at the scale of the dispersed particles. The most important factors are the particle volume fraction ϕ, the particle shape, the interactions between the particles, and the spatial arrangement Clearly, at a solids concentration of 51 wt.% (29% by volume), the behavior appeared to involve weak shear-thinning at shear rates less than 300 s -1 that tended to be considerably reduced at shear rates larger than 300 s -1 . At even higher particle concentrations, the flow curve became more convex upward, indicating shear-thinning in a wider range of shear rates (<400 s -1 ). This was due to large interparticle distances between the solid particles in a dilute slurry such that attractive forces existing between them had no effect. A larger concentration fraction of particles implied a higher concentration of particle-particle bonds wherefore greater shear energy was required to disrupt them. The higher energy dissipation was reflected by a higher viscosity and yield stress.

Effect of Particle Size Distribution

The particle size distribution is an important factor that needs to be taken into account to understand the rheology of certain industrial suspensions. During mineral processing, phosphate rock is subjected to various grinding operations that result in a wide particle size range, from nearly one to few hundreds of micrometers. Thus, to better understand the effect of particle size distribution changes on the viscosity of the CPS samples and to keep the size of this chapter reasonable, additional polydisperse systems were used alongside the six samples discussed previously. The total solids concentrations ranged from 51 to 56 wt.%. The values presented in Fig. 4.7 were obtained using the same methodology described previously in this work. The correlation of the apparent viscosity at a shear rate of 400 s -1 , η 400 = τ 400 / . γ 400 , with particle size distribution is provided in Fig. 4.7. The d L /d s ratio is the particle size ratio (d 80 /d 20 ) and is expected to correlate with the apparent viscosity η. This seems to be borne out by the data: η decreased as the particle size expanded. This can be predicted using the theory of Farris when applied to suspensions of single-sized fractions of spherical particles [START_REF] Farris | Prediction of the Viscosity of Multimodal Suspensions from Unimodal Viscosity Data[END_REF]. For a given solids concentration, a wide particle size distribution will decrease the apparent viscosity of the mixture. The maximum packing fraction of a monomodal suspension is not influenced by the particle size. However, the maximum packing fraction increases when the particle size distribution (PSD) gets wider because the small particles can fill the void between the large particles. This relation between ϕ m with the PSD can be used to predict the Krieger-Dougherty equation, and hence, to predict a slurry's apparent viscosity. This justifies further research on the viscosity of multimodal suspensions to better understand the role of particle size distributions. From an industrial perspective, finding the best size proportions would be useful to optimize the transport process by reducing the resistance to the flow. In this context, air classification might be used before mixing the solids with water in the primary classification of mineral processing.

Laminar Pipe Flow

In the mineral-processing industry, it is essential to know the rheological properties of pumped slurries over extended distances in order to evaluate the key design parameters (flow rates, pressure gradients, etc.). In this subsection, mathematical models of the flow of time-independent and non-Newtonian fluids through a straight pipe in steady-state laminar flow regimes are described in order to highlight the influence of rheological parameters on flow behavior. In the following, it was assumed that the liquidsolid mixture was homogeneous and that there was no slippage between the suspension and the pipe wall.

When considering a unidirectional and axisymmetric flow of a fluid through a circular tube with radius R, the relation between the volumetric flow rate Q and the pressure gradient ∆p/L is obtained via integration of the Rabinowitsch-Mooney equation:

Q πR 3 = 1 τ 3 w τw 0 τ 2 . γ dτ, (15)
where Q is the total volumetric flow rate and τ w is the wall shear stress. At r = R, the relationship between the wall stress and the pressure gradient is as follows:

τ w = R 2 ∆p L . (16)
When the relationship between the shear stress and the shear rate raw data from the rheometer can be accurately modeled using explicit algebraic expressions, such as the rheological models that we proposed in Section 4.1.2, we can use them as a substitute for . γ in Equation ( 16), allowing the integral to be evaluated analytically. Similarly, cross-sectional velocity profiles for a given pressure drop can be found. The shear rate is expressed as . γ = du/dr. Assuming a no-slip condition u(R) = 0, this integrates to give the flow velocity distribution in the axial direction u(r ) with respect to the radial distance:

u (r) = R r . γ dr. (17)
In Fig. 4.8, the Herschel-Bulkley model was used to plot the flow curves of samples 4 and 6 in the laminar regime. As can be seen, S6 required a greater pressure gradient (or equivalently, τ w ) in order to reach the same mean flow velocity as S4 through a cross-section. For a mean velocity of 0.8 m/s, sample 6 theoretically needed nearly three times more pressure per unit length for a given transfer duty. The required pressure gradients per length for S4 and S6 were 336.84 Pa•s/m and 1107.26 Pa•s/m, respectively. Indeed, as seen in Table . 4.3, the curve fitting indicated that sample 6 had a higher value of dynamic yield stress. Moreover, Fig. 4.3 revealed a higher value of apparent viscosity over the range of shear rates where the Herschel-Bulkley model was applied. Thus, the design of slurry pipelines relies heavily on the rheological behavior of the transported suspensions, which needs to be optimally controlled. Fig. 4.9 shows the velocity distributions for different rheological models during the laminar flow of phosphate ore slurry through a 0.1 m pipe evaluated at a constant pressure gradient ∆p/L. The rheological parameters obtained via curve fitting using the raw data of sample 4 served as the basis to illustrate the effect of different rheological models on the velocity profile. In these plots, the value of the pressure gradient was adjusted to maintain a laminar flow regime [START_REF] Khataniar | Evaluation of Rheological Models and Application to Flow Regime Determination[END_REF][START_REF] Metzner | Flow of Non-Newtonian Fluids-Correlation of the Laminar, Transition, and Turbulent-Flow Regions[END_REF][START_REF] Bailey | A Generalized and Consistent Pressure Drop and Flow Regime Transition Model for Drilling Hydraulics[END_REF][START_REF] Peixinho | Laminar Transitional and Turbulent Flow of Yield Stress Fluid in a Pipe[END_REF].

As shown in Fig. 4.9, the velocity profile of a yield stress fluid presented a plug core flowing in the center of the pipe. In this zone, the applied stress was less than the yield stress. When the threshold was exceeded, the sheared fluid presented a curved velocity profile. For the Bingham and Casson models, the flat region occupied a major part of the pipe section. This could be explained by the fact that the curve fitting based on these two models resulted in high values of yield stress in comparison with the Robertson-Stiff and Herschel-Bulkley models, which enlarged the zone of the plug flow. Lower values of yield stress imply a reduction in the size of the plug core and a higher velocity magnitude [START_REF] Ohen | Golden Section Search Method for Determining Parameters in Robertson-Stiff Non-Newtonian Fluid Model[END_REF]. Fig. 4.9 clearly demonstrates that using different models to describe fluid flow may lead to completely dissimilar outputs. The accuracy of the mean velocity or pressure loss prediction is then closely related to the accuracy of the data and the closeness with which the slurry flow properties fit any chosen model. Performing a benchmark test with pipe flow experiments is crucial for validating the obtained results.

Conclusions

As a starting point, we proposed an experimental protocol for measuring the shear rheology of a concentrated phosphate slurry. The flow curves appeared to be non-Newtonian and exhibited shear-thinning behavior in all cases and could be described as yield pseudoplastics. Then, well-known rheological models were considered for the modeling. These equations were more suitable for describing the behavior of our mixture over the range of shear rates. Based on the results of this work, the following conclusions were reached. The Robertson-Stiff model was found to give the best description of the flow curves, yielding a value of R 2 higher than 0.97 in all cases. The Casson model, as a two-parameter model, offered a good alternative for this suspension. The apparent viscosity increased as the particle concentration increased. It was also observed that an expanded particle size distribution yielded a lower apparent viscosity for all cases. The re-sults show that apparent viscosity can be reduced at a constant loading by up to 20%. Later, it was found that the accuracy of the pressure loss prediction relied on the accuracy of the rheological model. Despite the fact that a significant amount of rheological data for the concentrated phosphate slurry were evaluated in this study, the data did not cover the complete range of the shear rates. Hence, we are aware that some conclusions are limited to the range of the data considered in this study. These findings can be exploited for better control of phosphate dispersion in biocomposites polymeric-based systems, which allow for better control of the filler dispersity and avoid harmful aggregation. Overall, shear rheology seems to be a suitable tool for controlling the phosphate filler properties and especially their polydispersities upstream of their processing. This part is the subject of an ongoing investigation that is to be published in our future works.

Shear Rheometry of Phosphate suspensions

Introduction

The objective of this study is to conduct a thorough analysis of the rheological characteristics of phosphate suspensions. The understanding of these properties is crucial for several industrial applications and the optimization of various processes. To achieve this goal, it is essential to establish an appropriate rheometry protocol that accurately reflects the real-world conditions of these suspensions. One important aspect of the rheometry protocol is to ensure that the suspensions are properly agitated either before or during the tests. This is necessary to avoid any potential settling or formation of clusters that could affect the test results and lead to inaccurate conclusions. Another crucial aspect is to perform the measurements over brief time intervals, as phosphate suspensions often exhibit rheological changes that need to be captured in the analysis. Finally, the geometry of the rheometer must be adapted to suit the specific requirements of phosphate suspensions.

The present section aims to provide a comprehensive understanding of the rheological properties of phosphate suspensions and to establish an effective rheometry protocol for their analysis.

Physical stability of phosphate suspensions

Materials and methods

Two types of materials were considered for testing:

• Suspensions collected directly during the industrial process.

• Extracted dry phosphate powders, allowing to prepare reconstituted suspensions.

These materials correspond to extractions in 3 different locations (F1 -F2 -F3)

The aim here is to describe the behavior of homogeneous water-solid mixtures. Given these conditions, to prepare the collected samples, we shake the container several times so that the supernatant water mixes well with the paste formed at the bottom due to sedimentation. To prepare the reconstituted suspensions, we mix a certain volume of water with a volume of powder in order to homogenize the suspension. Thus, the concentration of the prepared mixture is the ratio between the solid mass and the total mass of the material. It is to be noted that sedimentation is a major problem: reconstituted suspensions sediment quickly at rest. However, sampled suspensions sediment much less quickly than collected ones.

The shear rheology is measured using a couette geometry (coaxial cylinders) with an inner diameter of 25 mm and an outer diameter of 34 mm, i.e. a wide gap of 4.5 mm, which is about 15 times greater than the maximum grain size, which allows us to verify the assumption of the continuous medium. The surfaces of both inner and outer cylinders are roughened in order to avoid the wall-slip effects (roughness of about half a millimeter). A large gap (about 25 mm) is left between the bottom of the inner cylinder and the bottom of the outer cylinder to avoid the effects of sedimentation that may block the rotation of the cylinder.

In view of the problem of sedimentation, the rotational tests are carried out rapidly after pouring the homogenized suspensions into the geometry. To determine the behavior of the reconstituted suspensions, we imposed an increasing shear rate ramp followed by a decreasing shear rate ramp. For the collected suspensions, this protocol was refined and we finally settled on imposing a fairly fast stress ramp, typically 1 min 30s, so as to cover quickly enough a wide range of velocity gradients, before sedimentation or migration effects are too important.

It should be noted that it will be possible to study these effects in greater detail in the future, and/or to develop even more effective experimental procedures. For example, one can consider imposing creep at different stresses, using a sample renewed between each measurement, to be sure to have a homogeneous sample at each stress level, at least at the beginning of the test. We can also consider quantifying the migration by crossing the data obtained over time at different velocity levels. Finally, one can consider MRI tests to directly measure the local concentration and the local velocity gradient during the flow. In this context, it is important to characterize migration and sedimentation, as these are effects that are likely to occur significantly within the real flow. However, all these aspects are left aside in the present study which aimed at determining the main characteristics of the behavior of phosphate suspensions. The reconstituted suspensions sediment very quickly, which makes it very difficult to assess the behavior of such a suspension in a homogeneous state. Using the above measurement procedure, we obtain a first global overview of the behavior of this material as a function of its concentration, illustrated by Figure 4.10. We can observe that, at low solid concentration, the measurements are not conclusive at low speed as the applied constraints are weak. On the other hand, at higher speed, the material exhibits a Newtonian behavior, as expected with a diluted suspension of non-colloidal particles. As the concentration increases, a threshold appears, but this is most likely due to the effects of sedimentation. Indeed, in this case, the particles seem to disappear from the upper layers very quickly during the flow, and at the end of the test, we see that the concentration of the suspension is very high at the bottom of the Couette apparatus. As the inner cylinder is in contact with this compact area, it gives the impression that the suspension has a stress threshold. This effect is particularly noticeable when the concentration is increased, as the thickness of this concentrated paste increases. Then, at very high concentration (>65 wt.%), another effect appears: the suspension can no longer sediment as the particles are immediately in contact with each other. The behavior of this suspension is then that of a granular paste, which gives rise to shear localization. This induces a higher yield stress at the start of the test compared to the real flow threshold.

Assessment of rheometry tests

These measurements give us a glimpse of all the possible behaviors and experimental problems. Taking into account these challenges, it would be necessary to set up adapted procedures in order to determine the exact behavior of these suspensions in a homogeneous state. However, one may wonder what the interest of this would be with respect to the real industrial processes, since under real flow conditions these artefacts should be widely developed.

Collected samples

The collected suspensions also sediment significantly at rest, but much more slowly than the reconstituted samples. It is therefore possible to prepare a suspension that remains homogeneous for a period of time sufficient to accurately measure its rheological behavior. It is to be noted that, given the sedimentation effect, the obtained behavior does not describe all aspects of the real mixture behavior (sedimentation, migration, aggregation) and their impact on the flow characteristics in real conditions.

In Figure 4.11, the observed behavior is essentially that of a yield stress fluid. However, various problems (migration -sedimentation) tend to be developed during the tests, which affect more or less quickly the material, and thus its apparent behavior. These effects can be underlined during repeatability tests: by keeping the same sample, one can see that the apparent behavior of the material evolves in an unexpected way (see Figure 4.11). The main effect in this case is a reduction in concentration along the inner cylinder. Again, further studies would be needed to understand the origin of these effects and their impact on the behavior. In this study, we have limited ourselves to a first assessment of the behavior by applying a ramp of increasing stress on a freshly prepared (homogenized) material, considering that this material remains sufficiently homogeneous during a single measurement. As shown in Figure 4.12., the apparent behavior of the collected samples is very different from that of the reconstructed ones. Thus, regardless of the different artifacts, the yield stress of the collected mixtures is much higher than that of the reconstructed mixtures. This explains why sedimentation is much slower in the present case.

Questions may arise about to the origin of such difference in response to the applied stresses. Presumably, in collected sample, there is a network of interactions (in the form of "links") between the particles, resulting from the presence of the additive. This network is at the origin of the exhibited yield stress, and allows to reduce sedimentation. But it is also apparently very fragile: As soon as the bonds have been broken, the particles are simply dispersed in the liquid, which gives this simple viscous (Newtonian) behavior with a low viscosity (when compared to the high yield stress value). The addition of the additive is therefore particularly successful, as it slows down sedimentation, while leaving the material very liquid during the flow. However, it is still necessary to be able to impose sufficient pressure to exceed the threshold, and it is necessary to manage the migration or sed- imentation tendencies, which will be developed more and more as soon as the threshold is exceeded.

Improving the rheometry of phosphate suspensions

The above observations suggest that a more adequate rheometry protocol may be considered for the rheological characterization of phosphate suspensions, which would take into account sedimentation and migration effects.

For this purpose, we have tested the "vane" geometry, i.e. formed by a number of blades. In a first attempt to evaluate the rheological measurements obtained using this geometry, the tests were carried out with a vane consisting of 6 blades of 5 cm in height and 25 mm in diameter, and an outer rough cylinder of 37 mm diameter. This geometry has the advantage of having a high roughness, which prevent any wall slip effects along the moving geometry. Additionally, for a yield stress fluid, the generated flow is very close to that generated by the standard rough Couette geometry. This has been demonstrated by MRI (cf. Ovarlez et al, J. Rheology, 2011), and stems from the fact that the yield stress fluid is essentially "stuck" between the blades during the flow, whereas the material outside the blade envelope is sheared due to its friction with this trapped material. Experiments using the vane geometry with incressing and then decreasing stress ramps were performed, and the flow curves were compared to those obtained with the usual rough geometry. It can be seen from figure 4.13 that, contrary to all our previous tests, the downward and upward ramps overlap. This is a major result, which proves in principle that there is no significant evolution of the material during the flow. In addition, these results are repeatable with the same sample; and they are reproducible with different samples. This further confirms the absence of development of heterogeneities within the material. The plots in Figure 4.13 also show the apparent yield stress observed with the simple couette gemetry was indeed underestimating the actual threshold, which is, in turn, more accurately given by the stress plateau observed with the vane. However, it is also observed that, at a certain high shear rate value, the slope of the apparent flow curve is close to 2 in logarithmic scale. At a sufficiently high speed, the level of the suspension in the annular space between two coaxial cylinders increases abruptly. It is difficult to capture the flow characteristics in this case. However, we can think that this phenomenon is related to the Taylor-Couette instability, which leads to the formation of vortices. The Taylor-Couette instability criterion is expressed as:

ρ µ Ωr 0.5 1 (r 2 -r 1 ) 1.5 < 41 (4.1)
by introducing the apparent shear rate ( γ = Ωr 1 /(r 2 -r 1 )), we obtain:

ρ µ γr -0.5 1 (r 2 -r 1 ) 2.5 < 41 (4.2)
In the present case, we found that this instability occurs at an apparent shear rate of 160s -1 , and the apparent viscosity deduced from the flow curve in figure 4.14 is 0.135P a.s. Thus, the left term of this equation is 44, which is almost exactly the expected theoretical critical value. In practice, we wish to be able to control the flow by varying the velocity. This seems to be more suitable for systematic industrial measurements, as stress control requires to have an idea of the value of the yield stress, otherwise there is a risk of shearing the material too quickly, which can then be ejected from the geometry. However, when a decreasing shear rate ramp is applied, the obtained results are similar to the increasing stress ramp curves (see Figure 4.14), but the apparent flow curve no longer seems to have a threshold and is offset from the stress ramp curve. This phenomenon is attributed to a migration effect when very low velocities are imposed for too long.

In order to explain the difference between the obtained flow curves, it must be borne in mind that the behavior of our material is that of an unstable yield stress fluid : when the mixture is initially homogeneous at rest, we must apply a sufficiently high stress to make it flow, i.e. its yield stress value, but when it is slightly exceeded the material starts to flow and liquefies very quickly. On the other hand, if a very low shear rate is applied, which leads in theory to a shear rate value lower than the critical value reached at the yield stress, the fluid can not flow homogeneously. In this case, usually, the deformation is localized, so that the sheared zone is at the critical shear rate, while the rest of the material is not sheared. When the applied shear rate is lowered, only a part of the material within the couette gap is sheared, and this zone's thickness is reduced as the threshold is approached. At decreasing stress ramp, when the yield stress value is reached, the material stops flowing, which gives the observed stress plateau, similar to the plateau obtained throughout the increasing stress ramp. During a shear rate ramp test, rotation is sustained, even if the stress has reached the threshold, at apparent shear rates below the critical value. This implies that the deformation is localized. But we recall that as the stress was near the threshold shortly before, the sheared zone was already relatively thin. As a result, the localization becomes catastrophic, and the sheared layer of the material has a thickness of the order of the size of the particles. Finally, a slight migration effect leads to shearing just a very thin layer of liquid along the solid wall. This explains the drop in applied stress well below the threshold when the imposed shear rate is further lowered.

In view of these challenges, the measurement procedure may be improved by reducing the adverse effects. On one hand, in order for the Taylor-Couette instability to occur at the highest possible rotational speed, we can reduce the gap to 2.5 mm. Indeed, according to the above criterion, the Couette gap plays a critical role in this phenomenon. This value seems to be the smallest that can reasonably be used while preserving the validity of the continuity assumption. On the other hand, as the localization effect is mainly due to particles migration which is caused by the high velocities during the decreasing ramp (see Figure 4.14), it is preferable to impose an increasing shear rate ramp. By applying the above adjustments to the measurement protocol, we obtain a curve which is sufficiently close to the curve in increasing stress ramp, in which we have confidence as the disturbing effects have not had enough time to develop significantly (See Figure 4.15). In addition, The Taylor-Couette instability now seems to occur at much higher velocity gradients, above 1300 s -1 . We can therefore fit adequate models on the obtained curve.

Reproducibility of rheometry tests

Here, we assess the reproducibility of results for samples collected from different locations. The results are shown in Figure 4.16. We note a negligible variation of the results from one location to another. The apparent flow curves overlap reasonably well considering the usual uncertainty on rheometry measurements due to variation, from one test to another, of the boundary conditions, of the volume of material, etc. We can therefore extrapolate this inference to all the samples, and conclude that a good reproducibility of the results can be achieved using this protocol.

From Figure 4.16, we can take note of a certain "ageing" of the suspensions on a scale of several months: the flow curve observed, after about 6 months, has significantly higher stress values than the measurements on the new samples. Also, we observe that the lower the threshold, the earlier the Taylor-Couette instability occurs. This underlines a problem that we will have to take into account, as it limits the range of exploitable shear rates. For the present measurements, we use a 6 blades vane of 5 cm height im-mersed in an outer cylinder of 6 cm with different submerged vane heights. The stress is calculated based on the measured torque according to the standard formula resulting from the conservation of momentum:

Impact of the height of submerged vane

τ = C 2πhr 2 i (4.3)
It can be seen that the flow curve is slightly influenced by the chosen height, and no particular trend is observed with this parameter. Therefore, as long as the submerged height is largely superior to the Couette gap, there is no impact of this height. However, from a practical point of view, there is a significant uncertainty about the actual submerged height because it is difficult to inject the exact volume of slurry to reach that height, or to clean accurately afterwards to set the level at a particular height, or to measure that height.

Impact of the bottom gap

These measurements are performed by immersing the 5 cm long vane more or less high in the outer cylinder. We observe a slight tendency for the stress to increase when the sheared height is reduced, but no real divergence as long as we remain at a distance of more than 1 mm from the bottom, and as long as the standard protocol is used (fast increasing ramp). We conclude that the impact of the bottom gap is negligible. For these tests, the material must be prepared very carefully. First, the sample bottle is shaken strongly enough to resuspend the sedimented particles and obtain a homogeneous material. This constitutes a reference condition. Then the sample is quickly placed in the outer cylinder and the vane is quickly immersed as well in order to minimize the duration of this procedure (about 15 seconds). The end of the immersion corresponds to the "zero-time". From this point on, the usual shear rate ramp is imposed. For this particular experience, we wait for different durations "t" before starting the ramp. For some measurements, we observe a jamming of coarse particles (of more than 1 mm), which results in stress jumps (see Figure 4.19). However, the general trend is not significantly influenced. As shown in Figure 4.19, there was no observable effect of the resting time, no particular trend and notably no increase in threshold stress with "t" as what would be obtained with a thixotropic fluid. There is therefore no major impact of the rest time from the rheological point of view. But, it must be noted that a longer resting time leads to potentially significant sedimentation effects.

Impact of the imposed shear rate range

For this study, rheological tests are performed by imposing different ranges of shear rate over the same duration of 40s. In Figure 4.20, we observe ap- parent flow curves with similar trends, exhibiting a first phase of decreasing stress, then a plateau stress region and finally a significant rise at high shear rate values. It is also shown that the stress value determined at the plateau does not vary as the range of shear rates varies. However, the obtained flow curves indicate that starting the shear rate ramp at lower values leads to a growing difference between the initial yield stress and the plateau yield stress (up to 30%). The particle size distribution at the sheared zone is probably responsible for this phenomenon. This suggests that the yield stress determined at the plateau does not necessarily provide a perfectly representative value, as the migration during the increasing ramp may have disturbed the result. Thus, it seems preferable to start the test at a relatively high shear stress value, between 20 and 50 s -1 . This obviously has the disadvantage of reducing the observed range, but it reduces the gap between the first yield stress value and the minimum value reached afterwards and ensures that the fluid remains homogeneous. Note that this adverse effect disappears when starting directly at 50 s -1 (see Figure 4.20).

Impact of the gap width

Experiments following an identical procedure were performed with a vane using different gap widths. For all the measurements, the stress applied along the inner cylinder was plotted as a function of the apparent shear rate calculated by dividing the rotational velocity by a fixed fictitious gap width of 2.5 mm. It is worth saying that a gap of 1.5 mm seems too thin considering that there are some particles of the same size (1 mm) within the material. In these conditions, we can have at some moments jamming of these particles in the gap which yields important stress jumps. As shown in Figure 4.21, the obtained flow curves are almost perfectly overlapped. Obviously, under these conditions, if we had drawn the curve as a function of the apparent shear rate, assuming that the entire gap is sheared, we would have obtained curves that are significantly offset from each other. For example, the flow curve of a 1.5 mm gap would be offset horizontally by a factor of almost 4 compared to the flow curve of a 6.5 mm gap). This is an interesting result since it suggests that the gap width has no significant influence on the yield stress measured according to our approach, as the level of the threshold does not change when the gap width varies. On the other hand, the superposition of all the flow curves, even in the region where the stress increases significantly with the shear rate, is very surprising. In fact, one would expect a superposition as long as the material present within the gap is not completely sheared, as the stress increases. Indeed, in this case, the flow is identical for the different geometries, since the flowing suspension does not yet make contact with the outer wall, and the sheared thickness is the same for the different gap widths. Thus, the characteristics of the flowing zone of the material are the same for the different gaps as long as the sheared thickness is smaller than the smallest gap width tested. This phase eventually ends when the shear stress reaches the yield stress value at the outer wall, which corresponds to the case where the whole gap is sheared. It must be emphasised that, at higher shear rates, we can also notice a Taylor-Couette instability effect, which would as well prevent the curves from being offset even when the gap width is totally sheared.

Conclusions

In conclusion, the rheological characterization of phosphate suspensions presents several challenges due to sedimentation, migration effects, and the unstable nature of the yield stress fluid. To overcome these challenges, a vane geometry was employed, and various parameters were optimized to minimize adverse effects and ensure the accuracy of rheometry tests. The use of vane geometry with increasing stress ramps led to overlapping upward and downward ramps, demonstrating the absence of significant material evolution during flow. The apparent yield stress measured using this method was more accurate than with a simple Couette geometry. Adjustments to the measurement protocol, such as reducing the gap size and using an increasing shear rate ramp, further improved the results. Reproducibility of the rheometry tests was assessed, and the results showed negligible variation among samples collected from different locations. However, some aging of the suspensions was observed over time. Several other factors were investigated, including the submerged vane height, bottom gap, rest time prior to testing, imposed shear rate range, and gap width. The impact of these factors was found to be minimal, as long as certain conditions were met. Overall, this study demonstrates that the rheometry of phosphate suspensions can be improved by employing a vane geometry and optimizing the experimental protocol. This ensures the accuracy and reproducibility of the rheological measurements, which is essential for industrial applications.

Chapter 5

A New Rheological Model for Phosphate Slurry Flows

Introduction

Numerical modeling of multiphase flows along with complex rheology has particularly drew considerable attention during the last two decades and the study of their behavior has led to a large literature. Such complex phenomena are encountered in many industrial and engineering applications, including evaporation and condensation [START_REF] Bureš | Direct numerical simulation of evaporation and condensation with the geometric VOF method and a sharp-interface phase-change model[END_REF], flooding events [START_REF] Zhang | Machine learning strategy for viscous calibration of fully-nonlinear liquid sloshing simulation in FLNG tanks[END_REF][START_REF] Saghi | The influence of flexible fluid structure interactions on sway induced tank sloshing dynamics[END_REF], chemical and nuclear reactors [START_REF] Silvi | Understanding dry-out mechanism in rod bundles of boiling water reactor[END_REF], fluidized bed [START_REF] Vångö | Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds[END_REF], combustion and fuel atomization [START_REF] Giussani | A three-phase VOF solver for the simulation of in-nozzle cavitation effects on liquid atomization[END_REF][START_REF] Di Iorio | Investigation on sub-23 nm particles and their volatile organic fraction (VOF) in PFI/DI spark ignition engine fueled with gasoline, ethanol and a 30% v/v ethanol blend[END_REF]. Recently, researchers have devoted a considerable amount of research efforts into developing and implementing sophisticated mathematical techniques, relevant in the whole range of multifluid and multiphase flow problems. Being the two main families of solution methods, Level-Set (LS) and Volume-Of-Fluid (VOF) proposed respectively by Sussman et al. [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] and Hirt et al. [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] are the two most commonly employed approaches which have assisted in better understanding of the underlying physics governing the multiphase flows in different fields. Given the particularities of each model, VOF and LS have been widely used to investigate a broad range of engineering problems such as Rayleigh-Taylor Instability [START_REF] Kraposhin | A hybrid pressure-based solver for nonideal single-phase fluid flows at all speeds[END_REF], dam failure phenomena [START_REF] Yu | Numerical simulation of dam-break flow impacting a stationary obstacle by a CLSVOF/IB method[END_REF][START_REF] Jafari | Numerical simulation of wave interaction with porous structures[END_REF][START_REF] Garoosi | Experimental and numerical investigations of dam break flow over dry and wet beds[END_REF], bubble rising/nucleation [START_REF] Zhang | Dynamics of bubble formation on submerged capillaries in a non-uniform direct current electric field[END_REF][START_REF] Li | Numerical investigation of an underwater explosion bubble based on FVM and VOF[END_REF][START_REF] Das | Influence of surface tension on bubble nucleation, formation and onset of sliding[END_REF], droplet impact on both dry and wet surfaces [START_REF] Yang | Study on dynamic behavior of water droplet impacting on super-hydrophobic surface with micro-pillar structures by VOF method[END_REF][START_REF] Hanene | A numerical model of an immiscible surfactant drop spreading over thin liquid layers using CFD/VOF approach[END_REF][START_REF] De Lima | Numerical simulation of a water droplet splash: Comparison between PLIC and HRIC schemes for the VoF transport equation[END_REF][START_REF] Yun | Ellipsoidal drop impact on a single-ridge superhydrophobic surface[END_REF], metallurgical engineering [START_REF] Kalifa | Fluid flow phenomena in metals processing operations: Numerical description of the fluid flow field by an impinging gas jet on a liquid surface[END_REF], to mention only a few. In this study, we use the VOF method to investigate water-slurry modeling in a horizontal pipeline.

Viscoplastic or yield-stress fluids are materials which behave as solids as the shear rate tends to zero, and as liquids beyond a certain critical shear stress level (τ 0 ). In particular, simple yield-stress fluids are materials that are both non-thixotropic and inelastic, and characterized by a shear stress that depends only on the applied shear rate, and materials [START_REF] Frigaard | Simple yield stress fluids[END_REF]. The most popular simple viscoplastic model was first proposed by Bingham [START_REF] Bingham | Fluidity and Plasticity[END_REF] and is defined by a yield stress value and a constant plastic viscosity representing the slope of the shear stress versus shear rate curve. The Bingham constitutive equation can be expressed as follows:

   γ = 0 if τ ≤ τ y τ = µ + τ y γ γ if τ > τ y (5.1)
where τ and γ are the stress tensor and the shear rate tensor respectively. While the terms τ 0 , and µ indicate the yield stress and the plastic viscosity scalar values. τ = 1 2 τ : τ and γ = 1 2 γ : γ denote the second invariants of τ and γ, respectively. The strain rate tensor γ is defined as γ = ∇u + (∇u) T , u being the velocity vector. The flow domain of a yield stress fluid is divided into a yielded region (τ > τ 0 ) and an unyielded (τ ≤ τ 0 ) region, separated by the yield surfaces where τ = τ 0 . Two well known generalizations of the Bingham fluid model have been proposed by Herschel-Bulkley [START_REF] Herschel | Konsistenzmessungen von gummibenzollösungen[END_REF], and Casson [START_REF] Casson | A flow equation for pigment-oil suspensions of the printing ink type[END_REF], formulated respectively as:

Herschel-Bulkley: τ = K γ(n-1) + τ 0 γ γ if τ > τ 0 (5.2)
Casson:

τ = µ γ + √ τ 0 2 γ γ if τ > τ 0 (5.3)
where K is the consistency and n is the power-law index. The flexibility of the Herschel-Bulkley model to fitting various experimental data makes it very popular among other yield stress models. Moreover, the Herschel-Bulkley model can be reduced to power-law and Casson models when setting the yield stress to zero τ 0 = and the power-law index to n = 1, respectively. According to [START_REF] Tanner | Rheology of noncolloidal suspensions with non-Newtonian matrices[END_REF], the power law, the Bingham plastic, the Casson and the Herschel-Bulkley models remain the four widely used models for describing the viscous properties of suspensions.

Various research works have studied slurry flows as a continuum non-Newtonian fluid based on the equivalence assumption and have made use of either Bingham model or Herschel-[-35]Bulkley (H-B) model to represent the rheological properties [START_REF] Zhang | A semi-implicit discrepancy model of Reynolds stress in a higherorder tensor basis framework for Reynolds-averaged Navier-Stokes simulations[END_REF][START_REF] Dai | Viscometric functions for noncolloidal sphere suspensions with Newtonian matrices[END_REF]. Hamza et al. (2018) [START_REF] Belbsir | Study of the rheological behavior of the phosphate-water slurry and search for a suitable model to describe its rheological behavior[END_REF] have investigated in their work the rheological behavior of the phosphate-water slurry in an attempt to determine a model capable of describing its flow behavior. They came to the conclusion that the Herschel-Bulkley model is a suitable model for representing the rheological behavior of the phosphate slurry for low concentrations (less than 38.45 wt%). For higher concentrations in the range of 34.24-46.03 wt%, the Bingham model was more adequat. Finally, the Casson model was the best fitting model for calculating viscosity and yield stress for the range of concentrations (46.03-57.27 wt%). In addition, based on the experimental results, they conclude that the phosphate slurry rheological behavior tends to be a dilatant behavior beyond a concentration of 38.45 wt%.

The aim of this work is to develop a model that could describe the rheological behavior of the phosphate slurry over a wide range of concentrations. We will show that the model is capable of reproducing the rheological behavior that phosphate suspensions may exhibit in both low and high concentrations. This is achieved by comparing the output of the model with those of the models suitable for modeling the rheological behavior of phosphate slurry. A methodology for obtaining the model's parameters is presented in details. Whilst the term "multi-phase flows" covers the whole spectrum of gas-liquid, liquid-liquid, gas-solid, liquid-solid, gas-liquid-solid and gas-liquid-liquid systems, our concern in this work is to investigate the simultaneous co-current flow of a simple liquid and a non-Newtonian suspension of solid particles. The new model presented in this chapter should be of interest to broad and diverse areas of application. Furthermore, the subject may be relevant both to theoretical mathematicians and practising engineers with a wide range of backgrounds.

Materials and Methods

The flow behavior of the materials encountered in many chemical and process engineering applications is highly influenced by the characteristics of their components, including their non-Newtonian properties, shape, size and concentration of the suspended particulates, the applied shear rate and the geometry of the system. In general, the flow behavior of such systems is so complex that theoretical treatments, which tend to apply to highly idealized problems, have proved to be of little practical utility. Therefore, flow investigations rely heavily on analyses of the behavior of such systems in practice based on experimental work and well established assumptions.

The New Rheological Model

Over the years, to model the stress-deformation behavior, several empirical expressions have been proposed as a result of straightforward curve-fitting exercises and different yield criteria have been used. The majority of the available literature is empirical and still need to be developed, partially due to the difficulty of obtaining accurate, reliable data on yield stress materials over a wide range of shear rates.

In this study, we seek to investigate the rheological behavior of the Phosphate slurry, which is a mixture of insoluble particles and a continuous phase of water. The obtained rheograms of this material show a non-Newtonian flow behavior exhibiting a yield stress as the shear rate tends to zero. A two-branch model equation is presented to embrace the non-linear flow of phosphate slurry over a very large range of shear rates. In one dimensional steady shearing motion, it is written as:

τ = a( γ + γ0 ) b if γ < γc τ c + γη ∞ if γc ≤ γ (5.4)
where the parameters a and b can be considered similar to the the consistency and the power-law index of the Herschel-Bulkley model, respectively. On the other hand, the parameter γ0 has a different interpretation than the model parameters of both Herschel-Bulkley or Bingham model. In this model, γ0 can be regarded as a correction to the shear rate rather than the shear stress and the yield stress is defined as τ y = a γb 0 . The linear part of the curve is described using the parameters η ∞ and τ c , which represent respectively the consistency index [Pa•s], the yield stress [Pa]. The transition from one branch to the other depends on the critical value of shear rate which is experimentally defined as γc = 400 s -1 in the particular case of phosphate slurries considered in this work.

To determine the rheological coefficients of each model equation, we perform a curve fitting based on the least square method. Curved relationships between variables are not as straightforward to fit and interpret as linear relationships. Given M data pairs ( γm , τ m ) where m ∈ {1, ..., M }, the parameters of Equation (5.4) need to be determined. We consider the general case where the model can be formulated as τ = f (θ, γ) where θ is a vector of p parameters. In what follows, the least squares approach is applied to determine the parameters which to minimise the following expression:

min θ M m=1 [τ m -f (θ, γm )] 2 (5.5)
where the pairs ( γm , τ m ) are observed. Given the observed data pairs {( γ1 , τ 1 ), ..., ( γM , τ M )}, we may define the error associated to the second branch of Equation (5.4), τ = τ c + γη ∞ , by:

E(η ∞ , τ c ) = M m=1 (τ m -(τ c + γm η ∞ )) 2 (5.6)
Our goal is to define the values of η ∞ and τ c that minimize the error function. Following the least-squares method, we should find the values of (η ∞ , τ c ) such that:

∂E ∂η ∞ = 0, ∂E ∂τ c = 0 (5.7)
It is to be noted that, in this case, we do not have to worry about boundary points: as |η ∞ | and |τ c | become large, the fit will clearly deteriorate. Thus we do not need to check on the boundary. Differentiating E(η ∞ , τ c ) yields:

∂E ∂η ∞ = M m=1 2(τ m -(τ c + γm η ∞ ))(-γm ) (5.8) ∂E ∂τ c = M m=1 2(τ m -(τ c + γm η ∞ )) (5.9)
Setting ∂E/∂η ∞ = ∂E/∂τ c = 0 and dividing by 2 yields:

M m=1 (τ m -(τ c + γm η ∞ )) γm = 0 (5.10) M m=1 (τ m -(τ c + γm η ∞ )) = 0 (5.11)
We may rewrite these equations as:

M m=1 γ2 m η ∞ + M m=1 γm τ c = M m=1 τ m γm (5.12) M m=1 γm η ∞ + M m=1 1 τ c = M m=1 τ m (5.13)
We have obtained that the values of η ∞ and τ c which minimize the error (defined in Equation (5.6)) satisfy the following matrix equation:

M m=1 γ2 m M m=1 γm M m=1 γm M m=1 1 η ∞ τ c = M m=1 τ m γm M m=1 τ m (5.14)
We will show the matrix is invertible, which implies we find that,

η ∞ τ c = M m=1
detA = M M m=1 γ2 m -(M γm ) 2 = M 2 1 M M m=1 γ2 m -γm 2 (5.18)
where the last equality is derived from simple algebra. Hence, while γm are not all equal, detA will be nonzero and A will be invertible. Therefore, we note that as long as γm are not all equal, the best-fit values of η ∞ and τ c are achieved by solving a system of linear equations in the Equation (5.15). As for the first branch of Equation (5.4), it requires an initial estimate for one of the three parameters. We can use the following procedure for the estimation of these parameters. γ0 can be estimated first using the following equation:

γ0 = γmin γmax -γ * 2 2 γ * -γmin -γmax (5.19)
where γ * is the shear rate corresponding to the geometric mean of the maximum and minimum shear stresses:

τ * = (τ min τ max ) 1/2 (5.20)
Then, taking the logarithm of the first branch of Equation (5.4):

log(τ ) = log(a) + b log( γ + γ0 ) (5.21)
and using the same methodology described above, one can obtains the two remaining parameters by change of variables then calculation.

The proposed model, which is based on the formulation described earlier, has been implemented using the OpenFoam numerical code. This powerful open-source framework offers a wide range of options for modeling the behavior of fluids, including a library of viscosity models such as Bingham, Herschel-Bulkley, and Casson. These models are defined in terms of strain rate and can be easily customized by the user to suit the specific requirements of their application. Additionally, the transportProperties dictionary allows for further flexibility in defining the rheological properties of the fluids being modeled. The implementation of the proposed model using the OpenFoam framework thus enables a more accurate and comprehensive simulation of multiphase flows and complex rheology in various industrial and engineering applications.

Numerical Modelling

The new model is applied to multiphase flow to account for the true behavior of phosphate slurry and its deformations that occur during this engineering instance. The simplest way to solve a two-fluid flow is to track the interface between the two fluid phases with a simple α indicator, often referred to in the literature as a concentration function. This method is called VOF (Volume-Of-Fluid) method, other approaches include Level-Set, Front Tracking, etc. twoLiquidMixingFoam [START_REF] Zhang | Large eddy simulations of 45 inclined dense jets[END_REF][START_REF] Krpan | Simulation of turbulent wake at mixing of two confined horizontal flows[END_REF] is a well-established solver for the simulation of flows in which two incompressible fluids are present. A separate surface interface can be defined based on the VOF method in OpenFOAM, which is a free and open-source parallel processing software supported by a large user community [START_REF] Jasak | OpenFOAM: A C++ library for complex physics simulations[END_REF][START_REF] Zhang | A semi-implicit discrepancy model of Reynolds stress in a higherorder tensor basis framework for Reynolds-averaged Navier-Stokes simulations[END_REF]. The different phases are represented in the domain by their phase fractions using the VOF [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF][START_REF] Palmore | A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows[END_REF]. This method allows to reduce considerably the computational costs.

An algebraic VOF method is used in the twoLiquidMixingFoam solver, which is a modified version of the VOF method of [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] by adding a diffusion term. The phosphate slurry consists of solid particles and water. Hereafter, we considered the homogenenous regime of the slurry with constant density and viscosity across the pipe. Therefore the two-phase system is the slurry phase and the water batch.

Governing Equations

The sum of the volume fractions of all phases in a cell is equal to unity and is given by the Equation (5.22) below

α s + α w = 1 (5.22)
where α s and α w denote the volume fraction of the slurry and water phases respectively. Both phases are assumed to be incompressible. The continuity equation is given as follows

∇ • ⃗ v = 0 (5.23)
where ⃗ v represent the velocity. The momentum equation solved for both phases during the simulation is given by,

ρ ∂⃗ v ∂t + ρ⃗ v(∇ • ⃗ v) = -∇P + ∇ • τ + ρ⃗ g + ⃗ F cs (5.24)
where P and g are the pressure and acceleration due to gravity, respectively. τ represents the shear stress tensor; for Newtonian fluids, it is a linear function of the shear rate given by where α s is the cell fraction of the slurry phase, Γ t is the effective diffusivity, D AB represents the coefficient of molecular diffusion, ν t is the turbulent viscosity, and Sc t denotes the turbulent Schmidt number. The coefficient of diffusion D AB and the turbulent Schmidt number Sc t were set to 1.0 × 10 -6 m 2 /s and 1.0, respectively. The phase fraction in a cell is related to the RANS model by density as follows ρ = α s ρ s + α w ρ w = α s ρ s + (1 -α s )ρ w (5.32) where ρ s represents the density of the slurry that corresponds to the α s phase fraction, while ρ w is the density that corresponds to the water phase with the α w phase fraction. The calculated density ρ in the Equation (5.32) should not be considered as a real density of the fluid but as a density of the phase mixture that occupies a cell in the finite volume framework.

τ = µ γ = µ ∇⃗ v + (⃗ v) T - 2
The RANS model uses the density ρ in the conservation of momentum equation. When α s is equal to 1, the represented phase is the slurry while the water phase α w is defined as 0, and vice versa. Any value between 1 and 0 describes a mixture and can also be interpreted as a percentage of the amount of this or that fluid present in a cell. The Equation (5.30) can be interpreted as an advection-diffusion equation for the α s phase.

Turbulence Equations

However, we should mention that all the constants and wall function of the k -ωSST model have been derived experimentally by measurements on Newtonian fluids [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. Moreover, they are not adapted to non-Newtonian fluids. In the present work, we have used k -ωSST model to simulate the flow of water pushing the slurry batch, thus the turbulence is applied to the water batch. where Gk is the generation of turbulence kinetic energy due to the mean velocity gradients. Gω denotes the generation of the specific dissipation rate. Γ k and Γ ω are the effective diffusivity of k and ω respectively. Y k and Y ω represent the dissipation of k and ω due to turbulence. D ω is the crossdiffusion term. S k and S ω represent source terms set to be zero in this study. All the above terms are calculated and specified in [START_REF] Launder | The numerical computation of turbulent flows[END_REF][START_REF] Henkes | Naturalconvection flow in a square cavity calculated with low-Reynoldsnumber turbulence models[END_REF][START_REF] Yakhot | The renormalization group, the ϵ-expansion and derivation of turbulence models[END_REF][START_REF] Menter | Turbulence Modeling for Engineering Flows[END_REF].

Results and Discussions

Rheological Evaluation of the New Model

In the present study, we used phosphate ore samples of different grades and concentrations of 51, 54, and 56 wt% in water. These samples were labeled S1 to S6 and different primarily in their solids concentrations. The details on the rheological data and samples characteristics can be found in [START_REF] Maazioui | Rheological Characterization of a Concentrated Phosphate Slurry[END_REF].

The rheograms of phosphate slurries were obtained, at room temperature, by applying a ramp of predefined shear rate decreasing from 1000 s -1 to 112 s -1 , and measuring the corresponding shear stress. The diameters of the rotating bob and the cylindrical cup are respectively 38.713 mm and 44 mm. Prior to each measurement, the suspensions were stirred carefully in a vessel to wipe out material memory and obtain the same initial conditions for both samples and then rapidly filled into the external cylinder. The effects of wall slip were not completely prevented using this geometry, so any reproduction of the resulting rheograms should be carried out with caution. The experimental data are shown in Figure 5.1. The accuracy of the New model is checked by fitting the phosphate slurry data. Figure 5.1 shows a good agreement with the experimental data for six samples of the phosphate slurry.

Numerical Results

OpenFoam 9.0, the open source CFD framework was used to solve the flow equations, turbulence models and transport equations. It is a very flexible CFD code where each component can be personalized to meet the user's needs. In single-phase flow calculations, the SIMPLE algorithm is used to couple pressure with velocity, whereas the PIMPLE algorithm is used in bi-phase simulations. In the single-phase simulations, the relaxation factors used are 0.5 for U , k and ϵ and 0.3 for the pressure. While, in the two-phase flow case, the relaxation factor for all variables is 1. GAMG is the linear solver used for pressure, and the smooth solver was adopted for U , k, ϵ and α.

Model Implementation

The new model is implemented in the OpenFoam 9.0 source code. The test case "lock-Exchange", implemented in OpenFoam library, is used to test the new model. Since performing numerical simulations of a yield stress material flow is not a straightforward task, the regularization method was widely used in the last decades, see [START_REF] Maazioui | Numerical Study of Viscoplastic Flows Using a Multigrid Initialization Algorithm[END_REF][START_REF] Guadagni | On a Casson fluid motion: Nonuniform width symmetric channel and peristaltic flows[END_REF] and references therein. The regularization approach includes methods which approximate viscosity value by one regularized and smooth constitutive equation, which is well determined regardless of the shear rate magnitude. The regularized equation treats the whole material domain as a fluid of variable viscosity and locally assigns a large but finite value of viscosity to the unyielded regions. Similary, In OpenFoam, the material is modelled for low strain rates as a very viscous fluid with a high viscosity µ 0 defined by the user. Beyond a threshold in strain-rate corresponding to threshold stress, the viscosity is described by a the constitutive equations. Thus the implementation is formulated as: µ = min(µ 0 , µ( γ)). The physical characteristics of the two phases used water and sludge are represented in Table 5.1. The multiphase solver used twoLiquidMixingFoam is described in Table 5.2. For the purpose of comparison, three rheological models already present in the OpenFOAM library were employed-the Newtonian, Casson, and Herschel-Bulkley models. The S1 pulp sample was used in this comparison, with the values defined in Table 5.3 being assigned to the variables 'a' and 'b' respectively. An examination of Figures 5.2-5.5 reveals significant differences in the distribution of sludge concentration within the column when using the proposed new model in comparison to other models such as the Newtonian, Casson, and Herschel-Bulkley models. While the latter models tend to result in the mixing and formation of eddies between the water and sludge, the new model maintains a clear physical interface between the two fluids. This prevents mixing over time and instead results in a clear layering of the denser sludge at the bottom of the pipe and the less dense water at the top. This behavior can be clearly observed in the figures and highlights the potential advantages and unique properties of the new model in the simulation of multiphase flows and complex rheology in various industrial and engineering applications. In order to validate our new model for phosphate slurry flow visualization, we compared it to the well-established Newtonian model using high mesh resolution (2 million cells), due to the limited experimental data available on this subject. The simulations were performed in a pipe with a diameter of D = 0.9 m and a length of L = 50 m, with a mean flow velocity of U = 5 m/s to achieve high turbulent flow. This was carried out in order to replicate the conditions often found in industrial settings where phosphate slurry is transported. Using the least squares method, we obtained consistency index of 1.56 and flow index of 0.01, which correspond to a concentration of 56 wt%. These results demonstrate that our new model is a reliable and accurate representation of the rheology of phosphate slurry flow and can be used in industrial settings to better predict and control the flow of this fluid. Upon examination of the figure, it is clear that the new model provides significantly improved results in comparison to the Newtonian model. One of the most notable differences is the absence of mixing at the interface between water and slurry in the case of the new model. This indicates that the new model is able to capture the behavior of the slurry more accurately and can be used to make more accurate predictions about the behavior of the slurry in industrial settings.

Additionally, the new model has the ability to predict the slurry concentration distribution at different sections of the pipe, which can be used to optimize the design and operation of industrial systems that transport phosphate slurry. The new model can also be used to improve the control of the flow of the slurry in real-world applications.

In conclusion, the new model for phosphate slurry flow visualization has been validated and compared to the Newtonian model with high accuracy. The results of the simulation demonstrate that this new model can be used to better predict and control the flow of phosphate slurry in industrial settings.

Conclusions

A new rheological constitutive relation for modeling the phosphate slurry flows was presented in this study. This model is based on five constant pa-rameters, which were determined through rheological measurements. The validity of the model was assessed by comparing its outputs with data derived from stress test measurements. It was found that there was a good fit between the model and the measurements in the stress versus strain rate curves. Furthermore, it was shown that the new model gives good numerical results when implemented in the OpenFoam software, version 9.0, and tested with an existing tutorial case model. The model was also validated numerically by comparing its concentration distribution results with those of the Newtonian model for high mesh resolution. It was found that the new model provides improved results compared to the Newtonian model. Additionally, the new model was compared to the Herschel-Bulkley model for a low concentration of C = 38.45 wt% and to the Casson model for C = 56 wt%. It was found that the new model was able to simulate the rheological behavior of phosphate slurry flows accurately for both low and high concentrations.

The main advantages of this new model are that it can produce improved fits of the rheological data compared to other existing models, and that it provides a description of the rheological behavior of the phosphate suspension for a wide range of concentrations. This makes it a valuable tool for understanding and predicting the behavior of phosphate slurry flows in industrial settings and for optimizing the design and operation of systems that transport phosphate slurry.

Chapter 6

Concluding remarks and perspectives

This doctoral thesis, through an in-depth exploration of the rheology of phosphate slurry and the development of innovative numerical solutions, has contributed significantly to our understanding of non-Newtonian viscoplastic flows in the context of industrial applications.

The incorporation of the multigrid initialization algorithm into numerical solvers has enhanced the convergence rate, efficiency, and overall computational performance. This contribution can help overcoming the limitations of single-grid methods and providing a firm foundation for future studies exploring other non-Newtonian problems and complex geometries. The proposed experimental protocol for measuring the shear rheology of concentrated phosphate slurries is an equally substantial contribution, yielding important insights into the behavior of these mixtures and pointing toward more accurate, efficient solutions for their management.

The use of a vane geometry and subsequent optimization of rheometry tests further refined our understanding of phosphate suspensions and shed light on achieving more precise measurements of rheological characteristics. The robustness of our methodologies is clearly demonstrated in our rheological model, which not only aligns well with experimental findings and outperforms existing models in terms of predictive accuracy, but also successfully reproduces the flow dynamics of batches in a pipeline system.

Undeniably, the methodologies proposed in this thesis for addressing the intricacies of viscoplastic flows and rheology of phosphate ore slurry have yielded substantial and noteworthy outcomes. Nonetheless, the pursuit of knowledge is relentless, and so the research continues! As we look to the future, we envisage several potential extensions of the work presented here.

Future work may focus on extending the range of shear rates examined, with a particular focus on extremely low and high shear rates. Investiga-tions may also be expanded to other suspensions to verify further the advantages of the proposed rheological model. Additionally, with the increasing drive towards sustainable industrial processes, an interesting perspective could be the exploration of phosphate slurry interactions with environmentally friendly, biocomposite polymeric-based systems. This could pave the way for new strategies to improve filler dispersity and avoid harmful aggregation, advancing the utility and sustainability of phosphate ore processing and application.

Finally, this research encourages the continual adaptation and optimization of the rheological model proposed here to keep pace with evolving industrial needs and advancements in technology. Further computational enhancements and innovative algorithmic solutions may also be sought, improving the efficiency and the accuracy in the prediction and management of phosphate slurry flows.

In closing, this research has laid the groundwork for an innovative and constructive path in the study of phosphate slurry flows and viscoplastic behavior. The tools and methods we've developed hold promising potential to guide and shape ongoing and future research in this field.

2. 3

 3 Illustration of the motion and straining of a fluid element . 2.9 Effect of particle size distribution on the maximum packing density for spheres. Polydisperse suspensions have a higher maximum packing density ϕ m ax compared to monodisperse suspensions . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Papanastasiou regularization for the Bingham model: stress magnitude (left) and apparent viscosity (right) as functions of the magnitude of the shear rate tensor γ. For increasing m, we recover a closer approximation to the actual Bingham model. Note that logarithmic axes are employed in order to highlight the behavior in the low shear rate limit. . . . . . 3.2 Two-dimensional orthogonal grid. N , S , E and W correspond to neighbor cells of cell P; n , s , e and w denote cell P faces; ∆x and ∆y are cell P dimensions in the x and y spatial coordinates; δy n , δy s , δx e and δx w correspond to cell-center to cell-center distances from cell P to neighbor cells. . . . . 3.3 Illustration of different multigrid cycle strategies: (a) Vcycle, (b) W-cycle. . . . . . . . . . . . . . . . . . . . . . . 3.4 Geometry of a cavity with moving lid . . . . . . . . . . . . 3.5 Geometry of a two-dimensional pipe of circular cross section. Dimensions are height H and length L. . . . . . . . . . . . 3.6 Sections of velocity along the vertical (a) and horizontal (b) mid-planes, Re = 10 3 . . . . . . . . . . . . . . . . . . . . . 3.7 Principal vortex position for various Bn numbers at Re = 10 3 . Reference results retrieved from Vola et al. [73] and Prashant et al. [72]. . . . . . . . . . . . . . . . . . . . . . . 3.8 Streamline contours in Bingham flow for Bn = 2 . . . . . . 3.9 Streamline contours in Bingham flow for Bn = 10 . . . . . 3.10 Streamline contours in Bingham flow for Bn = 100 . . . . ix LIST OF FIGURES x 3.11 Numerical results: the developing velocity profile of a Bingham fluid entering the pipe, Re = 10 2 and Bn = 10. . . . . 3.12 Velocity profiles at x = 0.5L calculated with m = 100 (solid lines), m = 200 (dashed lines), m = 400 (dotted lines), and m = 800 ( red dotted lines). For Re = 10 2 , (a) Bn = 1, and (b) Bn = 10. . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Apparent viscosity profiles along the vertical section x = 0.5L calculated with m = 100 (solid lines) and m = 800 ( red dotted lines), Re = 10 2 and Bn = 10. . . . . . . . . . . 3.14 The L ∞ norm of u and v residuals as a function of the number of SIMPLE iterations on the fine grid, for Bn = 0, 1, and 10 (Re = 10 2 ). Results are shown for simulations with single grid (SG) in red lines and multigrid (MG) initialization in blue lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Cumulative particle size distribution of the concentrated phosphate slurry (CPS) samples in terms of the percentage of cumulative volume versus particle size. S1-S6 signify the sample labels. . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Searle principle: the instrument's motor rotates the measuring bob inside an immovable cup filled with slurry (α cyl.cone = 120 • ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Representative rheograms showing (a) the shear stress as a function of the shear rate of CPS samples at various solids concentrations and a fixed temperature of 25 • C and (b) the apparent viscosity as a function of the shear rate under the same conditions. . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Comparison of rheological model fits for (a) S1 and (b) S2. 4.5 Comparison of rheological model fits for (a) S3 and (b) S4. 4.6 Comparison of rheological model fits for (a) S5 and (b) S6.4.7 Apparent viscosity at a constant shear rate of 400 s -1 as a function of the D-value ratio of the large (d L ) and small (d S ) particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Pipe pressure loss analysis of the Herschel-Bulkley model for S4 and S6. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Axial laminar velocity profiles that were calculated based on different rheological models. . . . . . . . . . . . . . . . . . 4.10 Flow curves obtained during ramps of increasing then decreasing shear rate, with suspensions reconstituted at different solid concentrations. . . . . . . . . . . . . . . . . . . . 4.11 Repeatability tests (up-down ramps) using collected samples (F2). Successive tests on the same sample (left in place) show artifacts evolving with time (migration -sedimentation). . 4.12 Ramp of increasing then decreasing shear stress using collected samples and reconstituted samples at the same solid concentration (55 wt.%). . . . . . . . . . . . . . . . . . . . 4.13 Different tests performed on the two collected phosphate suspensions with a rough surface Couette geometry (described above) and a vane geometry (25 -37 mm). In all cases, increasing and then decreasing stress ramps are applied. The dotted curve has a slope of 2. . . . . . . . . . . . . . . . . 4.14 Ramp of increasing then decreasing stress, and ramp of decreasing shear rate. . . . . . . . . . . . . . . . . . . . . . . 4.15 Comparaison between ramps of increasing stress and shear rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Apparent flow curves for 4 different samples obtained by imposing a ramp of increasing shear rate, with a Vane geometry of 6 blades, 5 cm height immersed in a cylinder of 11 cm height and a Couette gap of 2.5 mm. . . . . . . . . . . 4.17 Apparent flow curves for different submerged vane heights.4.18 Apparent flow curves obtained by imposing a ramp of increasing shear rates, with a Vane geometry 6 blades -5 cm high immersed in a outer cylinder located at different heights from the bottom. . . . . . . . . . . . . . . . . . . . . . . . 4.19 Apparent flow curves (vane 5 cm -5 mm from the bottom -2.5 mm couette gap) obtained for different rest times prior to the tests. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Apparent flow curves (vane 5 cm -5 mm from the bottom -2.5 mm couette gap) obtained for different ranges shear rate imposed during a increasing 40 s ramp. . . . . . . . . . . . 4.21 Apparent flow curves (vane 5 cm -5 mm from the bottom) obtained for different gap widths. . . . . . . . . . . . . . .

1

 1

.

Figure 2 . 1 :

 21 Figure 2.1: The surface forces acting on a differential surface element expressed in terms of the stress tensor

Figure 2 . 2 :

 22 Figure 2.2: Stress components on three faces of fluid element

  2.3. After a time interval ∆t, the position of these vertices is marked as A', B', C' and D'. Fig. 2.3 illustrates the changes in the positions of the vertices resulted from of (a) motion without any deformation and (b) deformation of the fluid element.

Figure 2 . 3 :

 23 Figure 2.3: Illustration of the motion and straining of a fluid element

Figure 2 . 4 :

 24 Figure 2.4: Illustration of the deformation of the fluid element

Figure 2 . 5 :

 25 Figure 2.5: Illustration of simple shear flow: The bottom plate is fixed, while the upper plate is moving at a steady speed of V

  18) and has dimensions [L] 2 /[T ]. The kinematic viscosity ν represents a diffusion coefficient for momentum.

Figure 2 . 6 :

 26 Figure 2.6: Schematic flow curves for basic types of time-independent fluids. Reproduced from [46].

Figure 2 . 7 :

 27 Figure 2.7: Dependence of the shear rate on the stress (in semi-log coordinates) for a water-phosphate mixture. The solid triangles represent experimental data.

Figure 2 . 8 :

 28 Figure 2.8: "Phase diagram" for suspension rheology, based solely on a dimensional analysis. Reproduced from [47].

Figure 2 . 9 :

 29 Figure 2.9: Effect of particle size distribution on the maximum packing density for spheres. Polydisperse suspensions have a higher maximum packing density ϕ m ax compared to monodisperse suspensions

Figure 2 . 10 :

 210 Figure 2.10: Conservation of mass for a material volume of a fluid of mass m

Figure 2 . 11 :

 211 Figure 2.11: Pipe and flow parameters

Figure 2 .

 2 Figure 2.12: 2D flow in a pipe

Figure 3 . 1 :

 31 Figure 3.1: Papanastasiou regularization for the Bingham model: stress magnitude (left) and apparent viscosity (right) as functions of the magnitude of the shear rate tensor γ. For increasing m, we recover a closer approximation to the actual Bingham model. Note that logarithmic axes are employed in order to highlight the behavior in the low shear rate limit.

Figure 3 . 2 :

 32 Figure 3.2: Two-dimensional orthogonal grid. N , S , E and W correspond to neighbor cells of cell P; n , s , e and w denote cell P faces; ∆x and ∆y are cell P dimensions in the x and y spatial coordinates; δy n , δy s , δx e and δx w correspond to cell-center to cell-center distances from cell P to neighbor cells.

  H e -D e (∇P ) e(3.31) where D e = (δx e ∆y)/(a p ) e and (∇P ) e represents the pressure gradient in face e along the x-axis. According to Rhie and Chow interpolation, H e and 1/(a P ) e are obtained by a weighted averaging of the values at the neighboring cells as follows:   nb a nb u nb + b u . (3.31) becomes: u e = H e -1 a P e ∆y(p E -p P ) (3.34) where the over-bar refers to linear interpolation. Similar expressions can be formulated for remaining face velocities of cell P . By substituting Eq. (3.32) into Eq. (3.30) and considering Eqs. (3.28) and (3.29), Eq. (3.30) can be reformulated as: u e = [f + e u E + (1 -f + e )u P ] -∆y (p E -p P ) (a P ) e + f + e ∆y (p e -p w ) E (a P ) E + (1 -f + e )∆y (p e -p w ) P (a P ) P (3.35) The term [f + e u E + (1 -f + e )u P ] corresponds to the linear interpolation of cell velocity values. The remaining two terms incorporate the additional correction term that smooths the pressure field, and prevents the undesired checkerboard behavior. Assuming, D e ≈ D e and (D∆p) e ≈ D e ∆p e , and considering Eqs. (3.31) and (3.27), Eq. (3.35) can be re-written as: u e = u e -D e ((∇P ) e -(∇P ) e ) (3.36)

  and (3.39) from (3.25) and (3.26) and considering correction equations (3.40), (3.41) and (3.42) as:

nb a nb u ′ nb and nb a

 nb nb v ′ nb in Eqs (3.43) and (3.44). Thus, velocity correction becomes: u ′ P = -∆y (p ′ e -p ′ w ) P (a P ) P (3.45)

  e = ∆y (a P ) e , d n = ∆x (a P ) n , With known cell-face velocity corrections u * e and v * n , and following the same procedure for the rest of the faces, continuity equation is used to derive an equation for pressure correction p ′ . Eq. (3.1) can be integrated in the cell volume by using Green-Gauss theorem leading to: (ρuA) e -(ρuA) w + (ρvA) n -(ρvA) s = 0 (3.51) Introduction of face velocities from Eqs. (3.49) and (3.50) into the discretized continuity equation, Eq. (3.51) yields the pressure-correction equation:

7 ) 3 . 5 .

 735 The coefficients in the discretized momentum equations and the velocity components u * and v * are computed with a guessed or intermediate pressure field p * using Eqs.(3.38) and (3.39) 4. The interface velocity is calculated with a momentum interpolation method using Eq. (3.37) based on u * , v * and p * to determine the source term of the pressure-correction equation, deduced from the continuity equation The coefficients for the p ′ equation (3.52) are computed 6. The discretized p ′ equation (Eq. (3.52)) is solved 7. The pressure and velocities are corrected using equations (3.43), (3.44), (3.40), (3.41) and (3.42) 8. Other discretized transport equations are solved if necessary 9. Steps 2-8 are repeated until convergence condition is satisfied.

1 . 4 .

 14 Initialize the set of primitive variables and impose boundary conditions 2. At the finer multigrid level, execute few Simple iterations on the system (3.56) 3. Restrict the approximate solutions ϕ * h and the correponding residuals R h Compute the coarse grid correction using Eq. (3.59) 5. Check for convergence. If the solution is converged, prolongate back the velocities and pressures into the fine level 6. The algorithm returns to step 2 to perform few smothing iterations on the fine grid 7. If the solution is not converged yet at the fine multigrid level, repeat steps 2-6.

2 . 3 .a

 23 Few relaxation sweeps on the momentum equations are performed to get approximate solutions u * p and u * p The cell-face velocities are computed according to the momentum interpolation method. For the east interface of the volume cell, we can write: nb u nb + S u + f u

5 .

 5 Pressure equation (3.66) is solved, and the related pressure and velocity corrections are updated through p ′′ 6. The algorithm returns to step 2 to perform additional SIMPLE iterations.

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Geometry of a cavity with moving lid

  case, the boundary conditions are steady velocity uniform profile at the inlet and standard no-slip condition at the walls are implemented. The pressure is derived by extrapolation from the inner nodes. Test cases are solved for Reynolds numbers up to Re = 10 3 , and for Bingham numbers up to 200. The flow configuration and the system of coordinates in both cases are shown in Figs. 3.4-3.5.

Figure 3 . 6 :

 36 Figure 3.6: Sections of velocity along the vertical (a) and horizontal (b) midplanes, Re = 10 3

Figure 3 . 7 :

 37 Figure 3.7: Principal vortex position for various Bn numbers at Re = 10 3 . Reference results retrieved from Vola et al. [73] and Prashant et al. [72].

Figures 3 .

 3 Figures 3.8-3.10 display the flow field as the bingham number increases, for Re = 10 and Re = 103 , respectively. When Re = 10, the streamline contours are symmetric for all Bingham numbers, 2, 10, and 100. However, the unyielded zone is proportional to the yield stress magnitude. At Re = 10 3 , the primary vortex shifts slightly to the right for the case of Bn = 100, and moves further in the same direction for Bn = 10, then towards the center of the cavity for Bn = 2. We see that increasing the Reynolds number moves the vortex first toward the right, and then downwards and left to the center depending on the value of the Bingham number. This trend can be observed for higher Bingham numbers if the Reynolds number is large enough as greater stresses are required in order to make the material flow. In the works of[START_REF] Syrakos | Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method[END_REF][START_REF] Prashant | Direct simulations of spherical particle motion in Bingham liquids[END_REF], similar conclusions were drawn, the flow field at high Bn numbers follow the same trend of any lower Bn number if the Re number is also sufficiently decreased.

  (a) Re = 10 (b) Re = 10 3

Figure 3 . 8 :

 38 Figure 3.8: Streamline contours in Bingham flow for Bn = 2

Figure 3 . 9 : 3 Figure 3 . 10 :

 393310 Figure 3.9: Streamline contours in Bingham flow for Bn = 10

Figure 3 . 11 :

 311 Figure 3.11: Numerical results: the developing velocity profile of a Bingham fluid entering the pipe, Re = 10 2 and Bn = 10.

Figure 3 . 12 :

 312 Figure 3.12: Velocity profiles at x = 0.5L calculated with m = 100 (solid lines), m = 200 (dashed lines), m = 400 (dotted lines), and m = 800 ( red dotted lines). For Re = 10 2 , (a) Bn = 1, and (b) Bn = 10.

Figure 3 . 13 :

 313 Figure 3.13: Apparent viscosity profiles along the vertical section x = 0.5L calculated with m = 100 (solid lines) and m = 800 ( red dotted lines), Re = 10 2 and Bn = 10.

Figure 3 . 14 :

 314 Figure 3.14: The L ∞ norm of u and v residuals as a function of the number of SIMPLE iterations on the fine grid, for Bn = 0, 1, and 10 (Re = 10 2 ). Results are shown for simulations with single grid (SG) in red lines and multigrid (MG) initialization in blue lines.

Figure 4 . 1 :

 41 Figure 4.1: Cumulative particle size distribution of the concentrated phosphate slurry (CPS) samples in terms of the percentage of cumulative volume versus particle size. S1-S6 signify the sample labels.

Figure 4 . 2 :

 42 Figure 4.2: Searle principle: the instrument's motor rotates the measuring bob inside an immovable cup filled with slurry (α cyl.cone = 120 • ).

Figure 4 . 3 :

 43 Figure 4.3: Representative rheograms showing (a) the shear stress as a function of the shear rate of CPS samples at various solids concentrations and a fixed temperature of 25 • C and (b) the apparent viscosity as a function of the shear rate under the same conditions.

.

  implication than the yield stress of the Herschel-Bulkley or Bingham models. In the Robertson-Stiff model, .

Figure 4 . 4 :

 44 Figure 4.4: Comparison of rheological model fits for (a) S1 and (b) S2.

Figure 4 . 5 :

 45 Figure 4.5: Comparison of rheological model fits for (a) S3 and (b) S4.

Figure 4 . 6 :

 46 Figure 4.6: Comparison of rheological model fits for (a) S5 and (b) S6.

Figure 4 . 7 :

 47 Figure 4.7: Apparent viscosity at a constant shear rate of 400 s -1 as a function of the D-value ratio of the large (d L ) and small (d S ) particles.

Figure 4 . 8 :

 48 Figure 4.8: Pipe pressure loss analysis of the Herschel-Bulkley model for S4 and S6.

Figure 4 . 9 :

 49 Figure 4.9: Axial laminar velocity profiles that were calculated based on different rheological models.

Figure 4 . 10 :

 410 Figure 4.10: Flow curves obtained during ramps of increasing then decreasing shear rate, with suspensions reconstituted at different solid concentrations.

Figure

  Figure 4.11: Repeatability tests (up-down ramps) using collected samples (F2). Successive tests on the same sample (left in place) show artifacts evolving with time (migration -sedimentation).

  Figure 4.11: Repeatability tests (up-down ramps) using collected samples (F2). Successive tests on the same sample (left in place) show artifacts evolving with time (migration -sedimentation).

Figure 4 . 12 :

 412 Figure 4.12: Ramp of increasing then decreasing shear stress using collected samples and reconstituted samples at the same solid concentration (55 wt.%).

Figure 4 . 13 :

 413 Figure 4.13: Different tests performed on the two collected phosphate suspensions with a rough surface Couette geometry (described above) and a vane geometry (25 -37 mm). In all cases, increasing and then decreasing stress ramps are applied. The dotted curve has a slope of 2.

Figure 4 . 14 :

 414 Figure 4.14: Ramp of increasing then decreasing stress, and ramp of decreasing shear rate.

Figure 4 . 15 :

 415 Figure 4.15: Comparaison between ramps of increasing stress and shear rate.

Figure 4 . 16 :

 416 Figure 4.16: Apparent flow curves for 4 different samples obtained by imposing a ramp of increasing shear rate, with a Vane geometry of 6 blades, 5 cm height immersed in a cylinder of 11 cm height and a Couette gap of 2.5 mm.

Figure 4 . 17 :

 417 Figure 4.17: Apparent flow curves for different submerged vane heights.

Figure 4 . 18 :

 418 Figure 4.18: Apparent flow curves obtained by imposing a ramp of increasing shear rates, with a Vane geometry 6 blades -5 cm high immersed in a outer cylinder located at different heights from the bottom.

Figure 4 . 19 :

 419 Figure 4.19: Apparent flow curves (vane 5 cm -5 mm from the bottom -2.5 mm couette gap) obtained for different rest times prior to the tests.

Figure 4 . 20 :

 420 Figure 4.20: Apparent flow curves (vane 5 cm -5 mm from the bottom -2.5 mm couette gap) obtained for different ranges shear rate imposed during a increasing 40 s ramp.

Figure 4 . 21 :

 421 Figure 4.21: Apparent flow curves (vane 5 cm -5 mm from the bottom) obtained for different gap widths.

15 )

 15 We denote the matrix by A. The determinant of A is detA =

3 ∇

 3 • ⃗ vI(5.25) where µ is the dynamic viscosity of the Newtonian fluid. For non-Newtonian materials, a different relation is available in which µ ap is the apparent viscosity of the non-Newtonian fluid τ = µ ap γ (5.26)Several rheological models have been developed in the literature to describe the behaviour of these materials. In this work, the models used to describe the behaviour of the phosphate slurry are Herschel-Bulkley (Equation (5.2)), Casson (Equation (5.3)) and the proposed model (Equation (5.4)). ⃗ F cs represents the surface tension force which is modeled as a volumetric force by the continuum surface force (CSF) model of Brackbill et al.[START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. For two-phase system, the value of ⃗ F cs is expressed by Equation (5.27)⃗ F cs = σ ρk∇α s 0.5(ρ w + ρ s )(5.27)where σ represents the surface tension coefficient and k the interface curvature, which is further expressed by Equations (5.28) and (5.29) respectivelyk = ∇ • n (5.28) n = ∇α s |α s | (5.29)The interface between the two phases progresses thanks to the volume fraction Equation (5.30) ∂α s ∂t + ∇ • (⃗ vα s ) = ∇ • (Γ t ∇α s ) (5.30) Γ t = D AB + ν t Sc t (5.31)

Figure 5 . 1 :

 51 Figure 5.1: Comparison of experimental (red circles) and predicted (blue cross) values of shear stress as a function of shear rate for six phosphate slurry test samples.

Figure 5 . 2 :Figure 5 . 3 :Figure 5 . 4 :Figure 5 . 5 :

 52535455 Figure 5.2: Slurry concentration distribution inside an inviscid walls column, Newtonian Model.

Figure 5 . 6 :

 56 Figure 5.6: Flow configuration in the horizontal pipeline.

Figure 5 .

 5 Figure 5.7 illustrates a detailed comparison of the distribution of slurry concentration between the Newtonian and new models, along a horizontal pipe. The Newtonian model is used as a reference solution for validation, as it is a well-established model and there is currently no other reference solution available in literature to validate the new model.

Figure 5 . 7 :

 57 Figure 5.7: Slurry concentration distribution along the pipe for D = 0.9 m, U = 5 m/s, C = 56 wt%, a = 1.56 Pa•s, b = 0.01, Newtonian model (left) and the New model (right).

8

 8

  and 5.9 provide a comprehensive comparison of the distribution of slurry concentration along a horizontal pipe between the Herschel-Bulkley model, the Casson model, and the new model. The new model uses parameters a and b that are specific to the concentration of the slurry. For a concentration of 38.45 wt%, the values of a and b are respectively 0.4 Pa•s and 0.2, and for C = 56 wt%, the values are a = 1.56 Pa•s and b = 0.6. Upon examination of the figures, it can be seen that the new model is able to maintain the slurry batch compact until it exits the pipe, regardless of the concentration of the slurry. In contrast, the Herschel-Bulkley and Casson models exhibit mixing at the water-slurry interface and sedimentation of the slurry batch towards the bottom of the pipe. This suggests that the new model captures the rheological behavior of the phosphate suspension more accurately than the other models, especially for both low and high concentrations. In summary, the new model for phosphate slurry flow visualization has been compared to the Herschel-Bulkley model and Casson model and it has been found that it is able to reproduce well the rheological behavior of phosphate suspensions for both low and high concentrations. The new model can be used to improve predictions and control of the flow of phosphate slurry in industrial settings.

Figure 5 . 8 :Figure 5 . 9 :

 5859 Figure 5.8: Slurry concentration distribution across the pipe for D = 54.9 mm, U = 2 m/s, C = 38.45 wt%, a = 0.4 Pa•s, b = 0.2, Herschel-Bulkley model (left) and the New model (right).

  3.8-3.10 highlights the progressive growth of the unyielded zones. At constant Re number, as the value of the Bingham number increases, the main vortex center gets closer to the upper cavity side and the zone sheared by the upper boundary becomes thinner.

  3

Table 3 .

 3

		1: MG Initialization performance for different Bingham numbers: CPU
	time and speed-up ratio.
	Bingham number Number of cells	CPU time (sec)	Speed-up ratio
			Single Grid MG Initialization
	0	160 x 160	88.2	35.8	2.46
		320 x 320	872.9	143.7	6.07
	1	160 x 160	109.9	32.4	3.39
		320 x 320	990.1	153.4	6.45
	10	160 x 160	81.6	29.8	2.73
		320 x 320	797.3	141.5	5.63
	Table 3.2: MG Initialization performance for different Bingham numbers: Num-
	ber of iterations.
	Bingham number Number of cells	Number of Iterations
			Single Grid MG Initialization
	0	160 x 160		6076	1910
		320 x 320		14272	1964
	1	160 x 160		7643	1670
		320 x 320		16337	1956
	10	160 x 160		5424	1569
		320 x 320		13033	1856

Table 4 .

 4

		1: The elemental composition of CPS samples
	Element		Concentration wt.%
		S1	S2	S3	S4	S5	S6
	Bone phosphate of lime (BPL) 64.5 64.7 64.4 64.8 66.0 65.0
	CO 2	7.6	7.6	7.0	7.7	7.8	7.7
	SiO 2	3.23 3.35 3.62 3.16 3.34 3.71
	MgO	0.80 0.74 0.96 0.83 0.66 0.62
	Al 2 O 3	0.46 0.41 0.46 0.44 0.40 0.47
	Reactive SiO 2	1.43 1.24 1.19 1.29 1.27 1.06

Table 4 .

 4 2: Shear stress data.

	Rotational Speed			Shear Stress
	(rpm)			τ (Pa)
	Measuring Points S1	S2	S3	S4	S5	S6
	776.5	24.47 27.66 47.42 40.55 72.63 66.39
	713.8	23.26 26.48 46.19 39.94 70.31 65.07
	656.2	22.2	25.39 44.44 38.4	67.79 63.13
	603.2	20.47 24.58 42.74 37.07 65.63 61.48
	554.5	19.52 23.66 41.17 35.69 63.71 59.79
	509.7	18.78 22.71 39.77 34.42 61.9	58.22
	468.5	18.01 21.79 38.54 33.25 60.09 56.79
	430.7	17.27 20.91 37.35 32.21 58.39 55.35
	395.9	16.51 20.1	36.24 31.2	56.83 54.02
	364	15.79 19.35 35.13 30.24 55.3	52.66
	334.6	15.09 18.65 34.03 29.32 53.73 51.27
	307.6	14.39 18.02 32.87 28.5	52.21 49.9
	282.7	13.71 17.43 31.74 27.67 50.65 48.4
	259.9	13.09 16.89 30.54 26.84 48.95 46.86
	238.9	12.49 16.35 29.31 26.02 47.27 45.31
	219.6	12.17 15.84 28.12 25.15 45.69 43.69
	201.9	11.55 15.32 26.92 24.28 43.89 42.08
	185.6	11.18 14.8	25.74 23.41 42.34 40.48
	170.6	10.82 14.31 24.59 22.49 40.5	39.03
	156.8	10.39 13.83 23.5	21.63 38.56 37.13
	144.1	9.797 13.32 22.32 20.7	36.79 36.19
	132.5	9.117 12.85 21.34 19.79 34.93 33.85
	121.8	8.646 12.29 20.27 18.92 33.16 32.14
	112	8.036 11.84 19.12 18.09 31.42 30.41
	102.9	7.567 11.42 17.97 17.3	29.63 28.74
	94.61	7.102 10.99 16.9	16.52 27.83 27.02
	86.97	6.586 10.51 15.87 15.63 26.04 25.2

Table 4 .

 4 3: Constant parameters for rheological models.

	Model	Rheological			Test samples
		Parameters S1	S2	S3	S4	S5	S6
	Bingham	τ yB (Pa) η B (Pa•s)	5.54 0.024 0.02 3.81	6.75 0.037 0.03 4.47	16.8 0.051 0.061 18.92
	Casson	τ yc (Pa) η c (Pa.s)	4.42 0.01	2.56 0.012 0.017 0.017 0.019 0.022 5.29 3.27 12.86 14.8
		τ yH (Pa)	1.36	0.89	1.32	0.82	3.21	3.7
	Herschel-Bulkley	K (Pa•s n )	0.75	0.395 1.032 0.6	2.98	3.56
		n	0.51	0.6	0.52	0.58	0.42	0.42
		K (Pa•s n )	0.99	0.5	1.28	0.72	3.87	4.56
	Robertson-Stiff	n . γ 0 (s -1 )	0.47 6.63	0.56 3.84	0.49 1.32	0.56 0.82	0.39 3.22	0.39 3.7
		τ yRS (Pa)	2.40	1.06	1.66	0.65	6.1	7.5

Table 4 .

 4 4: Statistical indicators.

	Model	Statistical			Test samples
		Indicators S1	S2	S3	S4	S5	S6
		R 2	0.97	0.98 0.96	0.95	0.92	0.9
	Bingham	SSE	12.16 9.58 51.11 36.01 210.14 356.49
		RMSE	0.67	0.59 1.37	1.15	2.79	3.63
		R 2	0.99	0.99 0.98	0.98	0.95	0.93
	Casson	SSE	3.94	3.87 29.94 21.22 137.33 252.25
		RMSE	0.38	0.37 1.053 0.88	2.25	3.05
		R 2	0.99	0.99 0.99	0.99	0.98	0.97
	Herschel-Bulkley	SSE	0.48	0.93 6.28	5.41	44.9	93.8
		RMSE	0.13	0.18 0.48	0.44	1.28	1.86
		R 2	0.99	0.99 0.99	0.99	0.98	0.97
	Robertson-Stiff	SSE	0.52	0.97 5.25	4.79	40.4	86.74
		RMSE	0.14	0.19 0.44	0.42	1.22	1.79

of the particles. The results in Fig .4.3 and Table.

4.3 describe the rheological behavior of phosphate slurries with three solids concentrations at a temperature of 25 • C.

  The transport equations of the turbulent kinetic energy k and the specific dissipation rate ω are

			∂k ∂t	+	∂(kv i ) ∂z i	=	1 ρ	∂ ∂z j	Γ k	∂k ∂z j	+ Gk -Y k + S k	(5.33)
	∂ω ∂t	+	∂(ωv i ) ∂z i	=	1 ρ	∂ ∂z j	Γ k ω	∂ω ∂z j	+ Gω + Y ω + D ω + S ω	(5.34)

Table 5 .

 5 1: The two-phase system characteristics.

	Material	Density Kinematic Viscosity Rheological Model (kg/m 3 ) (m 2 /s)
	Water	990	1.00 × 10 -6	Newtonian
	Sludge	1000	6.00 × 10 -6 non-constant	Newtonian Herschel-Bulkley, Casson and New Model

Table 5 .

 5 2: Specification of the multi-phase solver.

	Term	Details
	Name of solver	twoLiquidMixingFoam
	Type of solver	Density-based, segregated solver
	Time dependency	Transient
	Pressure-velocity coupling Pimple
	nCorrectors	3
	nNonOrthogonalCorrector 0

Table 5 .

 5 3: Constant parameters of the new model.

	Model		Phosphate Slurry Samples
	Parameters S1	S2	S3	S4	S5	S6
	a (Pa•s b )	0.93	0.44	0.72	0.57	1.60	1.56
	η ∞ (mPa•s) 16.77 18.33 23.26 23.6	27.45 30.8
	b (-)	0.49	0.59	0.60	0.57	0.56	0.59
	τ c (Pa)	10.20 7.51	15.38 11.03 31.16 37.39
	γ0 [1/s]	12.54 6.22	2.42	22.78 18.46 15.26

Acknowledgements

In our simple example we have illustrated the main concepts of the multigrid methods. In practical CFD calculations, the multigrid transfer process is more sophisticated and different cycles of coarsening and refinement are used with special schedules of restriction and prolongation at different refinement levels. Common choices of multigrid cycles are the so-called V-and W-cycles, which are illustrated in Figure 3.3.

The simple V-cycle shown in Figure 3.3a consists of two legs. The calculation starts at the finer level. Iterations at any level are called relaxation. After a few relaxation sweeps on the finer level, the residuals are restricted to the next coarse level, and after relaxation on that level, the residuals are passed on to the next coarse level and so on until the coarsest level is reached. After final relaxation on the coarsest level, the prolongation steps are performed on the upward leg of the V-cycle until the finer level is reached. In the W-cycle, additional restriction and prolongation sweeps are used at coarser levels to obtain a better reduction of long-wavelength errors. A typical pattern is illustrated in Figure 3.3b.

In general, the multigrid idea requires more than two grids, otherwise its