Next Article in Journal
Research on Location Selection of Personnel Door and Anemometer Based on FLUENT
Next Article in Special Issue
A New Rheological Model for Phosphate Slurry Flows
Previous Article in Journal
Gas Dynamics of Micro- and Nanofluidic Systems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Shear Flows of Dilatant Fluids with Limited Shear Rates: Analytical Results and Linear Stability Analysis

Dipartimento di Matematica e Informatica “U. Dini”, Università degli studi di Firenze, Viale Morgagni 67/a, 50134 Firenze, Italy
Fluids 2023, 8(1), 25; https://doi.org/10.3390/fluids8010025
Submission received: 15 November 2022 / Revised: 3 January 2023 / Accepted: 7 January 2023 / Published: 9 January 2023
(This article belongs to the Special Issue Advances in Computational Mechanics of Non-Newtonian Fluids)

Abstract

:
In this paper, we study the simple shear flows of a class of dilatant fluids with a limited shear rate. This class of fluids is characterized by shear thickening behavior in which the apparent viscosity tends to infinity as the modulus of the stress approaches a finite threshold. The apparent viscosity function is a logarithmic type with two material parameters. We considered this specific form because it fits very well with the flow curves of some granular suspensions for specific values of the material parameters. Despite the nonlinearity of the constitutive law, it is possible to determine explicit steady-state solutions for a simple shear flow, namely (i) the channel flow; (ii) the flow between coaxial cylinders, and (iii) the flow down an inclined plane. We performed a two-dimensional linear stability analysis to investigate the onset of possible instabilities of the steady basic flow, putting into evidence the dependency of the critical Reynolds number on the material parameters.

1. Introduction

Shear thickening (or dilatant) fluids are part of a special class of non-Newtonian fluids in which the viscosity increases with the shear rate. One of the first investigations on shear-thickening fluids is probably the one by Metzner and Whitlock [1] who studied a suspension of titanium dioxide (TiO 2 ) particles in water. They showed that the relation between the shear stress and the shear rate is nonlinear with the apparent viscosity that grows with the increasing shear rate. For the elevated TiO 2 volume fraction, the stress curve appears to almost diverge as the shear rate increases. The most dramatic scenario of this type of behavior is that of discontinuous shear thickening (DST), where viscosity becomes infinite for a finite value of the shear rate and the material goes from fluid-like to solid-like, see [2,3,4].
In the recent past, many works have investigated the nature of shear-thickening and DST fluids, see [5,6,7,8]. In particular, the crucial effects of frictional contact forces (in the framework of jamming transition) on the shear-thickening behaviors of the fluid have been observed. As observed in [9], the scale of the stress (where the stress–shear rate constitutive law starts to show large growth) can be attributed to the interfacial tension forces, which increase abruptly due to the dramatic constraint of jamming.
In general, dense mixtures of granules and liquids behave as dilatant fluids, a classical example being cornstarch and water. Some fluids, such as high-concentration suspensions, insoluble polymer systems, clay, and quicksand display shear-thickening behaviors [10]. Although a large majority of fluids exhibit pseudo-plastic (or shear-thinning) behaviors [11], shear-thickening fluids are also ubiquitous in a variety of practical applications, from damping devices to machine mounts, as well as in the manufacturing of body armor [12], where the shear-thickening properties are exploited to dampen the impacts of a bullet.
Regarding the counterintuitive behaviors of dilatant fluids, one can search the web for a video of one running on a pool filled with cornstarch and water as if it was an actual rigid surface, but sinking into it when walking or standing still. At the microscopic level, experiments have shown that dense suspensions are basically frictionless for low confining pressure, while the transition from low to high friction occurs with normal stress near the critical stress, see [13].
The constitutive law of a shear-thickening fluid is commonly represented by a function in which the shear stress is expressed as a function of the shear rate and in which the first derivative of that function increases with the shear rate. The most common example is that of power law dilatant fluids, where the apparent viscosity is a given positive power of the shear rate modulus. In the latter model, any arbitrarily large value of the shear rate can be attained if a sufficiently large shear stress is applied (the strain is not limited). This type of model is not appropriate when one considers a dilatant fluid in which the flow behavior is similar to that of a discontinuous shear-thickening fluid. In this case, we may consider a constitutive law in which there is a critical value to which the shear rate tends when the stress becomes infinite (limited strain model), see [4,14].
A very basic limited strain model is the one considered in [4] for a system formed by a mass, a dashpot, and an inextensible string; see Figure 1a. If the string is not completely extended, the response of the system is that of a Newtonian fluid, but when the string becomes fully extended, the response is that of a rigid body (the applied force produces no deformation). The constitutive relation is, thus, the one shown in Figure 1b, which can be seen as the “dual” of the Bingham constitutive law (the fluid flow only if the applied stress is below a certain threshold). The one-dimensional unsteady channel flow of a fluid with the constitutive law, as shown in Figure 1b, was studied in [4].
Let us suppose that the fluid is incompressible so that the Cauchy stress tensor can be written (in the paper, the starred quantities are always dimensional) as T = p I + S , where p is the Lagrange multiplier due to the incompressibility constraint and S is the deviatoric part of the stress. The constitutive equation for a Newtonian fluid with a limited shear rate can be seen in [4,14].
S = 2 μ D , S τ c , D = d c , S τ c ,
where μ is the viscosity ( [ μ ] = Pa·s), τ c is the critical shear stress ( [ τ c ] = Pa), and d c is the critical shear rate ( [ d c ] = s 1 ). The quantities
S = 1 2 S · S 1 / 2 , D = 1 2 D · D 1 / 2 ,
represent the norm of the deviatoric part of the stress S and the symmetric part of the velocity gradient D , respectively. Looking at Figure 1b, it is evident that if we plot S vs. D , we obtain a graph (the stress is not determinate for D = d c ); hence, it is more natural to express D vs. S , so that the relation is an actual function (see [15] for a discussion on the material defined by implicit constitutive equations).
Relation (1) has a singularity in the first derivative at ( d c , τ c ) , which can be smoothed if one takes, for instance,
S = 2 μ eff ( D ) D ,
where the effective viscosity is given by
μ eff ( ξ ) = 1 2 ξ a ln b b 2 ξ .
The constitutive relation (2), (3) was selected because it fits the experimental flow curves of dense suspensions, see [1,16].
In Figure 2, we plotted the relation (2), taking the norms on both sides. We normalized S with 1 / ( 2 a ) and D with b / 2 . The parameter a has the dimension of the inverse of the pressure, while the dimension of b is that of a shear rate. In Figure 3, we plot the effective viscosities (3) for some values of a and with b = 18 s 1 . The ranges for these values are consistent with the experimental data of [16]. We notice that for sufficiently large values of a , the constitutive behavior of the fluid is similar to that of an inviscid fluid (the effective viscosity is nearly zero almost everywhere in the range ( 0 , b ) ) except for a small boundary layer near the critical value b , where the viscosity becomes arbitrarily large (rigid body behavior). The behaviors of the fluids for small values of D are similar to those of Newtonian fluids with constant viscosity
lim ξ 0 μ eff ( ξ ) = 1 a b .
Hence, parameters a and b determine the “almost constant” viscosity of the fluid for small shear rates, whereas b / 2 is the limiting value for the shear rate.
Figure 2. Constitutive law (2) in which the shear rate is normalized with b / 2 and the shear stress with 1 / 2 a .
Figure 2. Constitutive law (2) in which the shear rate is normalized with b / 2 and the shear stress with 1 / 2 a .
Fluids 08 00025 g002
The behaviors of parameters a and b with the particle concentration ϕ are related to the shear-thickening effects due to the particle concentrations. From [1,16], it is evident that the parameter b is a decreasing function of ϕ . This is physically consistent since the limiting shear rate increases if the particle’s volume fraction decreases because the frictional effects between particles become less pronounced. The same is true for a as the apparent viscosity (which is proportional to 1 / a ) must decrease as the volume fraction decreases. The exact dependencies of the parameters a and b on ϕ will surely involve some material constants that take into account the particular suspension considered. The definition of such a dependency is beyond the scope of this paper and certainly would require more experimental data.
The constitutive law (3) does not take into account the possibility that the slope of the constitutive curve may start to decrease after a certain value of the shear stress. In practice, it does not allow one to consider the S-shaped behavior exhibited by the suspension. Here, we are only interested in the range of stresses for which the slope of the shear stress shear rate curve is monotonically increasing and tends to be infinite, such as the ones of [1,16].
Figure 3. Plot of the effective viscosity μ eff given in (3) as a function of the shear rate 2 ξ for different values of a and b = 18 s 1 . The values of these parameters are taken from [16].
Figure 3. Plot of the effective viscosity μ eff given in (3) as a function of the shear rate 2 ξ for different values of a and b = 18 s 1 . The values of these parameters are taken from [16].
Fluids 08 00025 g003
In this work, we consider the responses of types (2) and (3) since we have seen that they describe quite well the experimental flow curves of some dense suspensions. In particular, we consider Figure 4 and Figure 5, in which a and b are selected to fit the experimental curves of two types of suspensions, namely a 47 % TiO 2 suspension in Figure 4 and a 50 % cornstarch suspension in Figure 5. The experimental data in Figure 4 and Figure 5 are taken from [1,16] respectively. This paper is devoted to the study of some simple shear flow of fluid with a constitutive law of types (2) and (3). In particular, we focus on three types of flow, i.e., (i) planar channel flow; (ii) the flow between coaxial cylinders; and (iii) the flow down an inclined plane. Despite the nonlinearity of (3), we shall see that in all geometrical settings considered in this paper, it is always possible to determine an analytical steady-state solution. For each problem, the flow equations are scaled appropriately so that the velocity field can be expressed in terms of the Reynolds number and other non-dimensional material (i.e., not depending on the flow) parameters. Because of the boundedness of the shear rate, in some geometrical settings, the Reynolds number happens to be bounded, meaning that the range of velocities that can be attained by increasing the force that drives the motion is limited. This is consistent with the shear thickening and limited strain nature of the fluid, where friction forces between layers (viscosity) grow with the shear stress, preventing the velocity from increasing unboundedly. Models of types (2) and (3) are useful because one can explicitly determine the velocity and stress profiles in the simple flow, but also because they allow approximating the DST models by avoiding the intrinsic singularities in the constitutive equations.
In addition to determining the steady-state solutions of the basic flow, we performed a modal stability analysis to investigate the onset of possible instabilities. In the flow considered here, the critical Reynolds number is an increasing function of both a and b . Looking at Figure 3, we may explain this behavior by observing that the increase of a (this is also true for the increase of b ) results in the formation of boundary layers with high viscosity near the rigid walls that tend to stabilize the flow. This type of behavior has also been put into evidence in a recent paper [17] for the stress power-law fluids.
The paper is organized as follows. In Section 2, we derive the mathematical formulation of the problem and focus on (i) planar channel flow; (ii) the flow between coaxial cylinders; and (iii) the flow down an inclined plane. In Section 3, we perform the linear stability analysis for the above-mentioned simple flow. The last section is devoted to conclusions and perspectives.

2. The Mathematical Model

Let us consider the following constitutive equation
T = p I + S , S = 2 μ eff ( D ) D ,
μ eff ( D ) = 1 2 D a ln b b 2 D ,
and let us rescale the dimensional variables with
x = L x , t = ( L U 1 ) t ,
p = Π p , S = a 1 S , D = ( U L 1 ) D .
Introducing the quantities b = L b U 1 and Π = Π a , we write the non-dimensional mass and momentum balance (neglecting body forces) as
div v = 0 , Re v ˙ = Π p + div ( S ) ,
where Re = ρ U 2 a is the Reynolds number. The non-dimensional constitutive equation becomes:
S = 1 2 D ln b b 2 D = μ eff 2 D .
We notice that b and Re are related by
b = Γ Re , Γ = L b ρ a
where Γ is a non-dimensional material parameter (it does not depend on the flow). From [1,16], we know that typically Γ = O ( 10 ) .

2.1. Channel Flow Driven by a Pressure Gradient

Let us consider problems (8) and (9) for a fluid flowing in a planar channel of length L and constant width 2 H driven by a pressure gradient Δ p = p in p out > 0 . For simplicity, here we assume that H = L , but minor changes allow one to consider the cases in which they are different. In this situation, it is natural to select Π = Δ p and Π = Δ p a . We look for a stationary solution of problem (8) in the form v = u ( y ) i , p = p ( x , y ) , so that the incompressibility constraint div v = 0 is automatically satisfied. The problem becomes
0 = Π p x + S 12 y , p y = 0 .
We immediately realize that
p x = 1 , D = u y 2 , S 12 = sign u y ln b b u y .
Integrating (11) 1 and exploiting (12), we find
u y = sign u y b 1 exp Π y .
Integrating once more with the no-slip boundary conditions u ( ± 1 ) = 0 , we have
u ( y ) = b 1 y 1 Π exp ( Π y ) exp ( Π ) ,
which provides the velocity profile in the whole domain. We may select the characteristic velocity U so that the velocity of the fluid in the centerline is 1. Imposing u ( 0 ) = 1 , we find
b = Π Π 1 + e Π Π = b b 1 + W b 1 b exp b 1 b ,
where W ( z ) is the zero branch of the Lambert function, see [18]. Relation (15) 1 is equivalent to setting
U = L b Δ p a 1 + exp ( Δ p a ) Δ p a .
As one can easily realize, the increase of the pressure drop Δ p is equivalent to the increase in the characteristic velocity U , but since the shear rate is limited, the velocity U is bounded from above (it cannot become arbitrarily large). Indeed, for Δ p ( 0 , ) , we have U ( 0 , L b ) , so that b > 1 , no matter how strong the applied pressure drop Δ p is. We notice that
U = Re ρ a , b = L b ρ a Re = Γ Re > 1 .
From the inequality in (17), we realize that the Reynolds number is bounded for this type of flow since from relation (16), the velocity U ( 0 , L b ) . In particular, Re ( 0 , Γ 2 ) and b = b ( Re ) are functions of the Reynolds number. Inserting (15) 1 into (14), and rearranging, we finally find the velocity field expressed in terms of the Reynolds number
u ( y ) = 1 Π y 1 + e Π y Π 1 + e Π , Π = b b 1 + W b 1 b exp b 1 b , b = Γ Re , Re ( 0 , Γ 2 ) .
In Figure 6, we show the velocity profile (18) for Re < Γ 2 , with Γ = 15 . We observe that the velocity profile assumes a wedge-like shape as Re approaches the limit value Γ 2 .
We observe that
Re = ρ a L 2 b 2 = Γ 2 Δ p a 1 + exp ( Δ p a ) Δ p a < 1 < Γ 2
Hence, Re is a monotonically increasing function of a . This is physically consistent since an increase in the value of a means a decrease in the apparent viscosity, i.e., less resistance to the flow and, hence, a larger value of Re , which is always bounded by Γ 2 . Notice that Γ 2 is a function of a . Hence, if a increases, Γ 2 also increases. For the other types of flow that will be studied in the next sections, we can make analogous considerations.

2.2. Flow between Coaxial Cylinder

In this section, we consider the flow between coaxial cylinders with inner and outer radii R 1 and R 2 , respectively. Without loss of generality, we may assume that the inner cylinder is fixed whereas the outer cylinder rotates anticlockwise with angular velocity Ω 2 . Distance is scaled with R 2 , velocity is rescaled with Ω 2 R 2 , and pressure is rescaled with a 1 . Exploiting cylindrical coordinates, we look for a steady solution of the type:
v = u ( r ) e θ , p = p ( r ) .
The only non-zero component of the deviatoric stress tensor is S r θ . Introducing the angular velocity ω ( r ) = u ( r ) / r , it is easy to see that
D = 1 2 r ω = r 2 ω ,
where = d / d r and where ω 0 , since the non-dimensional angular velocity w must be an increasing function of r ranging between 0 (inner cylinder) and 1 (outer cylinder). Therefore,
S r θ = ln b b r ω .
The mass balance is automatically satisfied, while the momentum balance is reduced to
Re u 2 r = p r 0 = r 2 S r θ = r 2 ln b b r ω .
Integrating the second of (22), we find
ω ( r ) = b 1 e c / r 2 r ,
with c being a positive constant to be determined. Integrating between R 1 = R 1 / R 2 < 1 and 1, we find (recall ω ( R 1 ) = 0 and ω ( 1 ) = 1 )
R 1 1 ω ( r ) d r = 1 = b ln R 1 1 2 E 1 ( c ) E 1 ( c / R 1 2 ) ,
where
E 1 ( ξ ) = ξ e z z d z ,
is the exponential integral. Equation (24) allows one to determine the constant c. Indeed let us rewrite (24) as
γ ( c ; R 1 ) = ln R 1 1 2 c c / R 1 2 e z z d z = 1 b .
The integrand function can be rewritten as
e z z = 1 z + n = 1 ( 1 ) n z n 1 n ! ,
so that
1 2 c c / R 1 2 e z z d z ln ( R 1 ) n = 1 ( 1 ) n 2 n ! n c n c R 1 2 n .
Hence, lim c 0 γ ( c ; R 1 ) = 0 and since
1 c 1 z R 1 2 c ,
we obtain
c c / R 1 2 e z z d z 1 c e c e c / R 1 2 > 0 0 , if c .
Therefore, lim c γ ( c ; R 1 ) = ln R 1 > 0 . Finally, we observe that γ / c > 0 , so that γ is a positive strictly increasing function of c bounded between 0 and ln R 1 > 0 . The existence of a solution to (24) (or equivalently (25)) is, thus, guaranteed if b 1 ( 0 , ln R 1 ) , i.e., if
Ω 2 < b ln R 2 R 1 .
Recalling that Re = ρ ( Ω 2 R 2 ) 2 a , we see that the existence of the solution is guaranteed if
Re < b R 2 ρ a ln R 2 R 1 = Γ ln R 2 R 1 = Γ ln R 1 .
In conclusion, the Reynolds number is limited for this type of flow. Once c is determined, the solution to the problem becomes
ω ( r ) = Γ Re ln r R 1 1 2 E 1 c r 2 E 1 c R 1 2 ,
where c is the unique solution of (24). From (23), we see that r ω < b , so the stress S r θ in (21) is well-defined. In Figure 7, we show the angular velocities ω ( r ) for different values of Re with Γ = 15 , R 1 = 0.1 , and with Re , satisfying (30).

2.3. Flow down an Inclined Plane

In this section, we focus on the flow down an inclined plane whose tilt angle with the horizontal direction is θ . We assume that the flow is uniform in the z direction so that v = ( u ( x , y ) , v ( x , y ) ) , and we denote with y = h ( x , t ) the upper free surface. We take a reference system { O , x , y } , where the x axis represents the bottom plane and y its normal direction. Rescaling as in the previous section, and assuming that the flow is driven by gravity, the motion equation is
Re v ˙ = p + div ( S ) + Re Fr 2 sin θ i cos θ j ,
where
Fr = U g L
is the Froude number. We suppose that the fluid is surrounded by a motionless ambient gas so that, neglecting the surface tension on h ( x , t ) , we can write the non-dimensional free boundary conditions as
h x 2 p S 11 + 2 h x S 12 + p S 22 1 + h x 2 = 0 , on y = h ,
S 12 h x 2 1 + S 11 S 22 h x = 0 , on y = h ,
v = h x u + h t , on y = h .
The first two relations express the continuity of the normal and tangential stresses, respectively (dynamic conditions), while the last is a consequence of the fact that h ( x , t ) is a material surface (kinematic condition). We look for a solution of the types v = u ( y ) i and h = c o n s t , so that S 12 = S 12 ( y ) , S 11 = S 22 = 0 and
0 = p x + S 12 y + λ sin θ , 0 = p y λ cos θ , λ = Re Fr 2 = ρ g L a .
The constant λ is a material parameter. The dynamic boundary conditions for the stress reduce to p = 0 (the external pressure is rescaled to zero), S 12 = 0 on y = h , while the kinematic boundary condition is automatically satisfied. The no-slip conditions on the bottom plane are u = 0 and 2 D = u y > 0 . Integration of linear momentum yields
p = λ cos θ h y , S 12 = ln b b u y = δ h y .
where δ = λ sin θ > 0 . After some calculation, we find
u y = b 1 e δ ( h y ) < b ,
and
u ( y ) = b y e δ ( h y ) e δ h δ
The thickness h of the fluid is unknown and can be determined by imposing the inlet non-dimensional flux Q = Q / ( U H ) , where Q is the dimensional discharge. Exploiting (38), we have
Q = 0 h u ( y ) d y = b h 2 2 1 e δ h ( 1 + δ h ) δ 2
Hence, Q = Q ( h ) with Q ( 0 ) = 0 , Q ( ) = , Q = h ( 1 e δ h ) 0 . As a consequence, for any positive given Q, there exists a unique h, satisfying (39). We may select U , so that the maximum velocity, which is attained at y = h , is equal to one. From (38), setting u ( h ) = 1 , we find
b = h 1 e δ h δ 1 h = 1 δ 1 + δ b + W e 1 + δ b ,
where W is again the zero branch of the Lambert function. Substituting (40) 1 into (38), we find
u ( y ) = e δ h + δ y e δ ( h y ) e δ h + δ h 1 .
Recalling that b = Γ / Re , we finally write the velocity profile for the flow down an inclined plane as
u ( y ) = e δ h + δ y e δ ( h y ) e δ h + δ h 1 , h = 1 δ 1 + δ Re Γ + W e 1 + δ Re Γ
In Figure 8 we show the velocity profile (42) for various Re with Γ = 15 , λ = 1.3 and θ = π / 4 .
We observe that, in this case, no limitation on Re is present. Indeed, from (40), we see that the increase of Re results in an increase in h, i.e., in practice, the height of the fluid layer can adapt to the imposed discharge.

3. Linear Stability Analysis

In this section, we consider two-dimensional perturbations of the basic solutions determined in the previous sections. In particular, let ψ b denote any of the steady solutions found in the previous sections. We write
ψ = ψ b + ψ ^ ( y ) e i α ( x c t )
so that ψ is the sum of the basic function ψ b and the perturbation ψ ^ ( y ) e i α ( x c t ) , where ψ ^ is the complex wave amplitude, α is the wavenumber, c is the complex wave velocity, and we assume that ψ ^ = O ( ε ) , ε 1 . The stress tensor can be decomposed as
S b + S ^ e i α ( x c t ) = ln b b 2 D b + D ^ e i α ( x c t ) D b + D ^ e i α ( x c t ) D b + D ^ e i α ( x c t ) .
The above can be linearized (neglecting terms of order ε and higher) to obtain
S ^ = ln b b 2 D b D ^ D b ( D b · D ^ ) D b 2 D b 3 + ( D b · D ^ ) D b D b 2 b 2 D b
so that S ^ is linear in D ^ .

3.1. Channel Flow

In the two-dimensional setting, we have
D ^ 11 = D ^ 22 = v ^ , D ^ 12 = 1 2 i α v ^ + α 2 v ^
In the context of channel flow, we find
S ^ 11 = S ^ 22 = 1 u b ln b b u b T ( y ) 2 v ^
S ^ 12 = 1 b u b M ( y ) ( v ^ + α 2 v ^ ) i α ,
where we indicated with u b ( y ) the basic solution (18). The perturbed momentum equation is
Re [ ( u b c ) v ^ α 2 v ^ ) u b v ^ = 2 i α S ^ 11 + S ^ 12 + α 2 S ^ 12 .
Substituting (47), (48) into (49), we obtain the perturbation equation for v ^
i α Re [ ( u b c ) D 2 α 2 ) u b v ^ = 4 α 2 D T D v ^
D 2 M D 2 + α 2 v ^ α 2 M D 2 + α 2 v ^
where M ( y ) and T ( y ) are defined in Equations (47) and (48), respectively. The boundary conditions are v ^ ( ± 1 ) = D v ^ ( ± 1 ) = 0 . Equation (50) is solved via a spectral collocation method based on Chebyshev polynomials. Notice that because of the particular constitutive equation considered here, the classical Orr–Sommerfeld equation cannot be recovered from (50). In Figure 9, we plot the neutral stability curves for some values of Γ , while in Figure 10, we display the critical Reynolds number Re c as a function of the parameter Γ . Notice that the critical Reynolds number is always below the limit Γ 2 , which is the upper bound of Re , see Section 2.1. In particular, we observe that the deviation of Re c from the limit value Re = Γ 2 increases with Γ .

3.2. Taylor–Couette Flow

Here, we study the linear stability of the flow between coaxial rotating cylinders. In this case, the perturbed velocity is v = v b ( r ) e θ + v ^ ( r ) e i α ( z c t ) , where
v ^ = u ^ e r + v ^ e θ + w ^ e z
is the complex amplitude of the perturbation of the velocity, where v b ( r ) = r ω b ( r ) is the basic tangential velocity and where ω b ( r ) is the angular velocity given in (31). We are, thus, considering an axisymmetric perturbation (there is no dependence on the azimuthal coordinate θ ). The components of the symmetric part of the perturbation velocity gradient are
D ^ r r = u ^ , D ^ θ θ = u ^ r , D ^ z z = i α w ^ .
D ^ r θ = 1 2 v ^ v ^ r , D ^ r z = 1 2 i α u ^ + w ^ , D ^ θ z = i α v ^ 2 ,
where the prime now stands for differentiation with respect to r, i.e., = d / d r . The linearized momentum equation for the perturbations ( u ^ , v ^ , w ^ ) are
Re i α c u ^ 2 ω b v ^ = Π p ^ + r S ^ r r r S ^ θ θ r + i α S ^ r z
Re i α c v ^ + ( r 2 ω b ) r u ^ = ( r 2 S ^ r θ ) r 2 + i α S ^ θ z
Re i α c w ^ = Π i α p ^ + ( r S ^ r z ) r + i α S ^ z z
Multiplying (53) by i α , differentiating (55) with respect to r, and summing up, we have
Re i α 2 ω b v ^ c w ^ α 2 c u ^ = i α r r S ^ r r
1 r 2 r S ^ r z + 1 r r S ^ r z + i α S ^ θ θ r + S ^ z z i α S ^ r z
From the continuity equation, we find
w ^ = i α u ^ + u ^ r , w ^ = i α u ^ u ^ r 2 + u ^ r
The linearized constitutive Equation (45) in the cylindrical coordinates becomes
S ^ = ln b b r ω b 2 D ^ r ω b 4 D ^ r θ D b ( r ω b ) 2 + 4 D ^ r θ D b r ω b ( b r ω b )
from which we derive
S ^ r r = T ( r ) 2 u ^ , S ^ θ θ = T ( r ) 2 u ^ r , S ^ z z = T ( r ) 2 u ^ + u ^ r
S ^ r θ = M ( r ) v ^ v ^ r , S ^ r z = T ( r ) i α u ^ u ^ r 2 + u ^ r + α 2 u ^ ,
S ^ θ z = T ( r ) i α v ^ , S ^ z z = T ( r ) 2 u ^ + u ^ r ,
where
T ( r ) = 1 r ω b ln b b r ω b , M ( r ) = 1 b r ω b .
Inserting (58)–(60) into (54), (56), we finally find the eigenvalue problem
c Re D 2 + D r 1 r 2 α 2 0 0 Re i α u ^ v ^ = L 1 M 1 L 2 M 2 u ^ v ^
where D = d / d r and where
L 1 = i T 2 α r 2 + 3 α r 4 α 3 + T α r 3 α r 3 + T 1 α r 2 α +
i T 2 α r 3 α r 3 + T 3 α r 2 + 2 α T 1 α r D +
i T 3 α r 2 + 2 α T 3 α r T 1 α D 2 +
i T 2 α r + T 2 α D 3 T 1 α D 4
M 1 = 2 Re i α ω b , L 2 = Re r ω b + 2 ω b
M 2 = M r + M r 2 + α 2 r 2 T M r + M D M D 2
The eigenvalue problem is coupled with boundary conditions u ^ ( 1 ) = D u ^ ( 1 ) = v ^ ( 1 ) = 0 and u ^ ( R 1 ) = D u ^ ( R 1 ) = v ^ ( R 1 ) = 0 . Problem (61) is solved by a spectral collocation method based on Chebyshev polynomials (the spatial domain [ R 1 , 1 ] is mapped in the interval [ 1 , 1 ] ). In Figure 11, we plot the neutral stability curves for some values of Γ . In Figure 12, we show the critical Reynolds number as a function of the parameter Γ . Notice that the critical Reynolds number is always below the limit ( Γ ln R 1 1 ) 2 . In particular, we observe that the deviation of Re c from the limit value Re = ( Γ ln R 1 1 ) 2 increases with Γ as in the channel flow.

3.3. Flow down an Inclined Plane

In this section, we study the linear stability of the downhill flow presented in Section 2.3. The perturbed variables are in the form (43) and the equation for the perturbation is still (50), the spatial domain being [ 0 , h b ] , with h b given in (42) 2 . Introducing the perturbed free surface
h = h b + h ^ e i α ( x c t ) , h ^ C ,
we check that the linearized boundary conditions (32)–(34) become:
h ^ = 2 b 1 v ^ ( h b ) p ^ ( h b ) λ cos θ ,
h ^ = 1 δ b v ^ ( h b ) + α 2 v ^ ( h b ) i α ,
v ^ ( h b ) = i α h ^ 1 c .
The pressure term in (62) can be written in terms of v ^ exploiting the first component of the linearized momentum equation for the perturbation, namely
Re v ^ ( c u b ) + v ^ u b = i α p ^ + i α S ^ 11 + S ^ 12 ,
with S ^ 11 , S ^ 12 given by (47), (48). Using (65), and eliminating h ^ from (62), (63) we find that
v ^ + α 2 v ^ 1 c + δ b v ^ = 0 , y = h b ,
v ^ + α 2 v ^ v ^ + α 2 v ^ i α cot θ + δ
v ^ 4 α 2 + i α b Re ( 1 c ) = 0 , y = h b .
Conditions (66), (67), together with the no-slip conditions
v ^ ( 0 ) = v ^ ( 0 ) = 0 ,
provide the boundary conditions for the perturbation v ^ . In conclusion, the eigenvalue problem for the stability of the flow down an inclined plane is given by Equation (50)—in the domain [ 0 , h b ] -and the boundary conditions (66)–(68). The complex values of the parameter c that produce a non-trivial solution give the relation between the basic flow and the evolution of the disturbance modes. In particular, when the growth rate of the disturbance σ = α I m ( c ) is positive, the amplitude of the disturbance becomes unbounded and the system is unstable. Here, we limit our analysis to the so-called long-wave analysis, i.e., we assume that α 1 , so that the solution can be sought in the form of expansion in powers of α ,
v ^ = v ^ o + α v ^ 1 + α 2 v ^ 2 + . . . . . . , c = c o + α c 1 + α 2 c 2 + . . . . . . .
Inserting the expansion above in the eigenvalue problem and grouping the terms of the same order, we obtain a sequence of eigenvalue problems that can be solved iteratively. In particular, it can be proved that the problem at the zero order is given by
M ( y ) v ^ o = 0 , v ^ o = v ^ o = 0 , y = 0 v ^ o ( 1 c o ) + δ b v ^ o = 0 , y = h b , v ^ o δ v ^ o = 0 , y = h b .
Problem (69) 1 is solved by
v ^ o = e δ y δ y 1 , c o = e δ h b ( δ + b ) e δ h b b ( δ h b + 1 ) δ
from which we immediately see that c o R . We notice that the general eigenvalue problems (50), (66), (67), (68) are homogeneous, so that the solution is defined up to a multiplicative constant. Following the procedure used in [19], we select such a constant so that v ^ ( h b ) = v ^ o ( h b ) , implying v ^ i ( h b ) = 0 for i = 1 , 2 , 3 , . . . . . The problem with the first order is, therefore,
i Re ( u b c o ) v ^ o u b v ^ o = M ( y ) v ^ 1 , v ^ 1 = v ^ 1 = 1 , y = 0 v ^ 1 ( 1 c o ) + δ b v ^ 1 c 1 v ^ o = 0 , y = h b , v ^ 1 δ v ^ 1 + i Re v o ( c o 1 ) i cot θ v ^ o = 0 , y = h b , v ^ 1 = 0 , y = h b .
Integration of (71) 1 gives functions v ^ 1 = v ^ 1 ( y ; A , B , C , D ) in which the dependencies on the integration constants A , B , C , D are linear. The impositions of the five BCs (71) 2 –(71) 5 lead to a linear system whose solution provides the constants A , B , C , D , and the eigenvalue c 1 . The calculation for the determination of the eigenvalue c 1 was checked using the symbolic software wx Maxima. We obtain
c 1 = i ( c o 1 ) δ 3 ( B 2 h b + B 1 B 5 ) 2 B 4 + b δ B 5 δ 3 δ 3 ( e δ h δ h 1 ) +
+ i ( c o 1 ) B 4 ( δ h b ) 2 + 2 ( δ h b + 1 ) + δ ( 1 + δ h b ) ( B 3 b B 5 δ ) δ 3 e δ h ( e δ h δ h 1 ) C
where B i , i = 1 , 2 , 3 , 4 , 5 are functions of the parameters b, h b , δ , Re , cot θ , c o . In particular, B i is a pure real number (the imaginary part is zero), so that c 1 is a pure complex number. For the sake of brevity, we did not report their expressions, which are quite lengthy. In the long wavelength expansion, of the first order, we take c c o + α c 1 , so that the growth rate of the amplitude disturbance is
σ = i α c 1 .
Recalling that h b = h b ( Re , Γ , δ ) , b = b ( h b , δ ) (see (40)), δ = λ sin θ we immediately realize that at the first order
B i = B i ( Re , Γ , θ , λ ) , c 1 ( R e , Γ , θ , λ ) .
and
σ = σ ( α , R e , Γ , θ , λ ) = i α c 1 ( R e , Γ , θ , λ ) .
In Figure 13, we plot the approximation of the marginal stability curves (valid for small α ) σ = 0 for θ = π / 4 , λ = 1 and with Γ = 10 , 15 , 20 , 25 , 30 . These are clearly vertical straight lines in the plane ( Re , α ) that provide the critical Reynolds above which we have instability. The dependency of the critical Reynolds number Re c as a function of Γ is depicted in Figure 14.
As observed in Section 2.3, here we have no limitation of the Reynolds number, since the quantity b = Γ / Re is not limited. We notice that the critical Reynolds for the flow down an incline is, in general, much smaller than in the other two cases considered earlier. This is due to the fact that the channel flow and the flow between concentric cylinders are wall-bounded, while the flow down an incline has a free surface. Indeed, it is well known that in the wall-bounded flow, we may have the development of boundary layers that tend to stabilize the flow, see [20]. Furthermore, the discrepancy between the critical Reynolds number of the channel flow (or Taylor–Couette flow) and the flow down an incline is evident even in the case of a Newtonian fluid. For instance, in the channel flow of a Newtonian fluid driven by a pressure drop Re c 5700 (and in the free flow down an incline), the critical Reynolds is given by the relation, see [21]
Re c = 5 4 cot θ ,
for an angle of π / 4 (as the one I have considered in the paper) gives Re c = 1.25 . Finally, we remark that the critical Reynolds numbers for the free surface flow are in the same range as those obtained in [22].

4. Conclusions and Perspectives

We investigated a simple shear flow of a class of dilatant fluids in which the modulus of the deviatoric stress diverges as the modulus of the shear rate tensor tends to a finite positive value. For convenience, these fluids have been termed dilatant fluids with a limited shear rate. The constitutive equation is given by a logarithmic function with two material parameters, 1 / a and b , representing the characteristic apparent viscosity and the shear rate threshold, respectively. The motivation for choosing such a specific form for the apparent viscosity is that it fits (for appropriate values of a and b ) the flow curves of granular suspensions. The adopted constitutive equation can be viewed as a regularization of a discontinuous shear-thickening constitutive law, such as the one depicted in Figure 1b. Despite the nonlinearity of the apparent viscosity, we are able to determine steady-state solutions for some simple shear flow types. In particular, we considered the (i) planar channel flow, (ii) the flow between concentric cylinders, and (ii) the flow down an inclined plane. In the first two cases, we have shown that the Reynolds number cannot exceed a certain threshold that depends on the material parameters. This seemingly weird behavior is explained; the increase in the velocity automatically increases the resistance opposed to the flow by the increasing viscosity. In this mechanism, the acceptable velocities experienced by the fluid cannot exceed a fixed threshold. In particular, we proved that the limit threshold for the Reynolds number is an increasing function of a and b . To investigate the response of the basic steady-state flow to disturbances, we also performed a two-dimensional modal stability analysis for each of the three simple flow types considered. In cases (i) and (ii), the corresponding eigenvalue problems were solved by means of a spectral collocation method based on Chebyshev polynomials, allowing us to plot the neutral stability curves and detect the critical Reynolds number. For the case of the flow down an inclined plane, we limited our analysis to long-wave disturbances (i.e., small wave numbers), explicitly determining the critical Reynolds number. For the three simple shear flow types considered, we found that the critical Reynolds number is an increasing function of a and b . This proves that when the limit shear rate threshold is increased, or when the apparent viscosity is decreased, the fluid can be stable for a larger range of Reynolds numbers.
A natural continuation of this work would be (i) the extension to some other simple shear flow; (ii) the analysis of the flow in small-aspect ratio geometries (e.g., lubrication flow); (iii) the investigation of the effects of three-dimensional disturbances and the verification of the validity of Squire’s theorem, which is not automatically valid for non-Newtonian flow. All of these issues are currently under investigation and will be the subject of a forthcoming paper.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

This research was partially supported by the GNFM of the Italian INDAM.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Metzner, A.; Whitlock, M. Flow behavior of concentrated (dilatant) suspensions. Trans. Soc. Rheol. 1958, 2, 239–254. [Google Scholar] [CrossRef]
  2. Sidaoui, N.; Arenas Fernandez, P.; Bossis, G.; Volkova, O.; Meloussi, M.; Aguib, S.; Kuzhir, P. Discontinuous shear thickening in concentrated mixtures of isotropic-shaped and rod-like particles tested through mixer type rheometry. J. Rheol. 2020, 64, 817–836. [Google Scholar] [CrossRef]
  3. Thomas, J.E.; Goyal, A.; Singh Bedi, D.; Singh, A.; Del Gado, E.; Chakraborty, B. Investigating the nature of discontinuous shear thickening: Beyond a mean-field description. J. Rheol. 2020, 64, 329–341. [Google Scholar] [CrossRef] [Green Version]
  4. Farina, A.; Fasano, A.; Fusi, L.; Rajagopal, K. The one-dimensional flow of a fluid with limited strain-rate. Q. Appl. Math. 2011, 69, 549–568. [Google Scholar] [CrossRef] [Green Version]
  5. Mari, R.; Seto, R.; Morris, J.F.; Denn, M.M. Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions. J. Rheol. 2014, 58, 1693–1724. [Google Scholar] [CrossRef] [Green Version]
  6. Clavaud, C.; Bérut, A.; Metzger, B.; Forterre, Y. Revealing the frictional transition in shear-thickening suspensions. Proc. Natl. Acad. Sci. USA 2017, 114, 5147–5152. [Google Scholar] [CrossRef] [Green Version]
  7. Singh, A.; Mari, R.; Denn, M.M.; Morris, J.F. A constitutive model for simple shear of dense frictional suspensions. J. Rheol. 2018, 62, 457–468. [Google Scholar] [CrossRef]
  8. Wyart, M.; Cates, M.E. Discontinuous shear thickening without inertia in dense non-Brownian suspensions. Phys. Rev. Lett. 2014, 112, 098302. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Brown, E.; Zhang, H.; Forman, N.A.; Maynor, B.W.; Betts, D.E.; DeSimone, J.M.; Jaeger, H.M. Shear thickening and jamming in densely packed suspensions of different particle shapes. Phys. Rev. E 2011, 84, 031408. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Chhabra, R.; Richardson, J. Non-Newtonian Flow in Process Industries 1999; Biddles Ltd.: Guildford, UK; King’s Lynn, UK, 1999. [Google Scholar]
  11. Ancey, C. Plasticity and geophysical flow: A review. J. Non-Newton. Fluid Mech. 2007, 142, 4–35. [Google Scholar] [CrossRef]
  12. Lee, Y.S.; Wetzel, E.D.; Wagner, N.J. The ballistic impact characteristics of Kevlar® woven fabrics impregnated with a colloidal shear thickening fluid. J. Mater. Sci. 2003, 38, 2825–2833. [Google Scholar] [CrossRef]
  13. Gálvez, L.O.; de Beer, S.; van der Meer, D.; Pons, A. Dramatic effect of fluid chemistry on cornstarch suspensions: Linking particle interactions to macroscopic rheology. Phys. Rev. E 2017, 95, 030602. [Google Scholar]
  14. Blechta, J.; Málek, J.; Rajagopal, K. On the classification of incompressible fluids and a mathematical analysis of the equations that govern their motion. SIAM J. Math. Anal. 2020, 52, 1232–1289. [Google Scholar] [CrossRef] [Green Version]
  15. Rajagopal, K.R. On implicit constitutive theories. Appl. Math. 2003, 48, 279–319. [Google Scholar] [CrossRef] [Green Version]
  16. Ozturk, D.; Morgan, M.L.; Sandnes, B. Flow-to-fracture transition and pattern formation in a discontinuous shear thickening fluid. Commun. Phys. 2020, 3, 119. [Google Scholar] [CrossRef]
  17. Fusi, L.; Saccomandi, G.; Rajagopal, K.R.; Vergori, L. Flow past a porous plate of non-Newtonian fluids with implicit shear stress shear rate relationships. Eur. J. Mech.-B/Fluids 2022, 92, 166–173. [Google Scholar] [CrossRef]
  18. Corless, R.M.; Gonnet, G.H.; Hare, D.E.; Jeffrey, D.J.; Knuth, D.E. On the LambertW function. Adv. Comput. Math. 1996, 5, 329–359. [Google Scholar] [CrossRef]
  19. Pascal, J. Linear stability of fluid flow down a porous inclined plane. J. Phys. D Appl. Phys. 1999, 32, 417. [Google Scholar] [CrossRef]
  20. Nouar, C.; Frigaard, I. Stability of plane Couette–Poiseuille flow of shear-thinning fluid. Phys. Fluids 2009, 21, 064104. [Google Scholar] [CrossRef]
  21. Yih, C.S. Stability of liquid flow down an inclined plane. Phys. Fluids 1963, 6, 321–334. [Google Scholar] [CrossRef] [Green Version]
  22. Darbois Texier, B.; Lhuissier, H.; Forterre, Y.; Metzger, B. Surface-wave instability without inertia in shear-thickening suspensions. Commun. Phys. 2020, 3, 232. [Google Scholar] [CrossRef]
Figure 1. String-dashpot system (a) and constitutive law (b).
Figure 1. String-dashpot system (a) and constitutive law (b).
Fluids 08 00025 g001
Figure 4. Flow curve (shear stress vs. shear rate) for TiO 2 suspension (47%), see Figure 2 in [1]. Experimental data (stars) and fitting curve (continuous). Fitting parameters a = 1.5 (ft 2 /lbr), b = 10 2 s 1 .
Figure 4. Flow curve (shear stress vs. shear rate) for TiO 2 suspension (47%), see Figure 2 in [1]. Experimental data (stars) and fitting curve (continuous). Fitting parameters a = 1.5 (ft 2 /lbr), b = 10 2 s 1 .
Fluids 08 00025 g004
Figure 5. Flow curve (shear stress vs. shear rate) for the cornstarch suspension (50%), see Figure 1 in [16]. Experimental data (stars) and fitting curve (continuous). Fitting parameters a = 0.2 (Pa 1 ), b = 18 s 1 .
Figure 5. Flow curve (shear stress vs. shear rate) for the cornstarch suspension (50%), see Figure 1 in [16]. Experimental data (stars) and fitting curve (continuous). Fitting parameters a = 0.2 (Pa 1 ), b = 18 s 1 .
Fluids 08 00025 g005
Figure 6. Velocity profiles (18) for various Re < Γ 2 , Γ = 15 .
Figure 6. Velocity profiles (18) for various Re < Γ 2 , Γ = 15 .
Fluids 08 00025 g006
Figure 7. Angular velocity (31) for various Re < Γ 2 , Γ = 15 , R 1 = 0.1 .
Figure 7. Angular velocity (31) for various Re < Γ 2 , Γ = 15 , R 1 = 0.1 .
Fluids 08 00025 g007
Figure 8. Velocity (42) for various Re , Γ = 15 , λ = 1.3 , θ = π / 4 . Flow down an incline.
Figure 8. Velocity (42) for various Re , Γ = 15 , λ = 1.3 , θ = π / 4 . Flow down an incline.
Fluids 08 00025 g008
Figure 9. Neutral stability curves in the ( α , Re ) plane, channel flow. Γ = 10 , 15 , 20 , 25 , 30 .
Figure 9. Neutral stability curves in the ( α , Re ) plane, channel flow. Γ = 10 , 15 , 20 , 25 , 30 .
Fluids 08 00025 g009
Figure 10. Critical Reynolds number vs. Γ : channel flow. Notice that the discrepancy between the Re c and the upper bound Γ 2 increases with Γ .
Figure 10. Critical Reynolds number vs. Γ : channel flow. Notice that the discrepancy between the Re c and the upper bound Γ 2 increases with Γ .
Fluids 08 00025 g010
Figure 11. Neutral stability curves, the flow between rotating cylinders. Γ = 10 , 15 , 20 , 25 , 30 .
Figure 11. Neutral stability curves, the flow between rotating cylinders. Γ = 10 , 15 , 20 , 25 , 30 .
Fluids 08 00025 g011
Figure 12. Critical Reynolds number vs. Γ : flow between rotating cylinders. Notice again the discrepancy between Re c and the upper bound as Γ increases.
Figure 12. Critical Reynolds number vs. Γ : flow between rotating cylinders. Notice again the discrepancy between Re c and the upper bound as Γ increases.
Fluids 08 00025 g012
Figure 13. First-order approximation of the neutral stability curves, downhill flow. Γ = 10 , 15 , 20 , 25 , 30 , θ = π / 4 , λ = 1 .
Figure 13. First-order approximation of the neutral stability curves, downhill flow. Γ = 10 , 15 , 20 , 25 , 30 , θ = π / 4 , λ = 1 .
Fluids 08 00025 g013
Figure 14. Critical Reynolds number vs. the Γ downhill flow. θ = π / 4 , λ = 1 .
Figure 14. Critical Reynolds number vs. the Γ downhill flow. θ = π / 4 , λ = 1 .
Fluids 08 00025 g014
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Fusi, L. Shear Flows of Dilatant Fluids with Limited Shear Rates: Analytical Results and Linear Stability Analysis. Fluids 2023, 8, 25. https://doi.org/10.3390/fluids8010025

AMA Style

Fusi L. Shear Flows of Dilatant Fluids with Limited Shear Rates: Analytical Results and Linear Stability Analysis. Fluids. 2023; 8(1):25. https://doi.org/10.3390/fluids8010025

Chicago/Turabian Style

Fusi, Lorenzo. 2023. "Shear Flows of Dilatant Fluids with Limited Shear Rates: Analytical Results and Linear Stability Analysis" Fluids 8, no. 1: 25. https://doi.org/10.3390/fluids8010025

APA Style

Fusi, L. (2023). Shear Flows of Dilatant Fluids with Limited Shear Rates: Analytical Results and Linear Stability Analysis. Fluids, 8(1), 25. https://doi.org/10.3390/fluids8010025

Article Metrics

Back to TopTop