On the closure problem of the effective stress in the Eulerian-Eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions

Suspensions of solid particles in liquids appear in several industries, such as pharmaceutical, health-care, food, ink, paint, oil, and mineral (Coussot and Ancey, 1999 and Stickel and Powell, 2005). They are also observed in nature, for instance, in blood, magma flows, lavas, mud flows, and landslides (Mueller et al., 2010). The fluid dynamics of these suspensions are strongly dependent on their rheology, which in turn is dictated by the fluid flow around the particles and by the interactions among the particles (owing, for instance, to lubrication forces, cohesive forces, and friction). Therefore, to design, optimize, and control industrial processing units treating these systems, one has to acquire detailed knowledge about their rheological behavior.

The fluid dynamic behavior of suspensions has been of interest to two fields of research which apparently have not interacted closely over the past decades. These fields are suspension rheology and multiphase fluid dynamics. In the first, researchers focus on the rheological behavior of liquid-particle suspensions, which are regarded as effective fluids; the systems of interest here are those in which all the phases rapidly relax to local dynamical equilibrium and the slip velocity between the phases is small. In the second, researchers focus on the dynamics of multiphase systems in which the slip velocity is significant, as in fluidized beds; in this case, each constituent phase has to be considered separately and constitutive equations for the stress tensor of each phase are necessary. In this second field, the mixture rheology is rarely studied since the mixture may not be regarded and modeled as an effective fluid.

In general, there are mainly three key phenomena observed by researchers in almost all suspension flows. The first is a monotonic increase in the mixture viscosity with solid concentration. Several theories, which date back to Einstein (1906), have been developed to explain this occurrence, which has been detected for all ranges of solid concentration. The most accepted hypothesis maintains that, when the particle concentration increases, the liquid streamlines distort more due to the presence of more particles, and therefore the viscous dissipation in the liquid rises. The second phenomenon is particle migration. This was first recognized in the work of Gadala-Maria and Acrivos (1980). In it, and subsequently in many other experimental investigations, it was observed that in dense suspensions the particles migrate from regions of high to regions of low particle concentration and/or shear. Although the physical interpretation of this phenomenon is still controversial, a few theories have been proposed to explain it. We will discuss this in detail later on. The third phenomenon is the change of rheology, from Newtonian to non-Newtonian, for suspensions in which the interstitial liquid is Newtonian. The particle concentration value at which the change appears is system-dependent, and the cause of the phenomenon is still unclear (Mueller et al., 2010).

To predict the motion of suspensions characterized by the phenomena just described, researchers in both fields of suspension rheology and multiphase fluid dynamics have tried to develop theoretical models based on balance equations for the mixture and for the individual constituent phases. Fundamentally, to understand the interactions between the interstitial fluid and the particles and among the particles themselves, one should solve concurrently the mass, linear momentum, and energy balance equations at every point occupied by the fluid and the Newtonian equations of motion for every particle. This approach, called Eulerian-Lagrangian, is extremely demanding (Pan et al., 2002; Balachandar and Eaton, 2010; Sardina et al., 2012; Oke et al., 2016; Feng et al., 2016; Gualtieri et al., 2017; and Battista et al., 2018) and hence unsuitable for industrial applications. This inadequacy is not just a consequence of how numerically expensive solving these models is; an additional drawback is that they do not provide directly the information that is usually sought by process engineers (e.g., volume fractions and granular temperature and pressure); to extract this, end-users have to filter the results. So, to reduce the computational cost and formulate models able to yield physical quantities of direct interest, researchers have started to treat these systems in terms of mean properties, regarding the constituent phases as continua occupying the entire domain and interpenetrating each other. This modeling approach is known as Eulerian-Eulerian (E-E) or multifluid (Anderson and Jackson, 1967; Jackson, 1997; and Ishii and Hibiki, 2010). For isothermal mixtures of identical particles, these models are based on five balance equations: two sets of mass and linear momentum balance equations (one per phase) plus a pseudointernal energy balance equation for the solid phase. Such equations are subjected to initial and boundary conditions, the latter being assigned on the boundaries of the flow domain.

In some systems, there is rapid relaxation between the liquid and the disperse phase; this means that the time that the particles require to reach local dynamical equilibrium is far shorter than the characteristic time scale of the flow. The ratio of the former, which represents the relaxation time of the particles, to the latter is referred to as Stokes number. When this number is much smaller than unity, a simplified modeling approach, which can be rigorously derived from the multifluid approach, can be adopted (Jackson, 2000 and Ishii and Hibiki, 2010). In this mixture modeling approach, the suspension is treated as one effective fluid, which is referred to as mixture. The model, consequently, comprises one continuity equation and one linear momentum balance equation, which provide the density and velocity fields of the mixture. To obtain the volume fractions of each phase, one needs an additional mass balance equation, which is usually written for the solid phase. Similarly, an additional dynamical equation is required for the velocity fields of the two constituent phases. The latter is a force balance equation, which holds if the fluid and solid velocity fields relax rapidly to local equilibrium. This is normally an algebraic equation, and this is one of the reasons for which the mixture model is simpler than its multifluid counterpart. The local equilibrium assumption behind this approach is usually satisfied for mixtures of micron-sized particles dispersed in common Newtonian liquids, which are the systems of interest in our study. Note that in the literature other names are also adopted for the mixture model, such as drift-flux (Zuber, 1964), algebraic-slip (Pericleous and Drake, 1986), and diffusion (Ungarish, 2013) models.

Both E-E and mixture approaches have been used extensively for modeling liquid-particle suspensions. In the chemical engineering literature, some articles focus on the E-E approach (e.g., Panneerselvam et al., 2007; Chen et al., 2009; Ekambara et al., 2009; Hosseini et al., 2010; Wang et al., 2010; Liu and Barigou, 2013; Qi et al., 2013; and Wang et al., 2013); in these, each phase is normally regarded as Newtonian and the (laminar) viscosity of the solid phase is derived via the kinetic theory of granular flows (e.g., Jung et al., 2006a; 2006b; Jung and Hassanein, 2008; Gidaspow and Huang, 2009; and Yilmaz et al., 2011). Other studies, conversely, are based on the mixture approach (e.g., Silva et al., 2015 and Nayak et al., 2015); these consider the mixture as a generalized Newtonian fluid (the behavior being Newtonian solely in the dilute limit), and the viscosity of the mixture is normally modeled with empirical or semi-empirical constitutive equations (e.g., Wu et al., 2014; 2015; Wen et al., 2015; and Kim et al., 2016). Only few authors have considered both modeling approaches, comparing their performance (see, for instance, Fletcher and Brown, 2009 and Kaushal et al., 2012). Our work belongs to this last category: we propose to investigate both approaches, discuss when they are equivalent, analyze critically the constitutive equations on which they are based, and highlight the limitations that such equations present. Our focus is on the closure problem of the effective stress (of the constituent phases as well as of the mixture), for this poses the most significant challenge. The constitutive equations currently used in many models reported in the literature—in particular, based on the E-E approach—present some shortcomings (both in the dilute and dense limits) that are seldom discussed but which are important. In this article, we address them. We show that, under certain assumptions, all of these continuum approaches should lead to equivalent results; the key point is considering, and modeling accurately, the physical mechanisms responsible for stress generation in suspensions. We thus report and discuss these mechanisms, showing that the deficiency of some closures in replicating some important experimental observations arises because they disregard some essential physical mechanisms, occurring in suspensions, which generate stress.

The article is organized as follows. In Sec. II, we briefly present the multifluid and mixture models. In Sec. III, we first review the mechanisms proposed in the literature as possible causes for stress generation in mixtures; this should guide in deriving and selecting suitable closures for the phasic and mixture effective stress tensors. Afterwards, we discuss some of these closures, commonly used for liquid-solid suspensions, highlighting their limitations in capturing some of the physical mechanisms discussed. As we shall see, some constitutive equations for the effective stress tensors of the individual phases cannot predict a non-Newtonian behavior for the mixture. Furthermore, for dilute suspensions, some closures are not successful at foreseeing the mixture viscosity increase with solid concentration, which is observed experimentally and was predicted theoretically by Einstein (1906) and Batchelor and Green (1972). Postulating that mixtures should behave as generalized Newtonian fluids is not theoretically justified (Massoudi et al., 1999 and Massoudi, 2008), but doing so and employing empirical equations for the mixture viscosity results to be the best way to describe the fluid dynamics of these complex systems. In Sec. IV, we use Computational Fluid Dynamics (CFD) simulations to test the closures previously discussed. In a first case study, we investigate the laminar flow, in a horizontal pipe, of a dilute suspension of equal particles dispersed in a Newtonian liquid, adopting the E-E and mixture models with different closures for the effective stress tensors. To validate the model predictions, we compare the results to experimental data available in the literature. In a second case, we simulate the laminar flow of a dense particle suspension through an abrupt expansion. Our aim is comparing the recirculation lengths of the vortexes generated after the expansion obtained via CFD with those observed experimentally. We show that to predict the recirculation length correctly, one has to account for particle migration owing to shear-induced diffusion in the upstream pipe. Since this is related to the effective stress tensor of the solid phase (not of the mixture), it is crucial that a suitable constitutive equation be used for this term.

In this section, we briefly present the evolution equations of the multifluid and mixture modeling approaches. We focus on suspensions of equal particles dispersed in an isothermal Newtonian liquid. Both models regard the liquid and solid phases as interpenetrating continua, which therefore occupy the entire flow domain at any given time. The mixture model, as said, derives from the multifluid model; the latter is obtained by averaging the equations of motion of the interstitial liquid and of the particles. To this end, statistical, volume, and time averaging methods can be adopted (Mazzei, 2016). The resulting evolution equations are always very similar and present common features. The most noteworthy is that the process of averaging leaves behind a number of indeterminate terms not directly related to the averaged variables but associated with details of the motion at the microscopic length scale. These terms are represented by the effective stress tensors of the phases and by the interaction force exchanged by the phases. A closure problem therefore arises, which normally cannot be solved analytically and has to be overcome by means of empirical expressions. This is the major challenge that this modeling method poses. The same problem, of course, affects the mixture model since the latter is based on the multifluid model, of which it is an approximation.

The relations for the effective stress tensors and fluid-particle interaction force that the averaging schemes yield have no direct use but are extremely important nevertheless; this is because they indicate the origins of these terms and allow understanding their physical meaning. These expressions, furthermore, can be used to derive closure equations via more fundamental modeling approaches, such as the Eulerian-Lagrangian. The mathematical form of these expressions depends on the averaging method used. In this article, we adopt the averaging scheme proposed by Jackson (1997), which is based on “soft” volume averages. Also, the average variables associated with the solid phase are derived by treating the particles as rigid bodies and by averaging their overall mass and linear and angular momenta (Zhang and Prosperetti, 1994; 1997). Henceforth, when we refer to the mean value of a (point) variable, we intend its volume-averaged value; depending on the context, the volume average can be of the fluid or solid type.

For liquid-solid isothermal suspensions of monodisperse (that is, of identical) particles, the multifluid model involves a mass balance equation and a linear momentum balance equation per phase, plus a pseudointernal energy balance equation for the solid phase (an energy balance equation for the fluid is not needed, being the medium isothermal). For the liquid, the equations are

where ε_e and ρ_e are the volume fraction and density of the fluid, respectively, $⟨u⟩e$ is the mean fluid velocity, $⟨S⟩e$ is the fluid effective stress tensor, g is the gravitational field, and $⟨f⟩p$ is the force exerted by the fluid on the particles per unit volume of suspension. This term, representing the interaction force between the phases, has four main contributors: the buoyancy, drag, virtual mass, and lift forces. For details about these forces, we refer to the literature (Jackson, 2000).

Equation (2.2) is unclosed because one does not usually know how the terms $⟨S⟩e$ and $⟨f⟩p$ relate to the average variables characterizing the flow. For the fluid-particle interaction force, the relation given by the averaging procedure is the following:

where σ_e is the point stress tensor of the fluid, ψ is the weighting function employed in the volume-averaging procedure, z_r is the (time-dependent) position of the center of the generic particle r, while k is the outward unit normal to the surface ∂Λ_r bounding particle r. The relation above is clearly unclosed since to evaluate its right-hand side, one would need to know the position of the center of each particle at any given time plus the point stress tensor of the fluid as a function of the time and space coordinates. Only Eulerian-Lagrangian simulations would be able to give this information.

The expression of the effective stress tensor of the fluid is more complex (and perhaps less intuitive) and reads as follows:

where

Here, $⟨σ⟩e$ is the mean value of the point stress tensor of the fluid, n is the particle number density, $⟨u^u^⟩e$ is the mean value of the dyadic product of the fluid velocity fluctuations, while a is the particle radius. Also these equations, of course, are unclosed.

The first term on the right-hand side of Eq. (2.4) is present because, before the averaging is carried out, the liquid phase is already a fluid and is therefore already endowed with a point stress tensor (as we shall see, no analogous term features in the averaged linear momentum balance equation of the solid phase). The last term on the right-hand side of Eq. (2.4) is a Reynolds stress type of contribution, not necessarily turbulence-related. The other two terms which contribute to the effective stress tensor of the fluid are related to the fluid dynamic interaction between the fluid and the particles (through the traction force exerted by the fluid on the surface of the particles). These terms, therefore, are not only related to the fluid phase but also directly involve the particles; for this reason, their combined contribution is sometimes referred to as particle-presence stress (Hwang and Shen, 1989). Note that we did not choose to ascribe the latter to the effective stress tensor of the fluid; this is a direct result of the volume averaging method adopted, which is just a formal mathematical procedure. The particle-presence stress, therefore, is not part of the effective stress tensor of the solid phase. It is also worth pointing out that this type of stress plays an important role and should not be neglected. For instance, as we shall discuss in more detail subsequently, the correction derived by Einstein (1906) for the mixture viscosity of dilute suspensions arises from the particle-presence stress tensor.

The particle phase is governed by the averaged mass and linear momentum balance equations. These can be derived via two averaging methods, which are often referred to as solid averaging and particle averaging (Zhang and Prosperetti, 1994 and Jackson, 1997). In general, solid and particle averages of equal point variables have different values; in particular, it is

in which ε_s is the volume fraction of solid, V_p is the particle volume, $⟨u⟩s$ and $⟨u⟩p$ are the solid and particle averages of the particle velocities, $⟨ω⟩p$ is the particle average of the particle angular velocities, while L is the length scale characterizing the gradients of the mean fields. If we assume that this length scale is much larger than the particle radius (this is normally the case; furthermore, if this were not true, the values of the averaged variables would depend on the form of the weighting function adopted), the last terms on the right-hand side of the equations above can be neglected. Moreover, if $⟨ω⟩p$ has the same order of magnitude as $∂x×⟨u⟩e$⁠, a common situation—in particular, for liquid-solid mixtures, the second term on the right-hand side of Eq. (2.8) can also be neglected, being O(a²/L²). Therefore, at this level of approximation, ε_s = V_pn and the solid and particle averages of the particle velocities, as well as the equations governing their evolution, coincide. Then, the averaged mass and linear momentum balance equations read

where ρ_s is the point density of the solid material, $⟨u⟩s$ is the mean velocity of the particles, while $⟨S⟩s$ is the solid effective stress tensor, defined as follows:

where

in which f_rs is the force exerted by particle r on the generic particle s and k_rs is the outward unit normal to the surface of particle r at its point of contact with particle s.

The first two terms on the right-hand side of Eq. (2.11), taken together, represent the stress contribution due to direct particle contacts; we may refer to it as particle-contact stress. Direct particle contacts can be of two kinds: (nearly) instantaneous contacts forming during particle collisions and enduring contacts that establish when mixtures become very dense. The first kind gives rise to the collisional stress, which captures the linear momentum transfer over the distance 2a present between the centers of two colliding particles. This physical phenomenon becomes important in dense mixtures, in which the total volume occupied by the particles is not negligible compared to the volume of the system containing them. The second kind of direct particle contacts brings about the so-called frictional stress, which, as said, is important only in very dense mixtures. The last contributor to the solid effective stress, which arises from the Reynolds decomposition of the advection term in the averaged dynamical equation, is called kinetic stress and is related to the solid velocity fluctuations $u^s$ (referred to as peculiar velocities in kinetic theory, Chapman et al., 1990).

The final balance equation that needs to be considered for the solid phase is the transport equation for the pseudointernal energy of the particles; this reads

Here U_s ≡ 3/2Θ_s is the pseudointernal energy per unit mass of solid, Θ_s is the granular temperature, and q_s is the pseudothermal heat flux. Equation (2.14) differs from the customary internal energy balance equation (Bird et al., 2007 and Deen, 2012) because of two sink terms S_c and S_v representing pseudointernal energy degradation by inelastic collisions and by viscous resistance to particle motion, respectively, and a source term G_d related to the generation of particle velocity fluctuations (and therefore of pseudointernal energy). Several physical phenomena contribute to this source term, such as Brownian motion, mean fluid dynamic shear, fluctuating fluid-particle forces, and turbulence, but rarely all of them coexist. In this work, in particular, we assume that the particles are not Brownian and that the flow regime is laminar; hence, the respective generation terms are disregarded. We shall discuss the other two source terms (related to mean fluid dynamic shear and fluctuating fluid-particle forces) in Sec. III. For briefness, we do not report the expressions that the averaging scheme yields for the terms featuring in Eq. (2.14); the reader can refer, for instance, to Gidaspow (1994).

To derive the averaged mass balance equation for the mixture, we need to add the continuity equations for the two phases, that is, Eqs. (2.1) and (2.9). Doing so yields

where

To derive the linear momentum balance equation for the suspension, we have to sum the dynamical equations for the two phases. This yields

where

In the equation above, ω_e and ω_s are the mass fractions of the fluid and solid phases, respectively. To derive Eq. (2.17) in the form reported, we used the relations

$⟨S⟩d$ is referred to as diffusion stress tensor (Manninen et al., 1996) and arises because each phase moves at a velocity that differs from that of the mixture.

In addition to the continuity equation for the mixture, another mass balance equation is required to track the volume fractions of the phases. Usually, this is written for the solid phase. Expressing the convective flux in Eq. (2.9) in terms of mixture velocity and using Eq. (2.19), we obtain

The term in square brackets on the right-hand side of this equation arises because the solid does not move at the same velocity as the mixture. Hence, it may be interpreted as a diffusive flux.

As observed, diffusive fluxes (of mass and linear momentum) arise in Eqs. (2.17) and (2.20). These terms are unclosed because they involve the fluid-particle slip velocity instead of the mixture velocity. The problem of closure, however, can be overcome if the fluid and solid velocity fields relax rapidly to local equilibrium. In this case, an approximate equation can be obtained for the slip velocity. In the literature, different derivations of this equation have been reported (see, for instance, Manninen et al., 1996 and Jackson, 2000). In  Appendix A of this article, we obtain the equation employing a perturbation method so that the relative importance of the various terms contributing to the slip velocity is clearly revealed. The equation may be expressed in different forms, which, at the order of approximation at which the equation holds, are all equivalent; in this article, we favor the following expression:

which clearly relates particle migration to the solid effective stress tensor; notice that in mixtures of neutrally buoyant particles, the term involving $⟨S⟩s$ is the only one present. Above, β is the drag force coefficient and $Dt⟨u⟩m$ represents the material derivative of the mixture velocity field.

As discussed in Sec. II, the effective stress tensors of the fluid and solid phases, and the force of interaction between the phases, are the unclosed terms in the averaged dynamical equations of the multifluid and mixture models. The closures for the various contributions to the fluid-particle interaction force have been discussed extensively in the literature (e.g., Jackson, 2000; Mazzei and Lettieri, 2007; and Marchisio and Fox, 2013); so, they are not considered in this work. In this section, we focus on the closures for the effective stress tensors and on the shortcomings that some of these present in predicting the rheology of liquid-solid mixtures. Before doing so, however, we briefly review the physical mechanisms which generate stress in the constituent phases; this insight is necessary in the discussion that then follows.

As Eq. (2.4) reveals, the effective stress tensor of the liquid phase has three contributors: the mean action of the point fluid dynamic stress, the particle-presence stress, and the stress generated by fluid velocity fluctuations. As discussed by Joseph et al. (1990), the first term is closed and is equal to

where ⟨p⟩_e is the mean value of the fluid pressure, I is the unit tensor, and $⟨γ̇⟩v$ is twice the deformation rate tensor of the following velocity field:

Note that this stress contribution does not involve only the mean velocity of the fluid since $⟨u⟩v$ is related to the mean velocity of both phases. However, for mixtures of nearly neutrally buoyant particles where the two velocity fields rapidly relax, the mean velocities of the phases are equal at leading order in the Stokes number (see  Appendix A), and so, at this approximation level, $⟨u⟩v$ and $⟨u⟩e$ coincide. Notice that these are precisely the systems for which the mixture modeling approach is usually used. If the particles are not nearly neutrally buoyant, the mixture model is still valid provided the magnitude of the mean settling velocity is far less than that of the phasic mean velocities; this is because—if this condition is met—the mean velocities of the phases are still equal at leading order in the Stokes number. In any other case, the mixture modeling approach should not be employed because the phasic mean velocities differ considerably, and therefore an effective viscosity for the mixture can no longer be introduced (Massoudi, 2008).

The particle-presence stress arises due to the viscous dissipation of the (point) kinetic energy in the liquid containing the particles. At any solid volume fraction, these distort the fluid streamlines and generate larger velocity gradients, which increase viscous dissipation and in turn the effective viscosity of the suspension. In dense mixtures, lubrication films form between a significant number of particles; the viscous dissipation that ensues is also captured by the particle-presence stress. This last contribution makes the effective viscosity of the suspension diverge when the solid volume fraction nears its maximum value; however, the rate at which this occurs is less than that observed experimentally (Marrucci and Denn, 1985), and so the dominant cause of the viscosity divergence must be another (for instance, the enduring contacts between particles forming close to maximum packing, which yield the frictional part of the solid effective stress).

As Eq. (2.4) reveals, the particle-presence stress is part of the effective stress of the fluid phase; as stressed, this is a result that the averaging scheme yields. Notice that the same conclusion was also reached by Zhang and Prosperetti (1997) via ensemble averaging. In the literature, however, several researchers have incorrectly ascribed this contribution to the effective stress of the solid phase (see, for instance, Nott and Brady, 1994; Morris and Boulay, 1999; and Palma et al., 2016). Interestingly, Buyevich (1999) allocated the term related to the viscous dissipation due to the fluid streamline distortion to the fluid effective stress and that related to the viscous dissipation in the lubrication films to the solid effective stress. Nevertheless, his repartition—as well as that used in many other studies—was not rigorously justified. For more details about this aspect, we refer to Nott et al. (2011), whose analysis is quite thorough and rigorous.

The last contributor to the fluid effective stress is related to the fluid velocity fluctuations. In suspensions where the only element of randomness is given by the particle positions (for instance, in mixtures of particles dispersed in liquids flowing in laminar conditions), these fluctuations generate because on the surface of the particles the fluid is forced to satisfy no-slip boundary conditions. The difference between the point velocity of the fluid and its locally averaged value is expected to be of the order of the point relative velocity between the fluid and the particles. Thus, in systems in which this relative velocity is small, the stress arising from the fluid velocity fluctuations is negligible; this is the case, in particular, for mixtures of nearly neutrally buoyant particles in which the phasic velocity fields rapidly relax to equilibrium.

Moving on to the solid phase, its effective stress has solely two sources: particle velocity fluctuations and particle contacts. The latter can be (nearly) instantaneous collisions or enduring contacts that establish when the suspension becomes very dense. In liquid-solid mixtures of particles of moderate density, collisions have negligible effect because the lubrication films forming between approaching particles are robust, preventing the particles from colliding. Consequently, far from the jamming state, the collisional stress is expected to be zero, the main source of stress being the particle velocity fluctuations. By contrast, in denser suspensions, and in particular near the jamming state, the dominant stress source is the frictional stress term, which is a part of the particle-contact stress (Jackson, 2000; Mazzei et al., 2010; Oke et al., 2015; and Townsend and Wilson, 2017). Near the jamming state, mixtures behave like dense, slowly deforming granular materials. In this article, we will primarily consider systems far from the jamming state.

As stated before, particle velocity fluctuations have various origins: Brownian motion, turbulence, mean fluid dynamic shear, and fluctuating fluid-particle forces. We focus on systems of non-colloidal particles, for which Brownian fluctuations are negligible, and mixtures for which turbulent effects are absent or effectively suppressed, such as in concentrated suspensions (Buyevich, 1999). Accordingly, the only sources of velocity fluctuations are the mean fluid dynamic shear (which yields the so-called shear-induced fluctuations) and the fluctuating fluid-particle forces (which yield the so-called concentrational or pseudo-turbulent fluctuations). The latter occur because in any mixture there inevitably appear random fluctuations in particle concentration that yield fluctuations in the fluid-particle interaction force, particularly in the drag force, and in the effective weight of the suspension, which in turn induce fluctuations in the velocities of the particles (and of the fluid). This source of stress, nevertheless, is important only in systems in which the fluid-particle relative velocity is considerable, such as gas-fluidized beds or mixtures of dense particles sedimenting in liquids (here particles are far from being nearly neutrally buoyant). In liquid-solid mixtures of nearly neutrally buoyant particles in which the velocity fields rapidly relax, the relative velocity is zero at leading order in the Stokes number, and so this source of stress is negligible. Shear-induced particle velocity fluctuations, therefore, remain the only significant source of solid stress in these systems.

Shear-induced particle velocity fluctuations are generated by particle interactions occurring in shear flow when adjacent mixture layers move past one another. These interactions yield random particle displacements of the order of the particle size, which result in a (non-Newtonian) particulate stress. One important effect of this specific stress is the so-called shear-induced particle migration, which we briefly discussed in Sec. I. This phenomenon was first observed and investigated by Gadala-Maria and Acrivos (1980), who reported that in dense suspensions particles migrate from regions of high to regions of low mean shear. As Eq. (2.21) shows, in mixtures of nearly neutrally buoyant particles (where the solid and mixture densities are almost identical), this diffusive transport of particles originates from gradients in the effective stress tensor of the solid, which, in turn, arises from the shear-induced fluctuations in particle velocity.

For the systems of interest in this work (suspensions of nearly neutrally buoyant particles dispersed in liquids where the phasic velocity fields rapidly relax to equilibrium), the closure problem of the fluid effective stress tensor requires deriving a constitutive equation for the particle-presence stress tensor [the term involving the tensors $⟨A⟩p$ and $⟨B⟩p$ in Eq. (2.4)]. This has been performed analytically for very dilute suspensions of Stokesian spherical particles (i.e., particles for which the Reynolds number is vanishingly small) by Jackson (1997) and Zhang and Prosperetti (1997). Even if the approaches they employed are quite different, they found equivalent results (at the level of approximation considered); these are

where $⟨γ̇⟩e$ is twice the deformation rate tensor of the fluid mean velocity field and ϵ is the Levi-Civita tensor. As reported, these relations are accurate to O(ε_s). For the systems of interest in this paper, the rotational slip velocity between the phases is very small (that is to say, $|∂x×⟨u⟩e−⟨ω⟩p|$ is very small); so, the third term on the right-hand side of Eq. (3.3) can be neglected. One can also disregard the contribution of the tensor $⟨B⟩p$ since we are dealing with suspensions in which also the linear slip velocity between the phases is very small. In particular, at leading order in the Stokes number, it is

So, at the level of approximation adopted, in Eq. (3.1), we can replace $⟨γ̇⟩v$ with $⟨γ̇⟩e$⁠. Now, using the relations above along with Eqs. (2.4) and (3.1) (more details can be found in  Appendix B), we find

In the limit of extreme dilution, shear-induced fluctuations are absent, and so the solid effective stress tensor vanishes. Also the diffusion stress tensor is negligible because at leading order in the Stokes number the slip velocity between the phases is zero. Accordingly, the effective stress tensor of the fluid coincides with that of the mixture, defined in Eq. (2.18). Hence, we obtain

recovering the result obtained by Einstein (1906) for mixtures of non-Brownian, neutrally buoyant, spherical particles dispersed in Newtonian fluids in laminar flow. Equation (3.7) is closed because, being the interstitial fluid incompressible, $⟨p⟩e$ requires no constitutive equation (the unknown fields associated with the fluid phase are the mean velocity and the mean pressure, the fluid volume fraction being equal to 1 − ε_s).

The result just obtained holds for very dilute mixtures, in which the solid volume fraction normally does not exceed 0.05. In these conditions, the fluid dynamic interactions between spherical particles are negligible. Non-spherical particles start interacting sooner, and so for them, the validity range of Eq. (3.7) is narrower. For semi-dilute mixtures, in which fluid dynamic particle interactions are no longer negligible, the solid effective stress tensor can still be taken to be vanishingly small; so, the fluid effective stress tensor is still equal to the effective stress tensor of the mixture. The latter is given by Eq. (3.7), but the effective viscosity of the mixture is no longer a linear function of ε_s. Usually, the following relation is used:

with A and B taken within the ranges [1.5, 5] and [7.35, 14.1], respectively (Mueller et al., 2010). Batchelor and Green (1972), in particular, derived an expression for η_m analytically, finding A to be consistent with the Einstein coefficient and B = 7.6. This expression, however, is valid solely for steady pure straining motions of the suspension. For general flows, the values of A and B are usually obtained empirically.

For spherical particles, Eq. (3.8) starts giving poor predictions when ε_s exceeds about 0.25. This threshold value can be far lower for non-spherical particles and is system-dependent. For denser mixtures, the problem of closure of the effective stress tensors of the phases and of the mixture is significantly more complex; this is because the mixture starts displaying non-Newtonian behavior (often shear-thinning at moderate shear rates and shear-thickening at large shear rates; Jeffrey and Acrivos, 1976), and the cause of such a rheological change is unclear. However, for non-Brownian particles, several researchers (see, for instance, Gillissen and Wilson, 2018) believe that the origin of this behavior is related to the fluid dynamics within the lubrication films between the particles, which is captured by the fluid effective stress tensor. Various mechanisms have been suggested to explain the rheological change. Mueller et al. (2010), for instance, state that the shear-thinning behavior stems from the large viscous dissipation rates of point kinetic energy into point internal energy taking place within the lubrication films present among the particles; these lead to localized overheating of the liquid, which in turn results into a localized drop in liquid viscosity. The overheating is larger when the average shear rate of the mixture is larger, and this explains the shear-thinning rheological behavior. This is a compelling suggestion, but it still has not been verified (neither experimentally nor numerically). Another mechanism, again related to the presence of lubrication films, has been recently put forward by Vázquez-Quesada et al. (2016). They claim that, even if at moderate shear rates the interstitial liquid is Newtonian, at the extremely large local shear rates present in the lubrication films, the liquid shear-thins, leading to the shear-thinning behavior of the mixture. These researchers validated this idea numerically, via Stokesian dynamics simulations.

Both mechanisms presented are based on the notion that the suspension effective viscosity is dominated by the kinetic energy dissipation process taking place within the lubrication films. This conviction is shared by many research groups and has been used often to derive closures for the mixture effective viscosity (e.g., Frankel and Acrivos, 1967; Goddard, 1977; Jarzebski, 1981; and Buyevich, 1999). If we accept this point of view, then we may conclude that for mixtures in the semi-dense regime (that is, in a range of solid volume fraction for which particle enduring contacts are either absent or negligible in effect) the effective stress tensors of the fluid and of the mixture must be nearly equal; in particular, the effective stress tensor of the fluid has to obey the constitutive equation that defines shear-thinning fluids. This does not mean that the solid effective stress tensor can be taken to be zero; experimental evidence shows that in the semi-dense regime particle migration does take place, which implies that the solid effective stress tensor cannot be vanishingly small. However, we may expect this stress to contribute negligibly to the effective stress of the mixture (of course, for very dense suspensions, this is no longer true because the frictional part of the solid effective stress is expected to play a dominant role). So, in the semi-dense regime, we suggest this closure,

where

This is the magnitude of $⟨γ̇⟩e$ (in a simple one-dimensional shear flow, $|⟨γ̇⟩e|$ is equal to the absolute value of the shear rate). Furthermore, in Eq. (3.9), K is the consistency and n is the power-law index; their values may be obtained experimentally in simple viscometric flows and are functions of the volume fraction of solid. Of course, η_m can be closed with other relations, such as that of Herschel and Bulkley (1926).

In the multifluid modeling approach, in many articles, the effective stress tensor of the fluid is assumed to be Newtonian; the constitutive equation used therefore reads

where η_e is the fluid effective viscosity. In some studies, the flow regime is turbulent and so an additional term is added to the fluid effective stress tensor; here we do not report it because in this article we are focusing on laminar flow conditions.

If Eq. (3.11) is employed, the closure problem reduces to finding a suitable constitutive expression for η_e. In the literature (refer, for instance, to Jung and Hassanein, 2008; Chen et al., 2009; Gidaspow and Huang, 2009; Fan et al., 2010; Yilmaz et al., 2011; Kaushal et al., 2012; Wang et al., 2013; and Ofei and Ismail, 2016), the following relation is frequently adopted:

We believe that this closure is incorrect and should be avoided. As we shall discuss in Sec. III C, in the multifluid modeling framework, the solid effective stress tensor is often closed through constitutive equations based on the kinetic theory of granular flows. For liquid-solid mixtures, these granular models yield a stress tensor that is considerably small, and so the effective stress tensors of the fluid and of the mixture are almost identical. Based on our previous discussion, this appears to be correct. However, Eq. (3.12) must be rejected because it predicts a decrease in viscosity as the solid concentration rises; also, in the semi-dense regime, it does not predict shear-thinning behavior. When modeling turbulent mixtures, usually a turbulent viscosity is added on the right-hand side of Eq. (3.12). This additional term is dominant, and this explains why, in articles where the laminar viscosity is expressed via Eq. (3.12), good numerical results can be obtained. But this does not justify the use of an incorrect model for the laminar part of the fluid effective stress tensor.

In some other articles (e.g., Jung and Hassanein, 2008 and Yilmaz et al., 2011), Eq. (3.12) is employed, but the effective stress tensor of the solid is assumed to be essentially equal to that of the mixture (these studies do not resort to granular kinetic theory to overcome the solid stress closure problem). Hence, the use of Eq. (3.12) does not significantly affect the mixture dynamics because the effective stress of the fluid is not the dominant contributor to the mixture effective stress. Nonetheless, in light of our previous considerations, this modeling choice appears to be incorrect. In the semi-dense region, if the non-Newtonian behavior of the mixture stems from kinetic energy dissipation in lubrication films, then, since this effect is captured by the particle-presence stress, the fluid effective stress must play the dominant role.

As previously discussed, in the semi-dense region, the solid effective stress is of kinetic nature and originates from particle velocity fluctuations. Thus, to derive a constitutive equation for the effective stress tensor of the solid, one has to identify the causes of these fluctuations and then correctly model their effect. In Sec. III A, we have concluded that, for the system of interest here, the dominant cause is the mean fluid dynamic shear, and so this is the mechanism on which to focus.

In the literature on fluid-particle flows, an approach that has been frequently adopted to tackle the closure problem of the particulate stress is that of kinetic theory for granular flows (Gidaspow, 1994; Nott and Brady, 1994; and Berzi and Fraccarollo, 2016). This modeling framework is based on a generalization of the Boltzmann equation (Cercignani, 2012), which was originally derived for molecular systems. For them, the stress is Newtonian, the fluid viscosity is a function of the temperature, which is related to the molecule velocity fluctuations, and the temperature field is governed by an internal energy balance equation. This equation has one source of internal energy, owing to the kinetic energy degradation induced by viscous dissipation, and no sinks of internal energy. Similarly, for fluid-particle mixtures, the stress is Newtonian, the solid-phase viscosity depends on a granular temperature, which is related to the velocity fluctuations of the particles, and the granular temperature field is governed by a (pseudo)internal energy balance equation. This equation was presented in Sec. II A; see Eq. (2.14). As we said, in addition to the usual source term mentioned for molecular systems, the internal energy of the particle phase has two sinks and one source. The sink S_c accounts for the internal energy degradation due to inelastic collisions; for the systems we are considering in this article, this term is negligible. The sink S_v accounts for the internal energy degradation due to viscous resistance to particle motion; this term cannot be disregarded. Finally, the source term G_d accounts for the generation of particle velocity fluctuations owing to interactions between the fluid and the particles. The principal contributors to this term are the fluctuations induced by the mean fluid dynamic shear and the concentrational fluctuations. For the systems we are considering, the latter is negligible, but the former plays a crucial role and so has to be correctly modeled.

Unfortunately, deriving an expression for G_d is very complex. This has been achieved for the part related to concentrational fluctuations (Koch, 1990; Buyevich, 1994; and Buyevich and Kapbasov, 1994), but doing so for the part related to shear-induced fluctuations is an open challenge. Accordingly, some researchers have opted to neglect this term; the model developed by Gidaspow (1994), which is frequently used in the literature, is a prominent example. But this model was developed primarily for gas-fluidized beds, in which shear-induced fluctuations play a minor role; for liquid-particle suspensions, the same approximation cannot be accepted. In spite of this, several studies in the literature have employed Gidaspow’s model (see, for instance, Chen et al., 2009; Gidaspow and Huang, 2009; Kaushal et al., 2012; Wang et al., 2013; and Ofei and Ismail, 2016). This leads to a considerable underestimation of the solid effective stress tensor because the dominant source of particle velocity fluctuations is not accounted for, which in turn results in an underestimation of particle migration. In some flows (dominated by convection), this may not pose a problem, but in others, one of which we consider in Sec. IV B, it affects significantly the model predictions.

Because of the complexity posed by the generation term G_d, some researchers have chosen to model the solid effective stress tensor by relying on simple, heuristic arguments. In this case, the pseudointernal energy balance equation is no longer included in the model. To the best of our knowledge, two constitutive equations for the shear-induced effective stress tensor of the solid have been advanced in the literature. The first is that of Buyevich (1996), which reads

Here C is a constant, $⟨γ̇⟩s$ is twice the deformation rate tensor of the solid mean velocity field, and $|⟨γ̇⟩s|$ is its magnitude, defined with a relation analogous to Eq. (3.10). Furthermore, it is

where $εs⋆$ represents the maximum value that the solid volume fraction can attain when the suspension reaches packing conditions. Buyevich obtained this equation by assuming that the shear-induced particle fluctuations scale with $|⟨γ̇⟩s|a$ (or equivalently with $|⟨γ̇⟩m|a$⁠, for at leading order solid and mixture velocities are equal), by invoking the principle of material frame-indifference (Truesdell, 1977) to identify the permissible form of the constitutive equation and by using heuristic arguments concerning the pure shear flow of a liquid-particle mixture. His closure predicted the fully developed solid volume fraction profile in various viscometric flows with satisfactory agreement with experimental data (Buyevich, 1996 and Buyevich and Kapbsov, 1999), but it has never been used to predict the time evolution of solid volume fraction profiles.

The second constitutive equation for the solid effective stress tensor was derived heuristically by Morris and Boulay (1999) in a study aimed at investigating the influence of normal stresses on shear-induced particle migration. This equation reads

The first term on the right-hand side represents the viscous shear-stress part of the solid stress tensor, which is assumed to be Newtonian. η_s is a function of the solid volume fraction representing the particle contribution to the suspension shear viscosity, and $⟨γ̇⟩m$ is twice the rate of strain tensor of the mixture, with $|⟨γ̇⟩m|$ being its magnitude [defined with an equation similar to Eq. (3.10)]. The second term on the right-hand side accounts for the presence of viscous normal stresses; η_n is a function of the solid volume fraction referred to as normal stress viscosity, whilst Q is a diagonal, anisotropic material property tensor whose components are constant fitting parameters. If these parameters are properly tuned, Eq. (3.15) can predict well steady-state solid volume fraction profiles in several viscometric flows and transient solid volume fraction profiles in wide-gap Couette flows (Morris and Boulay, 1999). This equation, however, presents a considerable limitation: even if written in general tensorial form, it was derived—and is valid only—for viscometric flows (the anisotropy of the normal stress tensor arises inasmuch as in this kind of flows the normal stresses are different in the flow, gradient, and vorticity directions). Hence, in contrast with Eq. (3.13), the validity of Eq. (3.15) is quite limited; the equation, in particular, cannot be employed to model generic three-dimensional flows. Miller et al. (2009) generalized it to two-dimensional flows, but doing so requires the introduction of additional fitting parameters and does not render the closure of general validity. Thus, in this article, we adopt the constitutive equation of Buyevich (1996), which can be used for the complex flow investigated in Sec. IV B.

Finally, we have to mention that in many articles the effective stress tensor of the solid phase is assumed to be Newtonian; therefore, it is closed as follows:

Here $⟨p⟩s$ is the solid effective pressure and κ_s and η_s are the solid effective dilatational and shear viscosities, respectively. Various constitutive expressions for these quantities have been advanced in the literature, most of them being based on a generalization of the mathematical theory of granular gases. The model developed by Gidaspow (1994), in particular, is often employed. Nevertheless, in the light of the above considerations, we conclude that, for the systems of interest in this work, a Newtonian closure is unsuitable. Furthermore, as said before, since Gidaspow’s model neglects the contribution of shear-induced particle velocity fluctuations to the effective stress of the solid, this is grossly underestimated.

Based on what we have discussed, if we subscribe to the view that, in suspensions far from the jamming state, the mixture effective stress is dominated by the kinetic energy dissipation process occurring in the interstitial fluid, then it is reasonable to assume that the effective stress tensors of the fluid and of the mixture should be nearly identical. We thus believe that the correct closure for $⟨S⟩e$ should be that given in Eqs. (3.7) and (3.8), for dilute suspensions, and Eq. (3.9), for semi-dense suspensions. For the solid effective stress tensor, no empirical closures are available, and, among the few theoretical ones, that of Buyevich (1996) seems to be the only one of general validity (that of Morris and Boulay, 1999, conversely, is valid solely for simple viscometric flows). This is, accordingly, the only viable closure for complex, non-viscometric flows; we stress, nevertheless, that more validation work is necessary to test its accuracy. In the literature, various authors have used Eqs. (3.11) and (3.12) to close $⟨S⟩e$ and Eq. (3.16) to close $⟨S⟩s$ (e.g., Chen et al., 2009; Gidaspow and Huang, 2009; Kaushal et al., 2012; Wang et al., 2013; and Ofei and Ismail, 2016). We believe that doing so is incorrect because (a) the dependence of the mixture viscosity on solid volume fraction is erroneously modeled, η_m decreasing with ε_s, (b) the non-Newtonian rheology observed in semi-dense suspensions is not captured, and (c) the effective solid stress closure does not account for the effect of shear-induced particle velocity fluctuations and accordingly underestimates shear-induced particle migration.

The constitutive expressions presented above overcome the closure problem for the effective stress tensors of the fluid and solid phases in multifluid models. These equations, however, solve also the closure problem of the mixture effective stress, insofar as this is given by the sum of the fluid, solid, and diffusion stress tensors. The first two are now closed, whereas the third does not require closure. Naturally, in this implementation of the mixture model, all the issues affecting the stress closures in the multifluid model affect as well the mixture model. The two models, consequently, are equivalent (provided, of course, that the conditions on which the latter is based are fulfilled—the Stokes number, in particular, has to be much less than unity) and accordingly yield the same predictions and share the same limitations.

In this section, we use a commercial computational fluid dynamics code to solve the multifluid and mixture model equations, with some of the constitutive equations just discussed, to describe the behavior of mixtures of particles (assumed to be monodisperse) for which both modeling approaches apply. To validate the results of the simulations, we use experimental data from the literature. To this end, we consider two model systems, well characterized experimentally, for which both rheological data and flow characteristics are provided. The first is the laminar flow of a dilute suspension in a horizontal pipe; here the validation is conducted in terms of pressure drop. The second is the flow of a dense suspension through an abrupt axisymmetric expansion; here the validation is performed in terms of recirculation length downstream of the expansion. In the first flow, the focus is on the closures for the fluid and mixture effective stress tensors because these significantly influence pressure drops through pipes; the second flow, conversely, focuses on the closure for the solid effective stress tensor, insofar as this affects the shear-induced particle migration and in turn the solid volume fraction radial profile at the expansion inlet, which strongly influences the resulting recirculation lengths.

As model system, we select the flow of phosphate slurries in horizontal pipes investigated experimentally by Fam et al. (1987). The suspension is assumed to consist of monodisperse particles with radius a = 5 μm and density ρ_s = 2650 kg/m³ dispersed in water (density ρ_e = 997 kg/m³ and viscosity μ_e = 0.001 Pa s). The flow is laminar. The pipe diameter is D = 2 mm, and in the experimental apparatus, the pipe was long enough to let the velocity profile fully develop. Because simulating the development of the velocity profile along the entire pipe is computationally too demanding (timewise), we simulate the flow in a pipe of length 10D, at the inlet of which we assume the velocity profile to be fully developed. The sedimentation velocity of the particles is so small that along the real pipe the vertical motion of the particles induced by gravity is negligible. Also, in the real pipe, shear-induced particle migration does not occur, insofar as the time scale of this phenomenon is far larger than the mean residence time of the mixture in the real pipe. As the suspensions considered in this work are dilute (solid volume fraction up to ca. 7%), this is expected, shear-induced particle migration usually being significant only at quite higher concentrations of solid. Consequently, the particle uniform distribution present at the inlet of the pipe is retained throughout the entire real pipe.

The rheological characterization of the mixture shows that, despite being dilute, this has non-Newtonian behavior with yield stress. As stated in the work of Fam et al. (1987), this is due to the irregular shape of the particles and to their polydispersity. In the simulations, we assume that the particles are spherical and all identical, taking into account their real properties solely through the mixture rheology. Fam et al. fitted the model of Herschel and Bulkley (1926) to the mixture viscosity, reporting the corresponding rheological parameters (consistency, power-law index and yield stress) as functions of the solid concentration. We can thus use this experimental data to model the rheology of the suspension.

Our objective is predicting the pressure drop through the pipe, or equivalently the Fanning friction factor, adopting both the multifluid and mixture modeling approaches and different closures for the fluid, solid, and mixture effective stress tensors and validating the results with the experimental data provided by Fam et al. (1987). The cases considered in the simulations are reported in Table I. The CFD code we use is Fluent 17.2, in which we implement the required closures by means of User Defined Functions (UDFs); these include the expressions of Batchelor and Green (1972) and Herschel and Bulkley (1926) for the suspension viscosity, that of Buyevich (1996) for the solid effective stress tensor, and the closure of Mazzei and Lettieri (2007) for the mean fluid-particle drag force per unit volume of suspension. As Table I indicates, in some simulations, we also use the model of Gidaspow (1994), based on an extension of granular kinetic theory, to express the solid effective stress tensor; this is built-in in the software, and so no UDF is required for it. The parameters and constitutive equations adopted in this model are listed in Table II. We use a three-dimensional computational domain and, after a preliminary mesh-independence study, we select a structured grid with 244 260 hexahedral elements for the simulations. A fully developed parabolic velocity profile and a uniform solid volume fraction profile are considered as boundary conditions at the inlet of the pipe. The walls are assumed no-slip, and the pressure at the exit of the pipe is assumed to be atmospheric.

The results of the simulations, for a solid volume fraction of 4.23%, are compared with the experimental data of Fam et al. (1987) in Fig. 1. This reports the Fanning friction factor, which is related to the pressure drop through the pipe, as a function of a generalized Reynolds number. The definitions of these parameters, particularly for the Reynolds number, are quite involved; so, for briefness, we do not report them. They can be found in the work of Fam et al. (1987).

First of all, we notice that cases E1 and M1, which are based on the same closures, yield the same results. This is expected, when the assumptions whereon the mixture model is derived are fulfilled. Cases E4 and M3 also give nearly equal results; we ascribe the small difference between their values to the different numerical implementation of the two models in the code. Cases E1 and M1 yield exactly the same results because these adopt default closures—and so one does not have to modify the code to implement the constitutive equations. The predictions of the multifluid and mixture models coincide, when default closures are used.

Cases E1, E2, and E3 also yield equal results. These simulations differ in the closure adopted for the solid effective stress tensor. However, the pressure drop is significantly affected solely by the fluid effective stress because the mixture stress is dominated by the kinetic energy dissipation process occurring in the interstitial fluid. The effective stress of the solid is mainly responsible for shear-induced particle migration, which in the flow being investigated does not take place (this is confirmed by the numerical results, which indicate that the solid concentration is uniform). Hence, in this flow problem, the solid effective stress plays a negligible role. This is why these three cases (E1, E2, and E3) yield identical predictions. In the light of this, in the other cases considered, we set the solid effective stress tensor to zero.

When the fluid effective stress tensor is closed with Eqs. (3.11) and (3.12), the pressure drop through the pipe is underestimated. This is because Eq. (3.12) underestimates the viscosity of the mixture, which is modeled as a decreasing function of the solid volume fraction. Albeit the suspension is dilute, for most points, the error is significant. In turbulent flow simulations, this issue would not reveal itself because turbulent viscosity would dominate the laminar one, but this does not justify the use of an incorrect closure.

Case M2 is based on the constitutive equation of Batchelor and Green (1972), where the mixture viscosity increase with solid concentration is captured. The predicted pressure drop is indeed higher than that obtained in the cases considered before. However, this theoretical expression does not correctly describe the rheology of the mixture, and consequently the numerical results are still incorrect. Experimental evidence reveals that the suspension is non-Newtonian and behaves like a shear-thinning fluid with yield stress. So, to recover the correct pressure drop values, one must employ an empirical constitutive equation that correctly describes this behavior. This is performed in cases E4 and M3, for which the results closely agree with the experimental data. In case E4, we use the equation of Herschel and Bulkley (1926), with the experimental values of Fam et al. (1987) for the parameters featuring in it, to close the fluid effective stress tensor, whilst in case M3 the same expression is used for the mixture effective stress tensor.

All the findings reported above refer to a single solid volume fraction (this is radially uniform and equal to 4.23% at the pipe inlet and remains uniform throughout the pipe). To investigate the effect of this variable, we run additional simulations for suspensions with a solid volume fraction of 6.73%, for which experimental data are also available. The cases considered in the simulations are those shown in Table I, with the exception of E2 and E3, because these give equal predictions to E1. Table III reports the results of the simulations, which refer to a mean inlet mixture velocity of 0.33 m/s; the Reynolds number values for ε_s = 4.23% and 6.73% are 158 and 68, respectively.

The experiments reveal that the pressure drop increases significantly with solid concentration. This trend is not captured in cases E1 and M1 because in these the fluid effective viscosity (which is essentially equal to the mixture effective viscosity) is modeled with Eq. (3.12), where it decreases with the solid volume fraction. Case M2, in which the fluid effective viscosity is closed with the equation of Batchelor and Green (1972), does capture the trend, but the results considerably underestimate the pressure drop, insofar as the closure does not reflect correctly the viscosity rise with solid volume fraction. Only cases M3 and E4 accurately capture trend and pressure drop values, for in these the constitutive equation for the fluid (and mixture) effective viscosity properly describes the rheology of the suspension, being based on experimental data. The deviation between numerical and experimental results is, in such cases, within 5%. Note also that, for both solid concentrations, the simulation results are equal for cases E1 and M1 and nearly so for cases M3 and E4. The reason for this was discussed before and has to do with numerical implementation.

This analysis just conducted reveals that in laminar flows of dilute suspensions with complex rheological behavior (caused, for instance, by the irregular shape of the particles or the polydispersity of the particle size) it is essential that the effective viscosity of the mixture be correctly modeled. This usually requires the use of empirical closures, which, in the multifluid modeling approach, can also be employed to model the effective viscosity of the fluid. In these systems, the solid effective stress has no prominent role. In dense suspensions, conversely, this term, responsible for shear-induced particle migration, may not be dominant but is no longer negligible, as we shall see in Sec. IV B.

As a more sophisticated system, we simulate the flow of a dense suspension through an abrupt expansion. An abrupt change in the flow cross section leads to high gradients in velocity and recirculation patterns; this makes such a flow a popular benchmark problem for comparing theory and experiment, in both laminar and turbulent regimes (refer, for instance, to Mollicone et al., 2017; 2018). To validate the model predictions, we compare our results with the experimental data available in the work of Moraczewski et al. (2005). They provide rheometry data and values of the recirculation length after the expansion for two types of mixtures of neutrally buoyant particles. The first consists of polymethyl methacrylate (PMMA) particles with average radius a = 42.5 μm and density ρ_s = 1180 kg/m³ dispersed in a liquid with viscosity μ_e = 0.023 Pa s, whilst the second consists of polystyrene particles with a = 128 μm and ρ_s = 1045 kg/m³ suspended in two different liquids, one with μ_e = 0.03 Pa s and the other three times as viscous. All of the suspensions are slightly shear-thinning, their viscosity being well described by Eq. (3.9). The tube diameters upstream and downstream of the expansion are D₁ = 0.476 cm and D₂ = 1.910 cm, respectively, resulting in a 1:4 expansion. The lengths of the upstream and downstream tubes are, respectively, 128D₁ and 80D₂. In the upstream tube, the characteristic velocity of the first mixture is in the range of 3-26 cm/s, while that of the second mixture is in the range of 11-22 cm/s. The characteristic velocities in the downstream tube are 1/16 of the values just reported. The Stokes number is in the range of 10⁻⁵–10⁻⁴ for the first mixture and about 10⁻⁴ for the second mixture.

Simulating the entire experimental setup is very demanding computationally. Note that, even if the flow investigated is steady, the numerical implementation in the CFD code of some closures (in particular, that of Buyevich) requires us to run a number of simulations in transient conditions; this was also true for the system considered in Sec. IV A. This increases even more the computational time. So, in the simulations, we reduce the length of both upstream and downstream tubes, setting them to 5D₁ and 10D₂, respectively. Nevertheless, doing so poses a problem regarding the boundary conditions because, whilst the axial velocity and solid volume fraction profiles are known at the inlet of the actual tube (the profiles are uniform), they are unknown at the location where our computational domain begins (that is, at the inlet section of the shortened upstream tube). These boundary conditions should be obtained by solving the flow along the entire upstream tube; having decided not to do this, we must find a different way forward. Since the Reynolds number in the upstream tube has unit order of magnitude, the real upstream tube is sufficiently long for the velocity profile to fully develop. Therefore, with reasonable accuracy, in our simulations, we can assume a parabolic velocity profile at the inlet of the domain. Doing the same for the solid volume fraction boundary condition, however, is not necessarily correct, insofar as the development length of the solid concentration profile is much longer than that for the mixture velocity. This is due to the large time scale of the shear-induced particle migration process in the inlet tube, which strongly depends on the mean solid concentration as well as the particle/tube radius ratio. Moraczewski et al. (2005) ran experiments over broad ranges for these parameters; in some, the solid concentration profile at the end of the upstream tube was fully developed, whereas in others it was still developing. To overcome this impasse, in the simulations, we consider two limiting cases. In one, we assume that no particle migration occurs so that the solid concentration profile at the inlet of the domain is uniform; in the other, we assume that shear-induced particle migration has enough time to let the concentration profile fully develop. The fully developed profile may be determined analytically, from the steady-state solution of the linear momentum balance equation for the suspension flow in a circular tube. The solution is reported in the work of Buyevich and Kapbsov (1999). This solves the boundary condition issue. By considering these limiting cases, we can investigate how shear-induced particle migration in the upstream tube affects the flow characteristics in the expansion region. As stated in the work of Moraczewski et al. (2005), this phenomenon dictates the major features of the flow after the expansion. For experimental conditions in which the concentration profile is only partly developed, the experimental results (in terms of recirculation length after the expansion) are expected to lie between the numerical predictions for the limiting cases just mentioned.

As shown in Sec. IV A, the multifluid and mixture modeling approaches are equivalent, provided they employ the same closures. Thus, here we use only the latter. To close the mixture effective stress tensor, we use Eq. (3.9) because, as revealed by experiments, the mixture behaves as a generalized Newtonian fluid. The values of the parameters in the closure are obtained by fitting Eq. (3.9) to the experimental data of Moraczewski et al. (2005). The model also requires a constitutive equation for the solid effective stress tensor. While in the upstream tube the closures of Buyevich (1996) and Morris and Boulay (1999) are both applicable (since here the flow is simple shear), in the region after the expansion, only the former can be used (since here the flow is no longer viscometric). So, in the simulations, we must use Buyevich’s relation. Conversely, to determine (analytically) the inlet boundary condition for the solid volume fraction, one may use either closure. The steady-state solid concentration profiles they yield are very similar, and the slight difference between them does not affect the simulation results significantly. Thus, to be consistent, here we report only the results of the simulations based on the boundary conditions obtained with the equation suggested by Buyevich (1996).

In the recirculation region after the expansion, the flow pattern is dictated by convection, because, as the scaling analysis of Moraczewski et al. (2005) indicates, the convective time scale is far shorter than the time scale of shear-induced particle migration. To verify this, we simulate each case study twice, once accounting for the solid effective stress tensor and once neglecting it (by setting it to zero). We remind that, as Eq. (2.21) reveals, for suspensions of neutrally buoyant particles this term is the only cause of particle migration. Thus, if in the expansion region shear-induced particle migration is really negligible, the results of each simulation pair should be identical.

In light of the considerations above, the cases considered in the simulations, in terms of combinations of inlet boundary conditions and constitutive equations, are reported in Table IV. These are used to simulate the flow of the suspensions described before for different solid volume fractions (from 40% to 50%) and Reynolds numbers, with the latter defined as ρ_e(D₂/2)U/⟨η_m⟩, where U represents the mean velocity of the mixture in the downstream tube and ⟨η_m⟩ represents the value of the mixture viscosity (at a given solid concentration) averaged over the shear rate variable. Since the mixtures are only mildly shear-thinning, these averages are quite close to the value of the consistency [denoted as K in Eq. (3.9)].

For the shortened geometry, after a grid study, we select a hexahedral grid with 258 880 cells. The walls of the tubes are assumed no-slip, and the pressure at the exit section is assumed to be atmospheric. The aim of the simulations is reproducing the recirculation length x_R after the expansion; in the simulations, we take this to be the axial distance between the expansion and the further edge of the last computational cell in which the axial mixture velocity component is negative. Figure 2 shows a schematic of the flow system, along with the velocity field (magnitude and streamlines) for one of the cases simulated.

Figure 3 reports the normalized recirculation length x_R/h, where h represents the difference between the large and small tube radii, as a function of the flow Reynolds number for a suspension of particles with a mean radius and volume fraction of 42.5 μm and 40%, respectively. The sets of simulation results presented refer to the cases given in Table IV. All the simulations are able to capture the increase in recirculation length with the Reynolds number. As expected, cases M1 and M2, as well as M3 and M4, give equal results, confirming that in the expansion region shear-induced particle migration has no effect on the flow pattern. The difference between the predictions of cases (M1, M2) and (M3, M4), conversely, clearly highlights the strong influence that shear-induced particle migration occurring in the upstream tube does have on the flow pattern developed after the expansion. The experimental results lie between those of these two pairs of cases, indicating that the flow investigated does not fully develop in the upstream tube. This is in agreement with the estimation made by Moraczewski et al. (2005) based on the experimental formula of Hampton et al. (1997).

As revealed theoretically by Nott and Brady (1994) and observed experimentally by Hampton et al. (1997), the length required by the solid concentration profile to fully develop decreases considerably when a mixture becomes denser. Thus, one expects that at higher concentrations the assumption of having a fully developed solid concentration profile at the inlet of the computational domain should be more likely met. To investigate this, we simulate the motion of the same mixture (particles with a mean radius of 42.5 μm) for two additional mean solid volume fractions, namely, 48% and 50%. Figure 4 presents the results, including those obtained for the volume fraction of 40%. The values of the Reynolds number corresponding to the three volume fractions, from the lowest to the highest, are 0.70, 0.02, and 0.02. At higher solid concentration, the simulations predict the recirculation length quite accurately, when the fully developed concentration profile is considered at the computational domain inlet—this indicates that in the related experiments this profile is fully developed (or nearly so). The results confirm that in the expansion region shear-induced particle migration is negligible; so, the solid effective stress tensor can be set to zero without altering the results.

In the cases considered in Fig. 4, we observed that for ε_s = 0.40 the predictions of cases M3 and M4 do not match the experimental result, and we ascribed this mismatch to the fact that in the related experiment the solid concentration profile is not fully developed. The length required by the solid concentration profile to fully develop is also a strong decreasing function of the particle/tube radius ratio; so, the mismatch should be smaller for a suspension, flowing in the same tube, of equal solid concentration but bigger particles. To investigate this, we simulate the flow of the second suspension described at the beginning of this section, which involves particles with a mean radius of 128 μm, at volume fractions of solid equal to 40%, 45%, and 50%. The corresponding values of the Reynolds number are 0.1, 0.04, and 0.004, respectively. The simulations are based on case M3: the inlet boundary condition adopts the fully developed solid concentration profile and the solid effective stress tensor is set to zero. This is justified in the light of the previous considerations: case M4 would give equal results to case M3, while cases M1 and M2 would give poor predictions, which would become poorer as the suspension becomes denser. The results are shown in Fig. 5. The model predicts the recirculation length quite accurately for all the cases considered—including that for ε_s = 0.40—confirming the validity of our expectations.

In this article, after having briefly described the multifluid and mixture modeling approaches, we focused on the problem of closure of the effective stress tensors of the fluid, solid, and mixture phases. We first presented the (unclosed) expressions of these tensors generated by the (soft volume) averaging scheme, considering the contributors featuring in them and discussing their origin and physical meaning. Subsequently, we reviewed the mechanisms yielding stress, examining their relevance for liquid-particle suspensions and debating their allocation between the phasic stress tensors. With this insight, we discussed a number of closures commonly adopted for these tensors, highlighting the shortcomings that some of these present in reproducing key traits of these systems, namely, the rise in mixture viscosity and the change in mixture rheology from Newtonian to non-Newtonian with solid concentration and the particle migration induced by nonuniform shear. Finally, in the light of the considerations made, we made some recommendations about how the closure problem should be overcome for suspensions far from the jamming state. In multifluid models, the effective stress tensor of the fluid can be taken to be equal to that of the mixture [in the sense clarified in Sec. III B when referring to Eq. (3.9)]. To describe the rheology of the latter accurately, one often has to use empirical closures, especially for suspensions that are concentrated, polydisperse, or made of highly non-ideal particles (e.g., non-spherical or cohesive). The solid effective stress tensor must also be modeled accurately, for it has a strong bearing on particle migration. Its closure poses a considerable problem because it is hard to obtain both experimentally and theoretically; the only constitutive equation of general validity appears to be that derived, through simple heuristic arguments, by Buyevich (1996). This has been successful at capturing the process of shear-induced particle migration in a number of simple stationary viscometric flows but still needs to be thoroughly tested in transient conditions. The equations just discussed solve the closure problem in the mixture model as well because this is merely an approximation of the multifluid model.

In the second part of the paper, to support the ideas presented in the first part, we simulated the behavior of dilute and dense mixtures using the multifluid and mixture modeling approaches and some of the closures formerly reviewed, validating the results against experimental data available in the literature. We showed the equivalence of the two modeling approaches, when these use the same closures (and, of course, provided the assumptions on which the mixture model is based are satisfied). We confirmed the inadequacy of treating the fluid and solid phases as Newtonian continua, regarding the fluid viscosity as proportional to the fluid volume fraction and closing the solid effective stress tensor with expressions based on a generalization—intended for gas-fluidized suspensions—of granular kinetic theory. Numerical results and experimental data agreed only when the complex rheology of the suspensions considered was accurately modeled empirically and the solid effective stress tensor was closed with relations that correctly captured the shear-induced particle migration process. In flows where this has sufficient time to act, its influence on the flow pattern can be considerable, as the second case study investigated clearly revealed.

The authors would like to acknowledge the EPSRC (Grant No. EP/N024915/1) for the financial support offered for this work. They are also grateful to Dr. Jurriaan J. J. Gillissen for the stimulating discussions.

In this appendix, using a perturbation method, we derive Eq. (2.21) of the main article; this equation yields the slip velocity between the fluid and solid phases in the mixture model.

The equations required in the derivation are those governing the evolution of the linear momentum fields of the fluid and solid phases, that is, Eqs. (2.2) and (2.10) of the main article. To lighten the notation, we express these equations as follows:

The angular brackets representing volume averaging have been omitted, and the densities have been regarded as constants. These equations are unclosed. In the passages below, to be concrete, we adopt specific closures for the phasic effective stress tensors and the fluid-particle interaction force, but the final result, Eq. (2.21), is of general validity. This is because what counts in the mathematical derivation are the relative magnitudes of the terms present in the equations, which we assume to be independent of the closures used.

To express the effective stress tensor of the fluid, we assume that the medium is Newtonian; accordingly, the constitutive equation is

where p_e is the fluid pressure, I is the identity tensor, τ_e is the deviatoric part of the fluid stress tensor, η_e is the effective viscosity of the fluid, and $γ̇e$ is twice the fluid rate of deformation tensor. This closure does not account for the bulk viscosity contribution, which is usually negligible.

To express the effective stress tensor of the solid, we also assume that the medium is Newtonian; hence, the constitutive equation is

where p_s is the solid pressure, τ_s is the deviatoric part of the solid stress tensor, η_s is the effective viscosity of the solid, and $γ̇s$ is twice the solid rate of deformation tensor; also in this instance, we have neglected the bulk viscosity contribution, its role being usually minor.

To model the fluid-particle interaction force, we employ the constitutive equation below, which accounts for the buoyancy, drag, virtual mass, and lift forces,

The first term on the right-hand side is the buoyancy force. Various definitions can be adopted for it (refer, for instance, to Jackson, 2000); the most popular considers only the isotropic part of S_e, the force being equal to −ε_s∂_xp_e, but we think that the definition used above should be favored.

The second term on the right-hand side of Eq. (A5) is the drag force, β denoting the drag force coefficient. Several expressions for it are available in the literature; for briefness, we do not report any of them here. The interested reader can refer, for instance, to Jackson (2000) and Mazzei and Lettieri (2007).

The final term on the right-hand side of Eq. (A5) is the sum of the virtual mass and lift forces; ζ represents their coefficients, which are assumed to be equal, so that the resultant force satisfies the principle of frame indifference (Drew and Passman, 1998). Often, ζ is assumed to be equal to 1/2. This value holds solely at low solid concentration; for moderately denser systems, Zuber (1964) suggested that ζ = (1 + 3ε_s)/2.

The substantial derivatives featuring in the virtual mass force expression represent the local acceleration of the two phases and are defined as follows:

Using the continuity equations, we may express the linear momentum balance equations in Lagrangian form; if Eq. (A5) is adopted as closure, we obtain

and

As we can see, this equation contains a term related to the fluid effective stress tensor; this enters the equation through the constitutive equation used for the buoyancy force.

Our aim is deriving an algebraic equation that relates the fluid-particle slip velocity to variables which are either known or that may be calculated using the equations of the mixture model. To this end, let us subtract Eq. (A8) from Eq. (A7). This gives

where

To derive the algebraic equation mentioned above, we use a perturbation method (Lin and Segel, 1988 and Bush, 1992). The first step is scaling Eq. (A9). To do so, we introduce the variables

These have been scaled by dividing the original variable by an appropriate characteristic quantity. The scales of the dependent variables are selected so as to render the order of magnitude of the scaled variables equal to unity, whilst the scales of the independent variables are selected so as to render the order of magnitude of the derivatives of the scaled dependent variables equal to unity.

The volume fractions of fluid and solid, along with α, β, γ, and ζ, are variables and thus should be scaled; however, to simplify the analysis, we treat them as parameters. This does not invalidate the results. We do the same also for the effective viscosities of the two phases. Using the dimensionless variables introduced above, we now write Eq. (A9) in dimensionless form as

where