A New Rheological Model for Phosphate Slurry Flows

Abstract

In this paper, a new rheological model for the flow of phosphate-water suspensions is proposed. The model’s ability to replicate the rheological characteristics of phosphate-water suspensions under different shear rate conditions is evaluated using rheometric tests, and it is found to be in good agreement with experimental data. A comprehensive methodology for obtaining the model parameters is presented. 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 be applied to simulate and monitor phosphate slurry flows in industrial applications.

1. 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 [1], flooding events [2,3], chemical and nuclear reactors [4], fluidized bed [5], combustion and fuel atomization [6,7]. Recently, researchers have devoted a considerable amount of research efforts into developing and implementing sophisticated mathematical techniques, relevant in the whole range of multi-fluid 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. [8] and Hirt et al. [9] 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 [10], dam failure phenomena [11,12,13], bubble rising/nucleation [14,15,16], droplet impact on both dry and wet surfaces [17,18,19,20], metallurgical engineering [21], 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 $\left({\tau }_0\right)$. 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 [22]. The most popular simple viscoplastic model was first proposed by Bingham [23] 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:

$$\left\{\begin{array}{cc}{\displaystyle \dot{\gamma }=0}\hfill & \mathrm{if} \left|\tau \right|\le {\tau }_0\hfill \\ {\displaystyle \tau =\left(\mu +\frac{{\tau }_0}{|\dot{\gamma }|}\right)\dot{\gamma }}\hfill & \mathrm{if} \left|\tau \right|>{\tau }_0\hfill \end{array}\right.$$

(1)

where $\tau$ and $\dot{\gamma }$ are the stress tensor and the shear rate tensor respectively. While the terms ${\tau }_0$, and $\mu$ indicate the yield stress and the plastic viscosity scalar values. $\left|\tau \right|=\sqrt{\frac{1}{2}\tau :\tau }$ and $|\dot{\gamma }|=\sqrt{\frac{1}{2}\dot{\gamma }:\dot{\gamma }}$ denote the second invariants of $\tau$ and $\dot{\gamma }$, respectively. The strain rate tensor $\dot{\gamma }$ is defined as $\dot{\gamma }=\nabla u+{(\nabla u)}^T$, u being the velocity vector. The flow domain of a yield stress fluid is divided into a yielded region $\left(\right|\tau |>{\tau }_0)$ and an unyielded $\left(\right|\tau |\le {\tau }_0)$ region, separated by the yield surfaces where $\left|\tau \right|={\tau }_0$. Two well known generalizations of the Bingham fluid model have been proposed by Herschel–Bulkley [24], and Casson [25], formulated respectively as:

$$\mathrm{Herschel-Bulkley}:\tau =\left(K|\dot{\gamma }{|}^{(n-1)}+\frac{{\tau }_0}{|\dot{\gamma }|}\right)\dot{\gamma }\mathrm{if} \left|\tau \right|>{\tau }_0$$

(2)

$$\mathrm{Casson}:\tau ={\left(\sqrt{\mu |\dot{\gamma }|}+\sqrt{{\tau }_0}\right)}^2\frac{\dot{\gamma }}{|\dot{\gamma }|}\mathrm{if} \left|\tau \right|>{\tau }_0$$

(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 ${\tau }_0=$ and the power-law index to $n=1$, respectively. According to [26], 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–Bulkley (H–B) model to represent the rheological properties [27,28]. Hamza et al. (2018) [29] 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 paper 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.

2. 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.

2.1. 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 paper, 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:

$$\tau =\left\{\begin{array}{cc}a{(\dot{\gamma }+{\dot{\gamma }}_0)}^b\hfill & \mathrm{if} \dot{\gamma }<{\dot{\gamma }}_c\hfill \\ {\tau }_c+\dot{\gamma }{\eta }_{\infty }\hfill & \mathrm{if} {\dot{\gamma }}_c\le \dot{\gamma }\hfill \end{array}\right.$$

(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 ${\dot{\gamma }}_0$ has a different interpretation than the model parameters of both Herschel–Bulkley or Bingham model. In this model, ${\dot{\gamma }}_0$ can be regarded as a correction to the shear rate rather than the shear stress and the yield stress is defined as ${\tau }_y=a{\dot{\gamma }}_0^b$. The linear part of the curve is described using the parameters ${\eta }_{\infty }$ and ${\tau }_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 ${\dot{\gamma }}_c=400{\mathrm{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 $({\dot{\gamma }}_m,{\tau }_m)$ where $m\in \{1,...,M\}$, the parameters of Equation (4) need to be determined. We consider the general case where the model can be formulated as $\tau =f(\theta ,\dot{\gamma })$ where $\theta$ 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:

$$\underset{\theta }{min}\sum _{m=1}^M{[{\tau }_m-f(\theta ,{\dot{\gamma }}_m)]}^2$$

(5)

where the pairs $({\dot{\gamma }}_m,{\tau }_m)$ are observed. Given the observed data pairs $\{({\dot{\gamma }}_1,{\tau }_1),\dots ,({\dot{\gamma }}_M,{\tau }_M)\}$, we may define the error associated to the second branch of Equation (4), $\tau ={\tau }_c+\dot{\gamma }{\eta }_{\infty }$, by:

$$E({\eta }_{\infty },{\tau }_c)=\sum _{m=1}^M{({\tau }_m-({\tau }_c+{\dot{\gamma }}_m{\eta }_{\infty }))}^2$$

(6)

Our goal is to define the values of ${\eta }_{\infty }$ and ${\tau }_c$ that minimize the error function. Following the least-squares method, we should find the values of $({\eta }_{\infty },{\tau }_c)$ such that:

$$\frac{\partial E}{\partial {\eta }_{\infty }}=0,\frac{\partial E}{\partial {\tau }_c}=0$$

(7)

It is to be noted that, in this case, we do not have to worry about boundary points: as $|{\eta }_{\infty }|$ and $|{\tau }_c|$ become large, the fit will clearly deteriorate. Thus we do not need to check on the boundary. Differentiating $E({\eta }_{\infty },{\tau }_c)$ yields:

$$\frac{\partial E}{\partial {\eta }_{\infty }}=\sum _{m=1}^{M}2({\tau }_m-({\tau }_c+{\dot{\gamma }}_m{\eta }_{\infty }))(-{\dot{\gamma }}_m)$$

(8)

$$\frac{\partial E}{\partial {\tau }_c}=\sum _{m=1}^{M}2({\tau }_m-({\tau }_c+{\dot{\gamma }}_m{\eta }_{\infty }))$$

(9)

Setting $\partial E/\partial {\eta }_{\infty }=\partial E/\partial {\tau }_c=0$ and dividing by 2 yields:

$$\sum _{m=1}^M({\tau }_m-({\tau }_c+{\dot{\gamma }}_m{\eta }_{\infty })){\dot{\gamma }}_m=0$$

(10)

$$\sum _{m=1}^M({\tau }_m-({\tau }_c+{\dot{\gamma }}_m{\eta }_{\infty }))=0$$

(11)

We may rewrite these equations as:

$$\left(\sum _{m=1}^M{\dot{\gamma }}_m^2\right){\eta }_{\infty }+\left(\sum _{m=1}^M{\dot{\gamma }}_m\right){\tau }_c=\sum _{m=1}^M{\tau }_m{\dot{\gamma }}_m$$

(12)

$$\left(\sum _{m=1}^M{\dot{\gamma }}_m\right){\eta }_{\infty }+\left(\sum _{m=1}^{M}1\right){\tau }_c=\sum _{m=1}^M{\tau }_m$$

(13)

We have obtained that the values of ${\eta }_{\infty }$ and ${\tau }_c$ which minimize the error (defined in Equation (6)) satisfy the following matrix equation:

$${\displaystyle \left(\begin{array}{cc}{\sum }_{m=1}^M{\dot{\gamma }}_m^2& {\sum }_{m=1}^M{\dot{\gamma }}_m\\ {\sum }_{m=1}^M{\dot{\gamma }}_m& {\sum }_{m=1}^{M}1\end{array}\right)\left(\begin{array}{c}{\eta }_{\infty }\\ {\tau }_c\end{array}\right)=\left(\begin{array}{c}{\sum }_{m=1}^M{\tau }_m{\dot{\gamma }}_m\\ {\sum }_{m=1}^M{\tau }_m\end{array}\right)}$$

(14)

We will show the matrix is invertible, which implies

$${\displaystyle \left(\begin{array}{c}{\eta }_{\infty }\\ {\tau }_c\end{array}\right)={\left(\begin{array}{cc}{\sum }_{m=1}^M{\dot{\gamma }}_m^2& {\sum }_{m=1}^M{\dot{\gamma }}_m\\ {\sum }_{m=1}^M{\dot{\gamma }}_m& {\sum }_{m=1}^{M}1\end{array}\right)}^{-1}\left(\begin{array}{c}{\sum }_{m=1}^M{\tau }_m{\dot{\gamma }}_m\\ {\sum }_{m=1}^M{\tau }_m\end{array}\right)}$$

(15)

We denote the matrix by A. The determinant of A is

$$detA=\sum _{m=1}^M{\dot{\gamma }}_m^2.\sum _{m=1}^{M}1-\sum _{m=1}^M{\dot{\gamma }}_m.\sum _{m=1}^M{\dot{\gamma }}_m$$

(16)

As

$$\overline{{\dot{\gamma }}_m}=\frac{1}{M}\sum _{m=1}^M{\dot{\gamma }}_m$$

(17)

we find that,

$$\begin{array}{cc}\hfill detA& =M\sum _{m=1}^M{\dot{\gamma }}_m^2-{\left(M\overline{{\dot{\gamma }}_m}\right)}^2\hfill \\ \hfill & =M^2\left(\frac{1}{M}\sum _{m=1}^M{\dot{\gamma }}_m^2-{\overline{{\dot{\gamma }}_m}}^2\right)\hfill \end{array}$$

(18)

where the last equality is derived from simple algebra. Hence, while ${\dot{\gamma }}_m$ are not all equal, $detA$ will be nonzero and A will be invertible. Therefore, we note that as long as ${\dot{\gamma }}_m$ are not all equal, the best-fit values of ${\eta }_{\infty }$ and ${\tau }_c$ are achieved by solving a system of linear equations in the Equation (15). As for the first branch of Equation (4), it requires an initial estimate for one of the three parameters. We can use the following procedure for the estimation of these parameters. ${\dot{\gamma }}_0$ can be estimated first using the following equation:

$${\dot{\gamma }}_0=\frac{{\dot{\gamma }}_{min}{\dot{\gamma }}_{max}-{\dot{\gamma }}^{*2}}{2{\dot{\gamma }}^{*}-{\dot{\gamma }}_{min}-{\dot{\gamma }}_{max}}$$

(19)

where ${\dot{\gamma }}^{*}$ is the shear rate corresponding to the geometric mean of the maximum and minimum shear stresses:

$${\tau }^{*}={\left({\tau }_{min}{\tau }_{max}\right)}^{1/2}$$

(20)

Then, taking the logarithm of the first branch of Equation (4):

$$log\left(\tau \right)=log\left(a\right)+blog(\dot{\gamma }+\dot{{\gamma }_0})$$

(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.

2.2. 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 $\alpha$ 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 [30,31] 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 [27,32]. The different phases are represented in the domain by their phase fractions using the VOF [9,33]. 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 Hirt and Nichols (1981) [9] by adding a diffusion term. The phosphate slurry consists of solid particles and water. In this paper, 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.

2.2.1. Governing Equations

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

$${\alpha }_s+{\alpha }_w=1$$

(22)

where ${\alpha }_s$ and ${\alpha }_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

$$\nabla ·\overrightarrow{v}=0$$

(23)

where $\overrightarrow{v}$ represent the velocity. The momentum equation solved for both the phases during the simulation is given by equation

$$\rho \frac{\partial \overrightarrow{v}}{\partial t}+\rho \overrightarrow{v}(\nabla ·\overrightarrow{v})=-\nabla P+\nabla ·\tau +\rho \overrightarrow{g}+{\overrightarrow{F}}_{cs}$$

(24)

where P and g are the pressure and acceleration due to gravity, respectively. $\tau$ represents the shear stress tensor; for Newtonian fluids, it is a linear function of the shear rate given by

$$\tau =\mu \dot{\gamma }=\mu \left[\left(\nabla \overrightarrow{v}+{\left(\overrightarrow{v}\right)}^T\right)-\frac{2}{3}\nabla ·\overrightarrow{v}I\right]$$

(25)

where $\mu$ is the dynamic viscosity of the Newtonian fluid. For non-Newtonian materials, a different relation is available in which ${\mu }_{ap}$ is the apparent viscosity of the non-Newtonian fluid

$$\tau ={\mu }_{ap}\dot{\gamma }$$

(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 (2)), Casson (Equation (3)) and the proposed model (Equation (4)).

${\overrightarrow{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. [34]. For two-phase system, the value of ${\overrightarrow{F}}_{cs}$ is expressed by Equation (27)

$${\overrightarrow{F}}_{cs}=\sigma \frac{\rho k\nabla {\alpha }_s}{0.5({\rho }_w+{\rho }_s)}$$

(27)

where $\sigma$ represents the surface tension coefficient and k the interface curvature, which is further expressed by Equations (28) and (29) respectively

$$k=\nabla ·\widehat{n}$$

(28)

$$\widehat{n}=\frac{\nabla {\alpha }_s}{|{\alpha }_s|}$$

(29)

The interface between the two phases progresses thanks to the volume fraction Equation (30)

$$\frac{\partial {\alpha }_s}{\partial t}+\nabla ·\left(\overrightarrow{v}{\alpha }_s\right)=\nabla ·({\mathsf{\Gamma }}_t\nabla {\alpha }_s)$$

(30)

$${\mathsf{\Gamma }}_t=D_{AB}+\frac{{\nu }_t}{S{c}_t}$$

(31)

where ${\alpha }_s$ is the cell fraction of the slurry phase, ${\mathsf{\Gamma }}_t$ is the effective diffusivity, $D_{AB}$ represents the coefficient of molecular diffusion, ${\nu }_t$ is the turbulent viscosity, and $S{c}_t$ denotes the turbulent Schmidt number. The coefficient of diffusion $D_{AB}$ and the turbulent Schmidt number $S{c}_t$ were set to $1.0\times {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

$$\rho ={\alpha }_s{\rho }_s+{\alpha }_w{\rho }_w={\alpha }_s{\rho }_s+(1-{\alpha }_s){\rho }_w$$

(32)

where ${\rho }_s$ represents the density of the slurry that corresponds to the ${\alpha }_s$ phase fraction, while ${\rho }_w$ is the density that corresponds to the water phase with the ${\alpha }_w$ phase fraction. The calculated density $\rho$ in the Equation (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. The RANS model uses the density $\rho$ in the conservation of momentum equation. When ${\alpha }_s$ is equal to 1, the represented phase is the slurry while the water phase ${\alpha }_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 (30) can be interpreted as an advection–diffusion equation for the ${\alpha }_s$ phase.

2.2.2. Turbulence Equations

However, we should mention that all the constants and wall function of the $k-\omega SST$ model have been derived experimentally by measurements on Newtonian fluids [35]. Moreover, they are not adapted to non–Newtonian fluids. In the present work, we have used $k-\omega SST$ model to simulate the flow of water pushing the slurry batch, thus the turbulence is applied to the water batch.

The transport equations of the turbulent kinetic energy k and the specific dissipation rate $\omega$ are

$$\frac{\partial k}{\partial t}+\frac{\partial \left(k{v}_i\right)}{\partial z_i}=\frac{1}{\rho }\frac{\partial }{\partial z_j}\left[{\mathsf{\Gamma }}_k\frac{\partial k}{\partial z_j}\right]+{\tilde{G}}_k-Y_k+S_k$$

(33)

$$\frac{\partial \omega }{\partial t}+\frac{\partial \left(\omega v_i\right)}{\partial z_i}=\frac{1}{\rho }\frac{\partial }{\partial z_j}\left[{\mathsf{\Gamma }}_k\omega \frac{\partial \omega }{\partial z_j}\right]+{\tilde{G}}_{\omega }+Y_{\omega }+D_{\omega }+S_{\omega }$$

(34)

where ${\tilde{G}}_k$ is the generation of turbulence kinetic energy due to the mean velocity gradients. ${\tilde{G}}_{\omega }$ denotes the generation of the specific dissipation rate. ${\mathsf{\Gamma }}_k$ and ${\mathsf{\Gamma }}_{\omega }$ are the effective diffusivity of k and $\omega$ respectively. $Y_k$ and $Y_{\omega }$ represent the dissipation of k and $\omega$ due to turbulence. $D_{\omega }$ is the cross-diffusion term. $S_k$ and $S_{\omega }$ represent source terms set to be zero in this study. All the above terms are calculated and specified in [36,37,38,39].

3. Results and Discussions

3.1. 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 [40]. 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 1.

The accuracy of the New model is checked by fitting the phosphate slurry data. Figure 1 shows a good agreement with the experimental data for six samples of the phosphate slurry.

3.2. 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 $\alpha$.

3.2.1. 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 [41,42] 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 ${\mu }_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: $\mu =min({\mu }_0,\mu \left(\dot{\gamma }\right))$. The physical characteristics of the two phases used water and sludge are represented in

Table 1. The two-phase system used.

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

. The multiphase solver used twoLiquidMixingFoam is described in

Table 2. Specification of the multiphase solver used.

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

.

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 3. Constant parameters of the new model.

Model	Phosphate Slurry Samples
Parameters	S1	S2	S3	S4	S5	S6
a [Pa·s${}^{\mathrm{b}}$]	0.93	0.44	0.72	0.57	1.60	1.56
${\eta }_{\infty }$ [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
${\tau }_c$ [Pa]	10.20	7.51	15.38	11.03	31.16	37.39
${\dot{\gamma }}_0$ [1/s]	12.54	6.22	2.42	22.78	18.46	15.26

being assigned to the variables ‘a’ and ‘b’ respectively.

An examination of Figure 2, Figure 3, Figure 4 and Figure 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.

3.2.2. Two-Phase Pipe Flow

We recall that the configuration of interest in our work consists of a water-slurry batch flow, represented by Figure 6. A three-dimensional (3D) horizontal pipes are established to investigate the New rheological model of the particulate slurry flow using OpenFOAM simulations. Internal diameters $D=5.49$ cm and $D=90$ cm are used with corresponding pipe lengths $L=3.3$ m and $L=50$ m > 60 D respectively. The lengths of the slurry batch are respectively $l=0.54$ m and $l=8.9$ m.

Table 4. Ranges of parametric values.

Parameters	Ranges	Unit
Pipe diameter	0.054–0.9	m
Pipe Length	3.3–50	m
Solid concentration by mass	38–56	wt%
Water density	1000	kg/m${}^3$
Water viscosity	${10}^{-3}$	Pa·s
Velocity	2–5	m/s
Water Reynolds Number	>${10}^6$	-

represents physical characteristics of the two-phase system water-slurry. According to results in [29], the Herschel–Bulkley model is suitable for concentrations lower than 38.45 wt%, where the phosphate suspension rheological behavior is dilatant. The Casson model is used to model the rheological behavior of the suspension and to calculate the viscosity and yield strength for the concentration range ($46.03$ to 57.27 wt%). Therefore, the new model is compared to the Herschel–Bulkley model for the concentration $C=38.45$ wt% by mass, and to the Casson model in the range of concentrations $C=$ 46.03–57.27 wt%.

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.

Figure 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.

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.

The Figure 8 and Figure 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.

4. Conclusions

A new rheological constitutive relation for modeling the phosphate slurry flows was presented in this study. This model is based on five constant parameters, 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.