Flow of a Dense Suspension Modeled as a Modified Second Grade Fluid

In this paper, a simple shear flow of a dense suspension is studied. We propose a new constitutive relationship based on the second grade fluid model for the suspension, capable of exhibiting non-linear effects, where the normal stress coefficients are assumed to depend on the volume fraction of the particles and the shear viscosity depends on the shear rate and the volume fraction. After non-dimensionalizing the equations, we perform a parametric study looking at the effects of the normal stress coefficients and the variable viscosity. The numerical results show that for a certain range of parameters, the particles tend to form a region of high and uniform volume fraction, near the lower half of the flow.


Introduction
Flows of (solid) particles suspended in a fluid present many interesting problems for mathematicians, engineers and physicist.Natural examples of such flows are mudslides, debris flow (Davies, 1986) [1], snow or ice avalanches.Industrial examples of suspensions are slurries, (wet) granular materials, fluidized beds, etc. (see Govier and Aziz (1972) [2], Zandi (1971) [3], Shook and Roco (1991) [4]).In general, these suspensions show non-linear behavior such as normal stress differences, yield stress, discontinuous changes in stresses, etc. [5].tomathematically model these flows, from a continuum mechanics perspective, there are at least two distinct methods; we can look at the suspension as a single continuum where a single constitutive relation for the stress tensor is proposed and the material properties are variables (function of time and space, and perhaps function of shear rate, volume fraction of the solids, etc.), or alternatively, we can look at the two components as two distinct but interacting continua, using the methods of mixture theory, where two equations are needed for the stress tensors and one equation for the interaction forces.For details of the second approach, namely the mixture theory, we refer the readers to papers by Johnson et al. (1991) [6], Massoudi (2003Massoudi ( , 2008Massoudi ( , 2010) ) [7][8][9].In this paper, we will use the (single component) suspension approach.
From a historical point of view, Reynolds (1885) seems to have been the first scientist who noticed that if a shearing motion in a bed of closely packed particles is to occur, the bed must expand to increase the volume of its voids [10].He named this phenomenon "dilatancy," which is a manifestation of the normal-stress effects.Reiner (1945Reiner ( , 1948Reiner ( , 1958) ) was one of the first who derived a non-Newtonian fluid model to predict dilatancy in wet sand [11][12][13], even though his model did not show how the volume fraction affects the stresses.Massoudi (2011) [14] modified and generalized this constitutive relation, also known as the Reiner-Rivlin [15] fluid model, by suggesting that the shear viscosity would depend on the shear rate and the volume fraction.Later, Wu et al. (2012) considered the shearing motion of such a material [16], and Wu et al. (2013) [17] studied the flow of granular materials down an inclined plane.Their numerical simulations predicted only one non-zero normal stress difference, as is expected in the Reiner-Rivlin fluid type models.
In this short paper, we continue this line of thinking and propose a new constitutive model for the flow of a dense suspension by generalizing the second grade fluid model and assuming that the viscosity depends on the shear rate and the volume fraction, while the normal stress coefficients are only functions of the volume fraction.In the Section 2, the governing equations are introduced.In Section 3, the constitutive equation for the stress tensor is discussed.In Section 4, we show the geometry and the kinematics of the problem and in Section 5, we solve the simple shear flow and perform a parametric study on several different dimensionless numbers.

Governing Equations
We assume that the motion and the behaviour of the suspension can be described using the traditional methods of continuum mechanics.In the absence of any thermo-chemical and electro-magnetic effects, the basic governing equations are the conservation laws for mass, linear momentum and angular momentum (see Slattery (1999) [18]).As we are only considering a purely mechanical system, the energy equation and the entropy inequality are not considered.

Conservation of Mass
The conservation of mass is: where ∂/∂t is the derivative with respect to time, div is the divergence operator, v is the velocity vector, and ρ is the density of the suspension given by: where ρ s is the density of an individual particle (assumed to be constant) and φ(x, t) is the volume fraction of the particles, where 0 ≤ φ < φ max < 1; φ is an independent kinematical field and is a continuous function of position and time.We will see in Section 4, that the introduction of this kinematical field will present new challenges to us, especially in the specification of the boundary conditions.Furthermore, this term couples the density field and the velocity field and allows for finding particle distribution.

Conservation of Linear Momentum
Let T represent the Cauchy stress tensor, then the balance of linear momentum is: where dv dt = ∂v ∂t + (gradv)v and b stands for the body force.

Conservation of Angular Momentum
In absence of couple stresses the Cauchy stress tensor is symmetric, that is, In the next section, we will discuss the constitutive modeling of the stress tensor.

Constitutive Equation for the Stress Tensor
In this paper, we model the dense suspension of particles as a single component non-linear fluid-like material.Rajagopal and Massoudi (1990) [19] and Rajagopal et al. (1994) [20] proposed a model, where the Cauchy stress tensor depends on the symmetric part of velocity gradient and another second order tensor related to the density gradient (see also Massoudi and Mehrabadi (2001) [21]).And more recently, Massoudi and Tran (2016) [22] proposed a modified third grade fluid model for a dense suspension of particles.
In this paper, we focus our attention on suspension models where the normal stress effects are present (see Massoudi (2004Massoudi ( , 2010) ) [23,24]).One of the simplest constitutive equations capable of describing both of the normal stress effects is the second grade fluid, or the Rivlin-Ericksen fluid of grade two (Rivlin and Ericksen (1955) [25], Truesdell and Noll (1992) [26]).For a second grade fluid, the Cauchy stress tensor is given by: where I is the identity tensor, p is the constraint due to incompressibility (pressure), and the kinematical tensors A 1 and A 2 are defined through where L is the velocity gradient, µ is the coefficient of viscosity, α 1 and α 2 are material moduli which are commonly referred to as the normal stress coefficients.We modify the traditional second grade fluid model and assume that the suspension can be modeled as where .γ = 1/2tr(A 1 2 ) is the shear rate, (p(φ) + β r )I is the spherical (isotropic part of the) stress tensor (which includes the pressure effects) (see Serrin (1959) [27]);In this paper, we assume p(φ) = β p where β r can be a function of the principal invariants of A 1 , A 2 and φ, and we also assume that the shear viscosity, µ, depends on the volume fraction and the shear rate and the normal stress coefficients, α 1 and α 2 , depend on the volume fraction.That is, we assume: These expressions can be viewed as the Taylor series approximation for the material parameters (Rajagopal, et al. (1994) [20]).Clearly, in Equation ( 7), the term 'p' does not have the same meaning or implications as the pressure term in Equation (5).More accurately, as pointed out by Rajagopal (2015) [28], 'p' should be referred to as the "mean value of the stress." The thermodynamics and stability of fluids of second grade have been studied in detail by Dunn and Fosdick (1974) [29].They showed that if the fluid is to be thermodynamically consistent, then µ ≥ 0 (9a) It is known that for many non-linear fluids which are assumed to follow Equation ( 9), the experimental values reported for α 1 and α 2 do not satisfy the restriction of Equations (9b) and (9c).In an important paper, Fosdick and Rajagopal (1979) [30] show that irrespective of whether α 1 + α 2 is positive, the fluid is unsuitable if α 1 is negative.For further details on this and other relevant issues in fluids of differential type, we refer the reader to the review article by Dunn and Rajagopal (1995) [31].In certain applications, where the fluid is known to be shear-thickening (or shear-thinning), then modified (or generalized) forms of the second-grade fluid have been proposed (see Man (1992) [32], Massoudi and Vaidya (2008) [33], Man and Massoudi (2010) [34], Massoudi and Tran (2016) [22]).In this paper, we assume that the viscosity changes with the shear rate and can be modeled as a power-law type relation, where µ r is given by: where µ r0 is a reference viscosity, assumed to be a constant, .γ m represents the effect of the shear dependency.We should note that when m < 0, the fluid is shear-thinning, and if m > 0, the fluid is shear-thickening.The above power-law relation, as is well-known, suffers from the fact that for shear-thinning fluids, the model predicts that the viscosity goes to infinity as the shear rate approaches zero.To overcome this short-coming, we can use an expression of the type µ r = µ r0 1 + tr A 1 2 1/2 m which predicts a constant value for the zero shear rate viscosity.
Based on previous studies (Rajagopal and Massoudi (1990) [19]), we assume see also [19,36].We also assume that β r , a parameter which accounts for the compacting (a measure of the rigidity) of the particles, can be given by, where φ c is the critical value of the volume fraction.In this paper, we assume g(φ) = φ.In a sense, C r is a material parameter related to how much the particles can be compacted.For very rigid particles, C r is large, ensuring that the local volume fraction of the particles cannot be larger than φ c .We assume that h(φ) in Equation ( 12) is a smooth step function, such that where S is a parameter related to the slope of the step function.For example, if S is chosen as 3000, then h(φ c − 0.001)/h(φ c ) < 0.01 and h(φ c − 0.005) ∼ o 10 −13 , therefore β r is negligible when φ is not close to φ c .Thus the equation for the stress tensor used in our paper, in expanded form is: where β p , C r , φ c , µ r0 , m, α 10 and α 20 are material constants.Note that in the above equation T goes to zero when φ → 0 , as is expected.Finally, we need to emphasize that Equation ( 15) being a generalization of the second grade fluid model, contains higher order derivatives of the velocity gradient and as a result, once Equation ( 15) is substituted in the linear momentum equation, the order of the ensuing partial differential equations is one order higher than the Navier-Stokes equations.Therefore, in general, depending on the kinematics of the flow, we may need additional boundary conditions.For a detailed discussion of this issue, we refer the reader to Rajagopal (1995) [37].

Geometry and the Kinematics of the Flow
The geometry and the kinematics of the flow are shown in Figure 1.Whenever a new model is proposed or developed, it is worthwhile to consider some simple flow situations, usually referred to as viscometric flows (see Truesdell (1974) [38]).This is to see if analytical solutions can be obtained and whether we can gain some physical insights into the nature of the model and perhaps even compare the results with ideal one-dimensional experiments.In this paper, we consider a simple shear flow between two horizontal flat plates (see [16]).The non-dimensional form of the conservation of linear momentum is: where  is a reference length, for example, the distance between the two plates, and  is a reference velocity, for example the (shearing) velocity of the upper plate,  .In Equation ( 16) the asterisks have been dropped for simplicity.Among the above dimensionless numbers,  and  are related to the normal stress coefficients,  is related to the compressibility factor (the pressure term),  is related to the compacting effect of the particles. and  are always less than zero, indicating that with compression there is a tendency for densification. is the Froude number,  is related to the viscous effects (similar to the Reynolds number), and  is the dimensionless shear rate parameter.For a simple shear flow, we assume: where  is the unit vector along the  direction, then using Equation ( 18), we have, The non-dimensional form of the conservation of linear momentum is: where we have used the following reference parameters: where H is a reference length, for example, the distance between the two plates, and u 0 is a reference velocity, for example the (shearing) velocity of the upper plate, U 0 .In Equation ( 16) the asterisks have been dropped for simplicity.Among the above dimensionless numbers, M 1 and M 2 are related to the normal stress coefficients, B p is related to the compressibility factor (the pressure term), B r is related to the compacting effect of the particles.B p and B r are always less than zero, indicating that with compression there is a tendency for densification.Fr is the Froude number, R is related to the viscous effects (similar to the Reynolds number), and Γ is the dimensionless shear rate parameter.For a simple shear flow, we assume: Fluids 2018, 3, 55 where e x is the unit vector along the x direction, then using Equation ( 18), we have, trD = divv = 0 (20) where D is the symmetric part of the velocity gradient.In the above equations, prime designates the derivative with respect to y. Substituting ( 18)-( 21) into ( 16), the governing equations are simplified and we obtain a system of non-linear ordinary differential equations.The momentum equation in the x-direction is: And in the y-direction is: From Equation ( 22) we can see that dimensionless number, R, can be canceled for this problem.The normal stress differences for this problem are, From Equations ( 22) and ( 23), it is clear that we need two boundary conditions for U and one boundary condition for φ.We use the no-slip condition at both boundaries for the velocity: It is possible that the particles may slip at the walls (see Kim and Rosato (1994) [39], Zhao and Massoudi (2014) [40]).For volume fraction, φ, the appropriate boundary condition may be given as an average value defined in an integral form: or φ could be prescribed at Y = −1: where Θ is the value of the volume fraction at the boundary.In this paper, the condition of the average volume fraction is used (Equation ( 27)).Note that since we are including the effect of gravity, we do not use the symmetry boundary condition for the volume fraction, to allow for the settling of the particles.

Results and Discussion
The system of the non-linear ordinary differential Equations ( 22) and ( 23) with the boundary conditions ( 26) and ( 27) are solved using the MATLAB solver bvp4c, which is a collocation boundary value problem solver [41].The step size is automatically adjusted by the solver.The default relative tolerance for the maximum residue is 0.001.The boundary conditions for the average volume fraction are numerically satisfied by using the shooting method.Since the model we have proposed is a new model, there are no experimental data available to use for the material properties.Therefore, we will perform a parametric study for a limited range of the dimensionless numbers specified in Table 1.

The Effect of the Compacting of the Particles
First, we consider the isotropic (spherical) part of the stress tensor due to the particle contact.Recall that B r is related to the level of the rigidity of the particles and φ c is related to the maximum packing of the particles.Figure 2 shows the effect of B r on the volume fraction profiles.The influence of this dimensionless number on the velocity profiles, for the range of dimensionless parameters selected in this paper, is negligible and as a result we do not show it.The value of B r is always equal to or less than zero which ensures that with compression, the density increases.We can see that due to the effect of the gravity, the particles tend to accumulate near the bottom plate.Figure 2 shows that as B r changes, near the bottom plate the particles cannot be further compacted once the volume fraction has reached φ c ; as a result, a high concentration layer/region seems to be formed.In this region, the volume fraction is uniform; while from the edge of this region, along the y-direction (see Figure 1 for the definition of the coordinates) the gradient of the volume fraction profile suddenly becomes large.When B r is small, near the lower plate the local volume fraction continues to increase; from Figure 2 we can see that the spherical stress, β r , is important even when the magnitude of B r is small, see the cases for B r = 0 and B r = −0.1.

The Effect of the Compacting of the Particles
First, we consider the isotropic (spherical) part of the stress tensor due to the particle contact.Recall that  is related to the level of the rigidity of the particles and  is related to the maximum packing of the particles.Figure 2 shows the effect of  on the volume fraction profiles.The influence of this dimensionless number on the velocity profiles, for the range of dimensionless parameters selected in this paper, is negligible and as a result we do not show it.The value of  is always equal to or less than zero which ensures that with compression, the density increases.We can see that due to the effect of the gravity, the particles tend to accumulate near the bottom plate.Figure 2 shows that as  changes, near the bottom plate the particles cannot be further compacted once the volume fraction has reached  ; as a result, a high concentration layer/region seems to be formed.In this region, the volume fraction is uniform; while from the edge of this region, along the y-direction (see Figure 1 for the definition of the coordinates) the gradient of the volume fraction profile suddenly becomes large.When  is small, near the lower plate the local volume fraction continues to increase; from Figure 2 we can see that the spherical stress,  , is important even when the magnitude of  is small, see the cases for  = 0 and  = −0.1.
From Figure 3, we see that as  decreases, near the bottom plate the thickness of the high concentration layer/region increases; in this region, the local volume fraction is uniform and close to  .As Figure 3a shows, when  decreases, near the bottom plate the velocity increases and the profile becomes more linear in the region with high volume fraction.This may be due to the fact that when  is small, the volume fraction distribution is more uniform, which leads to a more uniform distribution of the viscosity.From Figure 3, we see that as φ c decreases, near the bottom plate the thickness of the high concentration layer/region increases; in this region, the local volume fraction is uniform and close to φ c .As Figure 3a shows, when φ c decreases, near the bottom plate the velocity increases and the profile becomes more linear in the region with high volume fraction.This may be due to the fact that when φ c Fluids 2018, 3, 55 8 of 12 is small, the volume fraction distribution is more uniform, which leads to a more uniform distribution of the viscosity.

Effects of the Normal Stress Coefficients and the (Variable) Viscosity
In our model, the normal stress differences are related to the dimensionless parameters M 1 and M 2 .Their effect appears as the combined parameter (2M 1 + M 2 ), see Equations ( 24) and ( 25).From Figure 4b, we can see that as (2M 1 + M 2 ) increases, near the upper plate, the volume fraction increases; this indicates that more particles have moved to that region.While near the bottom plate, the high concentration region tends to disappear.The velocity increases as (2M 1 + M 2 ) increases.Recall that B p is a dimensionless parameter related to the pressure.As Figure 5b shows, when the magnitude of B p increases, the concentration profiles become more uniform and linear-like, which indicates that the effects of gravity (Fr), spherical stress (B r ), and the normal stress differences (2M 1 + M 2 ) are less important (see Equation ( 23)).As B p increases, the volume fraction distribution becomes more uniform and the velocity profiles become more linear.

Effects of the Normal Stress Coefficients and the (Variable) Viscosity
In our model, the normal stress differences are related to the dimensionless parameters  and  .Their effect appears as the combined parameter (2 +  ), see Equations ( 24) and (25).From Figure 4b, we can see that as (2 +  ) increases, near the upper plate, the volume fraction increases; this indicates that more particles have moved to that region.While near the bottom plate, the high concentration region tends to disappear.The velocity increases as (2 +  ) increases.Recall that  is a dimensionless parameter related to the pressure.As Figure 5b shows, when the magnitude of  increases, the concentration profiles become more uniform and linear-like, which indicates that the effects of gravity (), spherical stress ( ), and the normal stress differences (2 +  ) are less important (see Equation ( 23)).As  increases, the volume fraction distribution becomes more uniform and the velocity profiles become more linear.Another important quantity of interest is the (dimensionless) wall shear stress, | | which is related to the skin-friction coefficient [42].In this problem, | | has the following expression: Table 2 shows the effect of the normal stresses, (2 +  ) and  , on the wall shear stress.As the effect of the normal stresses increases, the wall shear stress also increases.Another important quantity of interest is the (dimensionless) wall shear stress, |τ w | which is related to the skin-friction coefficient [42].In this problem, |τ w | has the following expression: Table 2 shows the effect of the normal stresses, (2M 1 + M 2 ) and B p , on the wall shear stress.As the effect of the normal stresses increases, the wall shear stress also increases.In our model, m is the parameter related to the effects of variable (shear-dependent) viscosity.Recall that m > 0 and m < 0 imply shear-thickening and shear-thinning behavior, respectively.Figure 6 shows the effect of m on the velocity and the volume fraction distributions.As the magnitude of m increases, the volume fraction distribution becomes more non-uniform and the high concentration region near the lower plate becomes more significant when m is positive.As m increases, the velocity decreases.

Effects of Gravity and the Bulk Volume Fraction
Figure 7 shows the effect of Fr number, which is the ratio of the inertial forces to the gravitational force.As Fr number decreases, that is, lower flow rates, more particles seem to accumulate in the lower half of the flow, and therefore the volume fraction near the bottom plate increases.When Fr number is large, the velocity seems to be more linear.The effect of the bulk volume fraction on the flow is shown in Figure 8.The thickness of the high concentration region increases as N increases.Increasing N seems to lead to a more linear-like velocity distribution and more uniform volume fractions.
In our model,  is the parameter related to the effects of variable (shear-dependent) viscosity.Recall that  > 0 and  < 0 imply shear-thickening and shear-thinning behavior, respectively.Figure 6 shows the effect of  on the velocity and the volume fraction distributions.As the magnitude of  increases, the volume fraction distribution becomes more non-uniform and the high concentration region near the lower plate becomes more significant when  is positive.As  increases, the velocity decreases.

Effects of Gravity and the Bulk Volume Fraction
Figure 7 shows the effect of  number , which is the ratio of the inertial forces to the gravitational force.As  number decreases, that is, lower flow rates, more particles seem to accumulate in the lower half of the flow, and therefore the volume fraction near the bottom plate increases.When  number is large, the velocity seems to be more linear.The effect of the bulk volume fraction on the flow is shown in Figure 8.The thickness of the high concentration region increases as  increases.Increasing  seems to lead to a more linear-like velocity distribution and more uniform volume fractions.

Effects of Gravity and the Bulk Volume Fraction
Figure 7 shows the effect of  number , which is the ratio of the inertial forces to the gravitational force.As  number decreases, that is, lower flow rates, more particles seem to accumulate in the lower half of the flow, and therefore the volume fraction near the bottom plate increases.When  number is large, the velocity seems to be more linear.The effect of the bulk volume fraction on the flow is shown in Figure 8.The thickness of the high concentration region increases as  increases.Increasing  seems to lead to a more linear-like velocity distribution and more uniform volume fractions.

Conclusions
In this paper, we study the simple shear flow of a non-linear fluid (a dense suspension of particles in a liquid) between two horizontal plates.The fluid is modeled as a modified second grade fluid.Important rheological properties such as the shear rate and the volume fraction dependency of the viscosity and the normal stress coefficients are considered.Through numerical simulations, we

Figure 1 .
Figure 1.Schematic of the flow.

Table 1 .
The studied dimensionless numbers.

Table 2 .
The effect of the normal stresses (2M 1 + M 2 ) and B p on the wall shear stress.