Flow of a limited stress fluid between plates rotating about different axes

In this paper we study the flow of a particular type of non-Newtonian fluids generated by the rotation of parallel infinite plates about distinct axes. The constitutive law of these fluids mimics the response of a class of seemingly viscoplastic or “yield stress” materials in which the norm of the stress is bounded by a critical threshold (limited stress fluids). We assume that the plates rotate with the same angular velocity and we show that, in this case, the mathematical problem can be reduced to solving a BVP where the unknowns are the coordinates of the center of rotation. We solve the problem numerically (by means of a spectral collocation method), and we investigate the dependence of the locus of the center of rotation on the material parameters. We prove that, even for small Reynolds numbers, the flow may exhibit boundary layers depending on the particular choice of the parameters.


Introduction
The flow of fluids between rotating plates has received considerable attention since the seminal work of Von Karman [1] who studied the axisymmetric flow due to two infinite parallel discs rotating about a common axis.Here we are interested in studying the flow between infinite plates rotating about distinct axes at constant identical angular velocity.Such flows have relevance to those that take place in an instrument that is called orthogonal rheometer (see Maxwell and Chartoff [2]), a device that is widely used to determine the rheological properties of various types of fluids by measuring the forces and torque that act on the plates.
The class of flows engendered belongs to the class of pseudo-planar flows that were first studied systematically by Berker [3].The flows also fall into the class of motions with constant principal relative stretch histories introduced by Noll [4] (see also Rajagopal [5]).
Abbot and Walters [6] were the first authors to derive an exact solution of the flow in the case of a Navier-Stokes fluid, and this has been followed by studies on the flows of various non-Newtonian fluids (see Goldstein [7], Rajagopal [8], Rajagopal and Gupta [8], Bower et al. [9], Rajagopal and Wineman [10], Zhang and Goddard [11], Al Khatib [12], Fusi et al. [13] and others).The paper by Rajagopal [14] provides an exhaustive list of references for the problem under investigation.
Here we are interested in studying the flow in an orthogonal rheometer of a fluid whose constitutive law mimics that of a seemingly viscoplastic material.This particular class of fluids has been introduced first in Garimella et al. [15], where the fully developed flow in a cylindrical pipe was investigated.
The simplest viscoplastic model is the Bingham fluid that behaves like a rigid body when the applied stress is below a critical threshold and flows like a Navier-Stokes fluid when the threshold is exceeded (see Huilgol [16]).In the Bingham model the slope of the shear rate1 γ * vs shear stress τ * is infinite for zero shear rate, rendering the constitutive law a graph and not a function.To smooth this particular behavior (that can be a source of serious problems especially when performing multi-dimensional numerical simulations) Garimella et al. [15] have developed a new model, characterized by two material constant α * 1 , α * 2 , in which the slope of the stress strain relation can become arbitrarily large in a neighborhood of γ * = 0. Differently from the Bingham model, in this new model the constitutive law is an actual function and the stress is always determined for every value of the strain rate.Moreover, they assume that, for large strain rates, the shear stress tends to the asymptotic limiting value α * 1 that cannot be overcome (limited stress fluids).The model is an attempt to trying to approximate "yield stress" materials and seems to well describe the behavior of some viscoplastic materials such as kaolin-water solutions, foams and paints.We remark, however, that the response of what one means by "yield stress" is quite subjective and also depends on the time scale, length scale and force scale of observation, and it is further complicated as there are really no reliable data at very low shear rates.From this perspective, we may say that how well a particular fluid constitutive relation fits the data close to zero shear rate is, to some extent, in the eye of the modeler.
The paper is organized as follows: in Sect. 2 we present the mathematical model and we write the balance equations that govern the flow.The problem is initially written for dimensional variables and subsequently put in a non-dimensional form.We show that, under appropriate hypotheses, the problem can be reduced to a nonlinear boundary value problem.In Sect. 3 we perform an asymptotic analysis that allows one to get an explicit solution for the particular case of a small Reynolds number or small aspect ratio.In Sect. 4 we show how to solve numerically the BVP by means of a spectral collocation method based on Chebyshev polynomials.In Sect. 5 we discuss the results of the numerical simulations.Finally in Sect.6 we make some final considerations.

The mathematical model
We consider a class of incompressible fluids that are determined by the constitutive relation T * = −p * I + S * , where where D * is the symmetric part of the velocity gradient, − p * I is the indeterminate part of the stress due to the constraint of incompressibility, and α * 1 and α * 2 are positive constant material parameters.The dimension of α * 1 is that of a pressure, whereas the dimension of α * 2 is that of a time.Constitutive relation ( 1) is a special case of the classical regularized Bingham model introduced by Papanastasiou in [17].We also notice that model (1) can be included in the classification of incompressible fluids presented in [18].The sketch of the system is depicted in Fig. 1.The velocity field is of the form v * = u * e 1 + v * e 2 + w * e 3 .The governing equations are where ρ * is the constant density of the fluid and g * is the acceleration due to gravity.The upper plate, placed at z * = h * , rotates around the axis (0, a * ) with angular velocity * .The lower plate, placed at z * = −h * , rotates around (0, −a * ) with the same angular velocity.Assuming no-slip at the lower and upper plate we record the following boundary conditions We rescale the system as follows The non-dimensional constitutive equation becomes where The non-dimensional governing equations are where are the Reynolds number, the aspect ratio and the ratio between gravitational and viscous forces, respectively.The components of the traceless symmetric tensor D are Rajagopal [5] proved that the equations of motion for a simple fluid are compatible with the following ansatz for the velocity field: where f (z) and g(z) are functions to be determined.It is straightforward to verify that the incompressibility constraint is automatically satisfied by (9) and that Because of (9) the balance of linear momentum reduces to Integrating the last equation of ( 11) we find that p = p o (x, y) − φz with p o (x, y) unknown.Substituting The constants C 1 , C 2 are nonzero only if an horizontal pressure gradient is applied.For the sake of simplicity, here we assume where η( γ ) and γ are given in ( 6) and (10), respectively.Recalling (4), ( 5), the boundary conditions for f, g become After some manipulations, problem (13) can be rewritten as where ϕ = 2ξ 2 Re and where we have re-defined the apparent viscosity in the following way: On eliminating g from the first of ( 15) and f from the second of ( 15) we end up with: Our problem is thus reduced to solving second-order nonlinear BVP (18) with BC (14).We notice that the coefficient of the second derivatives in ( 18) is always strictly positive so that the system does not degenerate.Indeed, one can show that the coefficient is such that Therefore, the system degenerates only if θ = ∞, i.e., only if γ = ∞.We observe that the only nonzero components of the stress tensor are i.e., the shear stress components.The modulus of the stress becomes Therefore the system degenerates only if the modulus of the stress becomes 1, which is impossible because the upper bound for the stress can be reached only if γ = ∞.

Asymptotic analysis at low ϕ
We begin our analysis showing that, when ϕ is small (this can happen either for small Re or for small aspect ratio ξ ), we can look for a solution that can be written as a power series of ϕ, i.e., The boundary conditions are In order to determine the differential problems at the various orders it is convenient to consider the problem in form (13), i.e., ⎧ ⎨ ⎩ where we recall that Inserting ( 22) into ( 23) and grouping the terms of the same order we get a series of linear problems that can be solved iteratively.Here, we shall consider the zero, first and second order only.Higher approximations can of course be obtained if one desires more accuracy.It is easy to see that, at the leading order, we get At the first order the problem is given by whose integration leads to At the second order the system is yielding Therefore, in conclusion, at the second order our approximated solution is We shall use (27) for comparison with the numerical solution in Sect. 5

Numerics
In order to solve nonlinear system (18) numerically we make use of a spectral collocation method based on Chebyshev polynomials.System (18) is discretized using Chebyshev differentiation matrices, whose definition can be found in [19], [20].The space variable z is discretized on a grid of N + 1 Gauss-Lobatto points z k = cos(kπ/N ), k = 0, . . ., N clustered at the boundaries z = ±1.The functions f , g together with their first and second derivatives are discretized by where and where D is the (N + 1) × (N + 1) differentiation matrix with entries (see [19], [20]) where c i = 0 if i = 0, N and c i = 1 otherwise.The discretized formulation of ( 18) becomes where the quantities , μ, • μ are the vectors of components and where is the Hadamard product (element-wise or Schur product).The discretized system is thus a nonlinear algebraic system of the type We solve system (29) via a Newton-Rapson method in which at each step we calculate the Jacobian of (29).The entries of the Jacobian are reported in appendix.The boundary conditions are implemented at each step of the Newton-Rapson scheme, and the iteration procedure is terminated when an accuracy of 10 −6 is reached.

Results
In this section we show the results of the numerical simulations of system (29).In Fig. 2 we show the x−component of the center of rotation x = f (z) for α = 0.2, β = 0.3 and Re ranging from 10 −2 to 10.We recall that ϕ = (2β 2 ) −1 Re.In Fig. 3 we display the y−component of the center of rotation y = g(z) for the same values of the parameters.We observe that, as Re increases, the center of rotation coincides with the z−axis except for boundary layers adjacent to the rotating plates.In practice, as Re is increased, the fluids rotate around the z−axis practically everywhere in the gap with the exception of the boundary layers in which the center of rotation ( f, g) rapidly moves toward the points (0, ±1), where the latter represent the centers of rotation of the infinite rotating plates.In Fig. 4 we display shear stress modulus τ (21) for increasing Re.For small Reynolds τ is essentially constant and equal to 1 − e −αβ .Indeed, recalling ( 21) and (24), for small values of ϕ we have f ≈ 0, g ≈ 1 so that the shear stress modulus becomes τ = 1 − e −αβ .When Re is increased, the shear stress τ is drastically reduced in the central part of the fluid layer (where the center of rotation is the z−axis), while it increases within the boundary layers.Notice that, since the stress has a limiting threshold, the maximum of τ (which is attained at the plates) is always bounded by 1.In Figs. 5, 6 we show the dependence of the center of rotation upon the parameter α.The Reynolds number is taken equal to 1 to show that we can have the formation of boundary layers even for small Re.Notice that in this case ϕ is fixed.Looking at Figs. 5,6 we see that, for small α, the center of rotation is the z-axis in the central part of the fluid layer and moves to (0, ±1) in the boundary layers (whose thickness decreases with the decrease of α).This behavior is due to the particular constitutive equation we are considering.Indeed, looking at the non-dimensional form of apparent viscosity (16), we see that μ is bounded by α, i.e., μ α for all θ 0. So, a reduction of α means a reduction of the apparent viscosity everywhere in the fluid.Moreover, recalling that the shear stress modulus is given by τ = 1 − e −α γ , we have that τ ≈ 0 for α γ 1.This means that, away from the boundary layers, the fluid behaves as a rigid body that rotates around the z−axis, as shown in Figs. 5, 6. Finally in Fig. 7 we show the dependence of τ on the parameter α.As expected, the stress τ increases when α is increased and τ is always bounded by 1.In Figs. 8, 9 we show the dependence of f , g on the parameter β.The Reynolds number Re is always set to 1.The existence of boundary layers is evident when β is small, that is for a large aspect ratio ξ .When β = 5, i.e., ξ = 0.1 and ϕ = 0.02, we see that f ≈ 0 and g ≈ z.Indeed, in this case the solution is leading-order solution (24) that we have derived in Sect.3. Finally, in Fig. 10 we show the stress τ for various values of β.When β is small, the shear stress τ is zero in the central part of the layer and it increases in the boundary layers.When β is increased, the stress increases and for β > 1 the stress becomes almost uniform in the layer.As in the previous cases, τ is always bounded.We conclude by showing the relative error between the numerical solution and the asymptotic expansion obtained in (27).In Fig. 11 we show the relative error defined as follows where f a and g a are the approximated solutions given in ( 27) and f n and g n are the numerical solutions.We observe that the error for the component g is smaller because the approximation is of the second order, while the approximation for f is of the first order.

Conclusions
The flow in an orthogonal rheometer is a well-studied problem for a large variety of non-Newtonian fluids.
In this paper we have investigated the flow for a class of fluids recently investigated by Garimella et al. [15], for which the constitutive law is a nonlinear monotonically increasing bounded function.Assuming a specific form of the velocity field, the mathematical problem is reduced to a nonlinear BVP for the (x, y) coordinates of the center of rotation.The problem is discretized by a spectral collocation method and solved numerically via the iterative Newton-Rapson scheme.Boundary layers in the Navier-Stokes fluids are regions of concentration of vorticity and in the case of the Navier-Stokes fluids exhibit themselves at high Reynolds numbers.However, in non-Newtonian fluids, concentration of vorticity (boundary layers) can present themselves even when the Reynolds number associated with the flow is zero.That this is indeed the case was pointed out in [21], [22].Also, in non-Newtonian fluids, boundary layers with multiple structure, with the vorticity being concentrated in one layer, and elastic effects being concentrated in another layer are also possible as discussed in [23].Boundary layers are also possible in elasticity and thermoelasticity, wherein one has a concentration of strain gradients (see [24], [25], [26]).
Of course, the thickness of these boundary layers strongly depends on the material moduli characterizing the fluid and the appropriate non-dimensional numbers that stem from the material moduli.We have proved that, for small Reynolds numbers or for small aspect ratio, the solution can be sought in the form of an asymptotic expansion.We have determined the approximated solution up to the second order, and we have estimated the relative error between the approximated and the numerical solution.A possible development of the present work is certainly the extension to the case of different angular velocities.by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.Funding Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement.
7 Appendix A. The Jacobian of system (29) and the Newton-Rapson scheme For convenience, we report here the entries of the Jacobian of system (29).The Jacobian is given by where and where I is the identity matrix and rm[a] is the (N + 1) × (N + 1) matrix whose (identical) columns are formed by the vector a.The Jacobian J is thus a matrix formed by the (N + 1) × (N + 1) blocks defined above.The dimension of J is thus 2(N + 1) × 2(N + 1).Notice that μ(θ), • μ(θ ) are given by ( 16), (17), while The Newton-Rapson scheme at each step k consists in solving and get

Fig. 1
Fig. 1 Schematic of an orthogonal rheometer.The origin of the coordinate system is O

Fig. 11
Fig. 11 Relative errors (in percentage) between the numerical solutions and approximated solutions (27) as a function of ϕ