Squeeze Flow of Stress Power Law Fluids

: In this paper, we studied the squeeze ﬂow between circular disks of a new class of ﬂuids deﬁned by an implicit relation referred to as stress power law ﬂuids. The constitutive response of these ﬂuids was written expressing the symmetric part of the velocity gradient as a tensorial function of the Cauchy stress. We assumed that the aspect ratio between the gap separating the disks and the radius was small so that a lubrication expansion could be adopted. We wrote the general problem and looked for a solution that could be written in terms of the small aspect ratio parameter. We obtained a sequence of problems that could be solved iteratively at each order, and we focused on the leading and ﬁrst order, deriving explicit expressions for the velocity ﬁeld, stress, and pressure. of the expansion (calculations become more involved as the order of approximation is increased). For the zero- and ﬁrst-order problem, we determined the semi-explicit expression for the main variables, namely velocity, stress, pressure, and total force exerted by the ﬂuid on the squeezing plates. We plotted the contours of the solution for different choices of the physical parameters, investigating in particular the behavior at the reference time at the interface between the ﬂuid and the exterior. The graphical representations obtained can be used to compare the power stress law model with the Newtonian model that was obtained setting the ﬂow index n = 0.


Introduction
Recently, Malek et al. [1] studied incompressible fluids from a new perspective, namely that in which the constitutive equation is written expressing kinematical quantities as a function of the stress. This is in contrast with the general approach adopted by Noll [2] to describe the response of incompressible simple fluids, wherein the Cauchy stress is determined by the history of the deformation gradient. The new approach provides the appropriate way to express fluids such as the yield stress fluids or fluids with rheological parameters that depend on the pressure. The new class of fluids introduced in [1], that is termed stress power law fluids, includes the classical incompressible Navier-Stokes model as a special subclass. An important feature of stress power law models is that they automatically satisfy the incompressibility constraint without requiring the definition a Lagrange multiplier to enforce the constraint. In the stress power law model, the constitutive equation is defined by the symmetric part of the velocity gradient being a function of the deviatoric part of the Cauchy stress tensor.
This new class of fluids has many interesting applications, among which we must mention blood and colloidal dispersions; see [3][4][5][6]. Blood is indeed a very complex mixture of a liquid plasma, red blood cells, platelets, proteins, ions, and other components that in large vessels can be safely modeled as a single-component fluid. Many biochemical reactions that occur during blood flow result in significant mechanical effects that may lead, for instance, to a reduction of the shear rate for increasing shear stress (coagulation) or to an increase of the shear rate for increasing shear stress (disaggregation of red blood cells). This type of behavior may occur in different regions of the shear stress/shear rate curve, exhibiting thus a non-monotonic relation between the shear rate and the shear stress that can only be described by implicit constitutive laws.
The constitutive equation of a generalized Stokesian fluid is such that the Cauchy stress tensor is given as a tensorial function of the density ρ * and of the symmetric part of the velocity gradient (starred quantities denote dimensional variables). D * = 2 −1 (∇v * + ∇v * T ), namely: T * = F * (ρ * , D * ). (1) where D * is a function of T * , so that: In general, a constitutive equation of type (1) cannot be inverted to one of type (5) and the opposite is also true. In an even more general framework, one can assume that the constitutive relation is fully implicit, i.e., G * (T * , D * , ρ * ) = 0 In the paper [1] Malek et al. introduced a new class of fluid models wherein the tensor D * has a power law relationship with the tensor T * . They named this new class of fluids stress power law fluids. The constitutive equation for such fluids is of the type (5) and cannot be inverted in an explicit manner to obtain an expression such as (1). In particular, they proved that for special values of the constitutive parameters, the relation between the norm of the tensor D * and the norm of the tensor T * can be non-monotone. Simple flows of stress power law fluids were studied in [7][8][9][10][11][12].
Let us consider an incompressible stress power law fluid wherein the Cauchy tensor (here, p * is the Lagrange multiplier due to the incompressibility constraint (pressure); see [13] for a discussion on the notion of pressure) T * = −p * I + S * is such that S * is related to D * via: where α * > 0 has the dimension of the inverse of a viscosity and β * has the dimension of the inverse of a squared pressure. The flow index n is a real number, while γ 0 is non-dimensional. Notice that for (n, γ) = (0, 0) or for (β * , γ) = (0, 0), the classical Newtonian model with viscosity (2α * ) −1 is recovered. Taking the Frobenius norm on both sides of (7), we obtained: where |D * | 2 = D * · D * = trD * 2 . Relation (8) can be plotted in the first half-quarter of the R 2 -plane as a function |D * | = f (|S * |); see Figure 1. Following [11], we know that for n −1/2 or n < −1/2 and γ ≥ d n , with: the function |D * | = f (|S * |) defined in (8) is monotonically increasing and thus invertible as |S * | = f −1 (|D * |). On the other hand, if n < −1/2 and γ < d n , the relation (8) is no longer monotone and cannot be rewritten in the form |S * | = f −1 (|D * |). In Figure 1, we plot the constitutive relation (8) with α * = 1, β * = 2, and γ = 0.1 for different values of the index n.
In this paper, we studied the squeeze flow between parallel circular disks of a fluid whose constitutive equation is of the type of (7) under the assumption that the gap between the squeezing disks is much smaller than their radius (lubrication flow; see Figure 2). Squeeze flows are found in many engineering, biological, and rheological applications. In the context of biomechanical valves and diarthrodial joints are examples of squeeze flows relevant in biology and bioengineering. Other examples can be found in the context of food intake (chewing) where a material is squeezed between rigid surfaces (teeth) and where the gap between the rigid plates is smaller than the characteristic length of the squeezing surfaces. * * ℎ * ( * ) For a general review on squeeze flow between disks of both Newtonian and non-Newtonian fluids, we refer the reader to [14][15][16][17]. We formulated the problem using cylindrical coordinates and assuming azimuthal symmetry so that all the variables do not depend on the angular coordinate. We looked for a solution that can be expressed as power series of the small aspect ratio parameter ε that represents the ratio between the half-distance that separates the disks and their radius. We focused on the leadingand first-order approximations for which we determined semi-analytical solutions for the velocity, stress, and pressure fields. We shall finally provide a graphical representation of the velocity, stress, and pressure fields, and we discuss the obtained results.

Squeeze Flow
We considered a portion of fluid squeezed between parallel disks of radius R * placed at a distance 2h * (t * ); see Figure 2. We assumed that the plates squeezed the flow so that dh * /dt * < 0, and we also assumed azimuthal symmetry so that there was no dependence on the angular coordinate θ. We defined H * = h * (0) as the initial position of the upper plate and assumed that: The velocity field is given by: The constitutive equation is T * = −p * I + S * where S * is related to the symmetric part of the velocity gradient D * via the implicit relation (7). The balance of mass and linear momentum in cylindrical coordinates yields: Exploiting the symmetry with respect to z, we write the boundary conditions as: where Γ * represents the interfacial tension and R * is the sample radius, which is equivalent to the radius of curvature at the surface, assuming that the sample maintains a cylindrical shape. The quantity p * out represents the external pressure at r * = R * . On the upper plate, the coefficient λ ∈ [0, 1] discriminates among different boundary conditions. Indeed, we observed that λ = 0 corresponds to the no-slip condition, λ = 1 represents free-slip, and λ ∈ (0, 1) stands for partial-slip with a friction coefficient that is related to the coefficient K * . We rescaled the system with the following non-dimensional variables: so that (7) becomes: while (11)- (13): ∂v z ∂z The non-dimensional boundary conditions acquire the form: whereḣ(t) = dh/dt and where: is the capillary number representing the relative effect of viscous drag forces versus surface tension forces acting across the interface between the liquid and the atmosphere at r = 1.
In the following, we assumed that Ca = O(1). We looked for a solution that could be expressed as a power series of ε, namely: Using the expressions above, we could reformulate the mathematical problem for each order of approximation (leading, first, etc.). The system was therefore transformed into a series of simpler problems that could be solved iteratively. We focused on the leading-and first-order of approximation only, but it is easy to see that our procedure could be extended to any grade of approximation. We observed that from our scaling, the symmetric part of the velocity gradient is given by: implying that S rθ = S θz = 0. It is straightforward to prove that: From (25) and (26), we immediately see that S zz = 0, meaning that the normal stresses are negligible at the leading order and that the only non-zero component of the deviatoric stress is S (0) rz . Exploiting this fact and expansion (24), we found that the term in the square brackets in (26) could be rewritten expanding in Taylor series around ε: The components of the tensor D up to the first order thus become: Let us now see how to write the boundary conditions on z = h(t). To be as general as possible, we assumed that h depends on ε so that expanding h = h(t; ε) around ε = 0, we write: and : Expanding the right-hand side of the above in ε, we found: Hence, the first boundary condition in (21) 1 becomes: Using the same argument, we found that the second boundary condition (21) 1 can be rewritten as: λKS

Remark 1.
When h(t) is independent of ε, we obtained h(t) ≡ h (0) (t), and the boundary conditions on z = h (0) reduced to: The dimensional force per unit surface (stress) exerted by the fluid on the upper plate is given by: The normal force per unit surface acting on the plate is thus: Rescaling f * as the pressure, we found.
The resultant non-dimensional normal force is thus: We shall see that p (0) and p (1) depend only on r and t so that, making use of the asymptotic expansion and recalling that S (0) zz = 0, we found that:

Zero-Order Approximation
At the leading order, the problem is reduced to: The boundary conditions are: We immediately realized that p (0) = p (0) (r, t) and S rz so that from (31): The integration of (35) provides the radial component of the velocity field.
3.1. The Case λ = 1 We noticed that when λ = 1, the pressure gradient vanished, G (0) ≡ 0. This was due to the fact that perfect slip implies no change in the internal pressure as the fluid is squeezed.
Hence, p (0) = 0 because of the boundary condition (34) 4 where we used Condition (34) 2 . The imposition of Condition (34) 1 yields: where we exploited Condition (34) 3 . Finally, replacing v Notice that in this case, the solution did not depend either on n or β and, hence, on the particular constitutive equation we considered.
3.2. The Case λ ∈ [0, 1) Let us now consider the case λ ∈ [0, 1). Recalling that S (0) rz = G (0) (r, t)z, the boundary condition for v r on the upper plate is: The integration of (35) provides: or equivalently: for n = −1 and: for n = −1. We observed that: where: when n = −1, and: when n = −1. Equation (42) is useful for the evaluation of the solution at the first order.
Integrating the continuity equation between z and h (0) , we found: After some calculations, we obtained: Equations (40), (41), and (46) provide the expressions of the velocity components at the leading order. To determine the pressure gradient G (0) , we imposed the symmetry conditions (34) 2 so that: Integrating the above between zero and r and exploiting the symmetry condition G (0) (0, t) = 0, we obtained: Equation (48) is a nonlinear equation in G (0) whose solution provides the pressure gradient G (0) (r, t).

Remark 3. We observed that: lim
irrespective of n and G (0) , so that: Integrating between r and 1 and imposing p (0) = 0 on r = 1, we found: Equation (51) shows that our model made sense only if (dh (0) /dt)/h (0) 2 remained bounded in the time interval considered. If after some time, this was no longer valid, then the asymptotic expansion lost its validity since the pressure field within the fluid became too large. It is easy also to check that in the Newtonian case (β → 0 and γ = 0): The general case (47) does not provide an explicit expression for G (0) , and the pressure gradient must be determined numerically.

First-Order Approximation
At the first order, the problem is: implying that the pressure is again independent of z, i.e., p (1) = p (1) (r, t). The boundary conditions are: From the definition of the stress components (28)-(30), we found: showing that the normal stresses do not vanish at the first order. Let us now introduce: which is a function that is known from the leading-order problem. Notice that f (0) requires the knowledge of the time derivative of v (0) r , whose explicit form is given in (42). Integrating the momentum equation with the boundary condition S (1) rz = 0 on z = 0, we found: where G (1) = ∂p (1) /∂r. Again, we distinguished between the free-slip condition and partial-/no-slip condition on z = h (0) .

The Case λ = 1
When λ = 1, it is easy to check that S (1) rz = 0 on z = h (0) and that: As a consequence, S rz is linear in z, and since S (1) rz = 0 on z = h (0) and z = 0, we see that S (1) rz ≡ 0. From (59), we obtained: since g (0) (r, z, y) = 2(1 + γ) (recall that S (0) rz ≡ 0 when λ = 1). The pressure gradient at the first order does not vanish (as it does at order zero). Equation (60) provides the pressure gradient in the fluid in the case of perfect slip at the first order. Setting S (1) r (r, t), so that, integrating the mass balance at the first order between zero and z and recalling that v Now, we imposed the boundary condition (53) 1 that, recalling (38), can be rewritten as: Evaluating (61) on z = h (0) and integrating in r, we obtained: The solution at the first order is therefore completely determined by the solution at the leading order. Notice that, as it occurs at the leading order, the radial velocity is linear in r and does not depend on z, whereas the vertical velocity is linear in z and does not depend on r. Recalling (37), (38), and (54)-(57), we found: which shows that, when λ = 1, the normal stresses depend only on time at the first order and that the shear stress is identically null.

Remark 4.
From the definition of g (0) and from the boundary condition S (0) rz = 0 on r = 0 and on z = 0, we see that: Therefore, a necessary condition for the denominator of (69) to vanish is that g (0) (r, z, t) = 0 for some r ∈ (0, 1], t 0 and z ∈ (0, h (0) (t)]. We proved here that when the constitutive relation (8) is monotone, i.e., n −1/2 or n < −1/2 and γ ≥ d n , then the denominator of (69) is always strictly positive. Indeed, let us suppose that g (0) (r, z, t) = 0 for some (r, z) ∈ (0, 1] × (0, h (0) (t)] and t 0. Then: Recalling that γ is a non-negative parameter, the above equation has a solution only if the term in the square brackets of (70) is non-positive, i.e., only if: Hence, a first necessary condition for g (0) to vanish is that n < −1/2 (if n 1/2 the function g (0) (r, z, t) is always strictly positive). Now, let us suppose that n < −1/2 and set X = βS (0) 2 rz 0. Equation (70) becomes: It is straightforward to see that m(0) = −1 and m(∞) = 0 (recall that we assumed that n < −1/2). Differentiating the function m with respect to X, we obtained: m (X) = −n(1 + X) n−2 (2n + 1)X + 3 , so that: (the solution X M = −1 is clearly disregarded). Now, we notice that: Therefore, Equation (71) holds true only if n < −1/2 and γ < d n , i.e., when the constitutive relation is non-monotone. In conclusion, we can state that, whenever the constitutive relation is monotonically increasing, the denominator in (69) is always strictly positive and the pressure gradient at the first order does not blow up.

Results and Discussion
In this section, we show the contours/plots of the main variables involved in the problem for a given squeeze function h(t). We did not have in mind any practical application; our intent was to give the reader an idea of the behavior of the system for different choices of the model parameters. We also observed that it was possible to extend our analysis to higher orders (second, third, etc.) to get a finer approximation of the solution. The procedure remained the same with the calculations becoming more involved.
To investigate the dependence on the physical parameters, we considered the following squeeze function: so that:ḣ The plot of the function h(t; ε) together with its approximation is shown in Figure 3. Notice that h(t; ε) > 0 for all time t 0, so that the plates cannot come in touch in a finite time.
The values of the parameters appearing in the constitutive equation and in the boundary conditions are the ones of Table 1.
In Figure 4, we plot the contours of the velocity components, the contours of the tangential stress, and the flow pattern at the leading order. The parameters used in the simulations were the ones of Table 1. Time was set to t = 0.5 so that h (0) = 0.8944. As physically expected, the modulus of the shear stress was maximum in the vicinity of the upper right corner of the (r, z) domain, while it vanished on the axes r = 0 and z = 0. The radial velocity v (0) r decreased with z and increased with r, whereas the vertical velocity v (0) z was almost uniform on z = const. In Figure 5, we show the velocities v z and the stress S (1) rz still using the parameters of Table 1. We noticed that the radial velocity v (1) r was negative (directed towards the center r = 0) except in the proximity of the axis r = 0 and near the plate. As a consequence, the radial velocity at the zero order was slightly larger than that obtained at the first order (the correction at the first order was negative almost everywhere). Recall that the first-order radial velocity is given by v   Table 1.  Table 1.
The stress S (1) rz was positive almost everywhere and reached its maximum at the external surface r = 1, showing that the modulus of the shear stress at the first order was slightly smaller than the one at the zero order (recall that the shear stress at the leading order was negative everywhere). The normal stresses at the zero order were null, so we plot only S (1) θθ , and S (1) zz ; see Figure 6. It is interesting to notice that the components S (1) rr and S (1) θθ were identical. This was observed also for other choices of the parameters. Notice also that the moduli of the normal stresses were maximum on z = 0. zz . The contours were evaluated at time t = 0.5 with h (0) (0.5) = 0.8944 and using the coefficients given in Table 1.
In Figure 7, we plot the shear stress, velocity components, and normal stresses at r = 1 and t = 1 (so that h (0) (1) = 0.7071) for the indices j = (0) , (1) and for n = −5, −2, 0, 5. The parameters here were λ = 0 (no-slip), γ = 1, β = 3, and Re = 3. These plots allowed us to study the dependence of the main variables on the index n. Looking at the shear stress plot, we noticed that S (0) rz and S (1) rz exhibited opposite monotonicity (with respect to z) irrespective of n. In particular, we noticed that at the leading order, the increase of n resulted in a decrease in the modulus of the shear stress, so that the shear between adjacent layers was reduced. This was probably due to the fact that the fluid shear thinned as n was increased. The velocity components v zz increasing with z. We then illustrated the dependence of the pressure on n at time t = 1. In Figure 8, we show the pressures p (0) and p (1) and the pressure at the first order p (0) + εp (1) at time t = 1. The parameters are the same as Figure 7. We noticed that whereas p (0) decreased with n the pressure, p (1) actually increased with n except for n = 5 (where the variation with r was very small). The pressure at the leading order (green) decreased with n, showing that the correction p (1) did not vary the monotonicity of the pressure at the leading order. This behavior was again in accordance with the fact that shear-thinning effects were enhanced as n increased. Notice finally that in all the plots, the Newtonian-like behavior was the one in which n = 0.
zz . The parameters are λ = 0 (no-slip), γ = 1, β = 3, and Re = 3: "o and |" (n = −5); "+ and " (n = −2); " * and " (n = 0), " and x" (n = 5). Finally, in Figure 9, we show the total force exerted by the fluid on the plate; see Remark 2. We noticed that F (0) increased with time, while F (1) exhibited the opposite behavior. The total force at the first order F slightly increased at the beginning and then decreased with time. The correction at the first order in this case seemed to be quite relevant. Figure 9. Total forces F (0) (t), F (1) (t), and F (t) exerted by the fluid on the plate as a function of time for different n with λ = 0, β = 3, Re = 3, p out = 0.2: 'o and |' (n = −5); + and (n = −2); * and (n = 0), and x (n = 5). On the right, the total force at the first order, ε = 0.3.
The numerical results indicated that, within the context of squeeze flow, the deviation from the Newtonian behavior seemed to have a stronger effect on the stress distribution and on the pressure than on the velocity. This was probably due to the shear-thinning or shear-thickening nature of the stress power law model considered. In particular, we noticed that the pressure at the leading order could change by an order of magnitude with the flow index n. This is an important result if we think in terms of applications (think for instance at the synovial fluid between hyaline cartilage that can be modeled using lubrication approximation), since it demonstrated that the constitutive response of the fluid plays an important role in term of forces and stress acting within the fluid and on the squeezing surfaces.

Conclusions
We presented a mathematical model for the axisymmetric squeeze flow of a fluid whose constitutive equation is of an implicit type. We assumed that the geometrical setting was such that lubrication approximation could be used. We looked for a solution that was an asymptotic expansion in terms of the small aspect ratio parameter. We focused only on the zero-and first-order approximation, but the procedure adopted can be used to derive also higher order terms of the expansion (calculations become more involved as the order of approximation is increased). For the zero-and first-order problem, we determined the semi-explicit expression for the main variables, namely velocity, stress, pressure, and total force exerted by the fluid on the squeezing plates. We plotted the contours of the solution for different choices of the physical parameters, investigating in particular the behavior at the reference time at the interface between the fluid and the exterior. The graphical representations obtained can be used to compare the power stress law model with the Newtonian model that was obtained setting the flow index n = 0.
Author Contributions: Writing-original draft preparation, L.F. and A.B.; writing-review and editing, L.F. and A.B. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.