Squeeze flow of a Hershel--Bulkley fluid

TL;DR

This paper presents an asymptotic solution for the axisymmetric squeeze flow of a Hershel--Bulkley viscoplastic fluid, extending previous asymptotic techniques to this specific non-Newtonian fluid behavior.

Contribution

It introduces a novel asymptotic approach tailored for the axisymmetric squeeze flow of Hershel--Bulkley fluids, building on earlier methods for different fluid types.

Findings

Derived an asymptotic solution for Hershel--Bulkley squeeze flow.

Extended existing asymptotic techniques to viscoplastic fluids.

Provides analytical insights into complex non-Newtonian flow behavior.

Abstract

We develop an asymptotic solution for the axisymmetric squeeze flow of a viscoplastic Hershel--Bulkley medium following the asymptotic technique suggested earlier by Balmforth and Craster (1999) and Frigaard and Ryan (2004).

Equations

$$\varepsilon Re\Bigl{(}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}+w\frac{\partial u}{\partial z}\Bigr{)}=-\frac{\partial p}{\partial r}+\varepsilon^{2}\frac{\partial\tau_{rr}}{\partial r}+\frac{\partial\tau_{rz}}{\partial z}+\varepsilon^{2}\frac{\tau_{rr}-\tau_{\theta\theta}}{r},$$

$$\varepsilon^{2}Re\Bigl{(}\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial r}+w\frac{\partial w}{\partial z}\Bigr{)}=-\frac{\partial p}{\partial z}+\varepsilon^{2}\frac{\partial\tau_{rz}}{\partial r}-\varepsilon^{2}\frac{\partial\tau_{rr}}{\partial z}-\varepsilon^{2}\frac{\partial\tau_{\theta\theta}}{\partial z}+\varepsilon^{2}\frac{\tau_{rz}}{r},$$

$$\frac{\partial u}{\partial r}+\frac{u}{r}+\frac{\partial w}{\partial z}=0.$$

$${\sigma }_{ij}=-p{\delta }_{ij}+{\tau }_{ij},\{\begin{array}{ll}{\tau }_{ij}=({\dot{\gamma }}^{n-1}+\frac{B}{\dot{\gamma }}){\dot{\gamma }}_{ij},& \text{mbox}if\tau >B,\\ {\dot{\gamma }}_{ij}=0,& \text{mbox}if\tau \le B.\end{array}$$

$${\dot{\gamma}}_{rr}=2\frac{\partial u}{\partial r},\quad{\dot{\gamma}}_{rz}=\Bigl{(}\frac{\partial u}{\partial z}+\varepsilon^{2}\frac{\partial w}{\partial r}\Bigr{)},\quad{\dot{\gamma}}_{zz}=2\frac{\partial w}{\partial z},\quad{\dot{\gamma}}_{\theta\theta}=2\frac{u}{r},$$

$$\epsilon =\frac{\hat{H}}{\hat{R}},Re=\frac{\hat{\rho }\hat{W}\hat{R}}{\hat{\mu }},B=\frac{\widehat{{\tau }_0}\hat{H}}{\hat{\mu }\hat{U}}.$$

$$-\frac{\partial p}{\partial r}+{\epsilon }^2\frac{\partial {\tau }_{rr}}{\partial r}+\frac{\partial {\tau }_{rz}}{\partial z}+{\epsilon }^2\frac{{\tau }_{rr}-{\tau }_{\theta \theta }}{r}=0,$$

$$-\frac{\partial p}{\partial z}+{\epsilon }^2\frac{\partial {\tau }_{rz}}{\partial r}-{\epsilon }^2\frac{\partial {\tau }_{rr}}{\partial z}-{\epsilon }^2\frac{\partial {\tau }_{\theta \theta }}{\partial z}+{\epsilon }^2\frac{{\tau }_{rz}}{r}=0,$$

$$\frac{\partial u}{\partial r}+\frac{u}{r}+\frac{\partial w}{\partial z}=0.$$

$$\begin{array}{rl}& u^s=u^{s,0}+\epsilon u^{s,1}+{\epsilon }^{s,2}{u}^2\dots ,w^s=w^{s,0}+\epsilon w^{s,1}+{\epsilon }^2{w}^{s,2}\dots ,\\ & p^s=p^{s,0}+\epsilon p^{s,1}+{\epsilon }^2{p}^{s,2}\dots ,{\tau }_{ij}^s={\tau }_{ij}^{s,0}+\epsilon {\tau }_{ij}^{s,1}+{\epsilon }^2{\tau }_{ij}^{s,2}\dots .\end{array}$$

$$\mathcal{O}(1)$$

$$\mathcal{O}(\epsilon )$$

$${\int }_0^1{u}^0(r,z)\text{leavevmode}dz=\frac{r}{2},{\int }_0^1{u}^1(r,z)\text{leavevmode}dz=0.$$

$$p^{s,0}=p_0(r),{\tau }_{rz}^{s,0}=z{p}_0^{\prime }(r),$$

$$\tau_{rz}^{s,0}=-\Bigl{(}\Bigl{|}\frac{\partial u^{s,0}}{\partial z}\Bigr{|}^{n}+B\Bigr{)},\quad\tau_{rz}^{s,1}=n\Bigl{|}\frac{\partial u^{s,0}}{\partial z}\Bigr{|}\frac{\partial u^{s,1}}{\partial z}.$$

$$u^{0}(r,z)=\begin{cases}\Bigl{(}\frac{B}{z_{0}}\Bigr{)}^{1/n}\frac{1}{1+1/n}[(1-z_{0})^{1+1/n}-(z-z_{0})^{1+1/n}],&z\in[0,\,z_{0}],\\ \Bigl{(}\frac{B}{z_{0}}\Bigr{)}^{1/n}\frac{1}{1+1/n}(1-z_{0})^{1+1/n},&z\in(z_{0},\,1].\end{cases}$$

$$u_0(r)=\frac{B^{1/n}(1-z_0{)}^{1+1/n}}{z_0^{1/n}(1+1/n)}.$$

$$p^{s,0}=p_0(r),p_0^{\prime }=-\frac{B}{z_0(r)},{\tau }_{rz}^{s,0}=-\frac{Bz}{z_0(r)}.$$

$$\frac{(1-z_{0})^{2+1/n}}{(2+1/n)}-(1-z_{0})^{1+1/n}+\Bigl{(}\frac{z_{0}}{B}\Bigr{)}^{1/n}\frac{r}{2}(1+1/n)=0.$$

$$u^p=u^{p,0}(r)+\epsilon u^{p,1}(r,z)+{\epsilon }^2{u}^{p,2}(r,z)+\dots ,$$

$${\tau }_{rr}^{p,-1}=\frac{2B}{{\dot{\gamma }}^{p,0}}\frac{\partial u^{p,0}}{\partial r},{\tau }_{rz}^{p,0}=\frac{B}{{\dot{\gamma }}^{p,0}}\frac{\partial u^{p,1}}{\partial z},{\tau }_{\theta \theta }^{p,-1}=\frac{2B}{{\dot{\gamma }}^{p,0}}\frac{u^{p,0}}{r},$$

$$\mbox{where}\quad\dot{\gamma}^{p,0}=\sqrt{\Bigr{(}\frac{\partial u^{p,1}}{\partial z}\Bigl{)}^{2}+4\Bigl{[}\Bigl{(}\frac{\partial u^{p,0}}{\partial r}\Bigl{)}^{2}+\Bigl{(}\frac{u^{p,0}}{r}\Bigr{)}^{2}+\frac{\partial u^{p,0}}{\partial r}\frac{u^{p,0}}{r}\Bigr{]}}.$$

$$\mathcal{O}(1)$$

$$\mathcal{O}(\epsilon )$$

$$p^{p,0}(r)=p_0(r),p_0^{\prime }=-\frac{B}{z_0},{\tau }_{rz}^{p,0}=-\frac{Bz}{z_0},u^{p,0}=u_0.$$

$$z_0^{\prime }=-\frac{(n+1)(2n+1)z_0^{1+1/n}}{2B^{1/n}(1-z_0{)}^{1/n}(2n^2{z}_0^2+2n{z}_0+1+n)}<0.$$

$$p^{s,1}=p_1(r),{\tau }_{rz}^{s,1}=z{p}_1^{\prime }(r)+g(r),$$

$$u^{s,1}=\frac{p_{1}^{\prime}(r)}{n+1}\Bigl{[}(z-z_{0})^{\frac{1}{n}}(z+nz_{0})-(1-z_{0})^{\frac{1}{n}}(1+nz_{0})\Bigr{]}+g(r)\Bigl{[}(z-z_{0})^{\frac{1}{n}}-(1-z_{0})^{\frac{1}{n}}\Bigr{]},$$

$$({\tau }^{p,-1}{)}^2=({\tau }_{rr}^{p,-1}{)}^2+({\tau }_{\theta \theta }^{p,-1}{)}^2+{\tau }_{rr}^{p,-1}{\tau }_{\theta \theta }^{p,-1}+({\tau }_{rz}^{p,0}{)}^2$$

$$\displaystyle=\frac{4B^{2}}{(\dot{\gamma}^{p,0})^{2}}\Bigl{[}(u_{0}^{\prime})^{2}+\bigl{(}\frac{u_{0}}{r}\bigr{)}^{2}+u_{0}^{\prime}(r)\frac{u_{0}}{r}\Bigr{]}+\frac{B^{2}z^{2}}{z_{0}^{2}}=B^{2}.$$

$$\dot{\gamma}^{p,0}=\frac{\eta z_{0}}{\sqrt{z_{0}^{2}-z^{2}}},\quad\mbox{where}\quad\eta=2\sqrt{\bigl{(}u_{0}^{\prime}\bigr{)}^{2}+\bigl{(}\frac{u_{0}}{r}\bigr{)}^{2}+u_{0}^{\prime}\frac{u_{0}}{r}}.$$

Full text

1 Introduction

Viscoplastic or yield-stress fluids are materials which behave like a solid below a critical yield stress and flow like a viscous fluid for stresses higher than this threshold. The flow field is thus divided into unyielded (rigid) and yielded (fluid) zones. The surface separating a rigid from a fluid zone is known as a yield surface. The location and shape of the latter must be determined as part of the solution of the flow problem.

The squeeze flow of viscoplastic materials has been investigated many times, see the detailed reviews [1]– [5] for a comprehensive list of references to available analytical, numerical and experimental results [6]–[19] and [20]–[15] for some more recent contributions. Different constitutive equations have been used, in both theoretical and numerical studies, such as the original Bingham model by Covey and Stanmor [6], Lipscomb and Denn[7], Sherwood and Durban [8], Smyrnaios and Tsamopoulos [9], and Roussel et al. [10], the bi-viscosity model by O’Donovan and Tanner [11]] and Wilson [12], the Hershel-Bulkley model by Covey and Stanmor [6], Sherwood and Durban [13], an elasto-viscoplastic model by Adams et al. [14], Bingham fluid with a deformable core by Fusi et al. [15] and the regularized Papanastasiou model by Smyrnaios and Tsamopoulos [9], Matsoukas and Mitsoulis [16], [17], and Karapetsas and Tsamopoulos [18].

We would like to construct the consistent asymptotic solution for the squeeze flow of Hershel-Bulkley material. Walton and Bittleston [27] studied the axial flow in an eccentric annular duct and showed analytically that a true plug exists in the middle of the channel on both the wide and narrow part of the annulus, with a pseudo-plug placed between the two rigid zones. Balmforth and Craster [28], Frigaard and Ryan [29], [30] suggested the asymptotic technique that allows constructing the consistent solution for thin–layer problems. The asymptotic solution was developed for a fluid flowing down an inclined plane [28] and for the flow along a channel of slowly varying width [29], [30]. Recently, Muravleva [21], [22] has analysed the planar and axisymmetric squeeze flows of a Bingham fluid and the axisymmetric squeeze flows of a Casson fluid [22] exploiting the asymptotic technique introduced in [28]–[30]. In this article, we continue our research into the squeeze flow problem [21]-[22] and consider the axisymmetric squeeze flow of a Hershel-Bulkley material, following an approach developed in [28]–[30].

The paper is organized as follows. In Section 2 the dimensionless governing equations of the flow are presented. Section 3 is the Finally, in Section 5, a summary of the results is given.

2 Problem statement.

Squeeze Flow (SF) is the process in which a fluid is squeezed between two approaching parallel plates resulting in a radial flow, outward from the center. The geometry of the problem is shown in Fig. 1: we use an axisymmetric cylindrical polar coordinate system $(r;\theta;z)$ to describe the squeezing of a cylinder of incompressible Hershel-Bulkley fluid with radius $\hat{R}$ and height $2\hat{H}$. The fluid has density $\hat{\rho}$, yield stress $\hat{\tau_{0}}$ and plastic viscosity $\hat{\mu}$. The plates are squeezed together at a velocity $\hat{W}$, inducing a flow. We denote the dimensional variables with a hat symbol. We have scaled lengths in the $r$ and $z$ directions differently, with the disk radius $\hat{R}$ and with the half-distance $\hat{H}$ between the disks, respectively. $\hat{W}$ is taken as the characteristic velocity in the transverse direction, and the radial velocity component is scaled with $\hat{U}=\hat{W}\hat{R}/\hat{H}$. The pressure is scaled with $\hat{\mu}\hat{W}^{n}\hat{R}^{n+1}/\hat{H}^{2n+1}$, and time with $\hat{H}/\hat{W}$. The shear-stress components are scaled with $\hat{\mu}\hat{W}^{n}\hat{R}^{n}/\hat{H}^{2n}$, the extensional stresses with $\hat{\mu}\hat{W}^{n}\hat{R}^{n-1}/\hat{H}^{2n-1}$.

The flow is governed by the dimensionless conservation equations of momentum and mass:

[TABLE]

[TABLE]

[TABLE]

The Herschel-Bulkley plastic stresses are related to the strain rates through the constitutive equations

[TABLE]

where {\tau}=\sqrt{\frac{1}{2}\mbox{\boldmath\tau}:\mbox{\boldmath\tau}} and \dot{\gamma}=\sqrt{\frac{1}{2}\dot{\mbox{\boldmath\gamma}}:\dot{\mbox{\boldmath\gamma}}} denote the second invariants of $\tau$ and {\dot{\mbox{\boldmath\gamma}}}, and the strain rate tensor {\dot{\mbox{\boldmath\gamma}}} is given by

[TABLE]

The scalings introduce the small aspect ratio $\varepsilon$, and the Reynolds and Bingham numbers:

[TABLE]

We consider $Re\ll 1$. Neglecting the inertial terms, (1)–(3) are replaced by:

[TABLE]

[TABLE]

[TABLE]

Exploiting the symmetry of the axisymmetric flow about the plane $z=0$, we solve these equations over the domain $0\leq z\leq 1$, $0\leq r\leq 1$ subject to no-slip conditions, $u=u_{b}$, $w=-1$ on the surface of the disc $z=1$, symmetry conditions $\tau_{rz}=0$, $w=0$ on the plane symmetry $z=0$ and $u=0$, $\tau_{rz}=0$ on the axis symmetry $r=0$, and stress-free $\sigma_{rr}=-p+\varepsilon^{2}\tau_{rr}=0$, $\tau_{rz}=0$ at the edge $r=1$.

3 Asymptotic expansions.

We now solve the equations by introducing an asymptotic expansion. First, we consider shear flow near the plate for which we may find a solution through a straightforward expansion of the equations. This shear solution is denoted with a superscript (s).

We consider regular expansions in $\varepsilon$ of form:

[TABLE]

We substitute these expansions into the governing equations (7)–(9), and collect together the terms of the same order.

The lubrication equations for the first two orders are:

[TABLE]

[TABLE]

From the conservation of mass we have $\int_{0}^{1}u(r,z)\leavevmode\nobreak\ dz=\frac{r}{2}$. We require that

[TABLE]

3.1 Zero-order approximation

After the solution of the equations (13), (13) we have

[TABLE]

where $p_{0}$ is a function only of $r$ and the prime $()^{\prime}$ represents derivative ${\frac{d}{dr}}$. The leading order second invariants of strain rate and stress are given by $\tau^{s,0}=|\tau^{s,0}_{rz}|$, \dot{\gamma}^{s,0}=\Bigl{|}{\frac{\partial u^{s,0}}{\partial z}}\Bigr{|}. We are looking for a solution with $\frac{\partial u^{s,0}}{\partial z}<0$ in the domain $r>0$, $z>0$. Provided the yield stress is exceeded, the stress tensor components become

[TABLE]

We note that the leading order second invariant of stress $\tau^{s,0}=|\tau_{rz}^{0}|=z|p_{0}^{\prime}(r)|$, so $\tau^{s,0}$ attains its maximum at $z=1$ and vanish at $z=0$. Therefore there exists the point $z=z_{0}$ at which $\tau^{s,0}=B$ and $\dot{\gamma}^{0}=0$. Hence, at leading order, the yield condition holds at this point and $z_{0}(r)=\frac{B}{|p_{0}^{\prime}(r)|}$ is the position of the pseudo-yield surface. For $z\in[0,\,z_{0}]$ we have $\tau^{s,0}<B$ and $\dot{\gamma}^{s,0}=\frac{\partial u^{s,0}}{\partial z}=0$. The velocity field is now given by

[TABLE]

We denote the pseudo-plug velocity by $u_{0}(r)$:

[TABLE]

The expressions (18) can be written as

[TABLE]

Substituting (20) into the first equation of the flow rate constraint (17) leads to the following equation where the unknown is the pseudo-yield surface $z_{0}(r)$

[TABLE]

Evidently, $z_{0}(r)$ and therefore $u_{0}(r)$ depend on position $r$, highlighting how the flow in $0<z\leq z_{0}$ is only a pseudo-plug, which is an extension in the radial direction. Thus, the pseudo-yield surface $z_{0}(r)$ separates the entire area occupied by the material into subregions: fully yielded zones located near the plates (shear region), and pseudo-plug containing the central plane $z=0$. In the pseudo-plug, the asymptotic expansion (10) breaks down. Therefore, we look for a slightly different asymptotic expansion of the radial velocity component, where the property $\frac{\partial u^{p,0}}{\partial z}=0$ at $z_{0}$ is explicitly built in:

[TABLE]

The pseudo-plug solution is denoted with a superscript (p). The stresses are now given by

[TABLE]

[TABLE]

The lubrication equations for the first two orders are:

[TABLE]

Using symmetry about the center plane, we have $w^{p}\left.\right|_{z=0}=0$, $\tau_{rz}^{p}\left.\right|_{z=0}=0$, $\frac{\partial u^{p}}{\partial z}\left.\right|_{z=0}=0$. Integrating (29), (29) and enforcing continuity $p^{0}$, $\tau_{rz}^{0}$, $u^{0}$ at $z=z_{0}(r)$ leads to

[TABLE]

Fig. 2 shows the graphs of the pseudo-yield surfaces $z=z_{0}(r)$, for different values of $B$ and $n$.

We see that $z=z_{0}(r)$ is the decreasing function of radius. This property is easy to confirm analytically, differentiating $z_{0}$ (23) with respect to $r$ ($0<z_{0}\leq 1$)

[TABLE]

Graphs corresponding to large values of the Bingham number (for fixed values of power-law index $n$) are located above. For a fixed Bingham number and variable power-law index $n$, the pattern is more complicated: for smaller values of $n$, the graphs are located higher near the center of the plate, and quickly decrease with increasing $r$.

3.2 The first-order approximation.

In the shear region we integrate (16), (16), (19), using $\tau_{rz}(r,0)=0$, and receive:

[TABLE]

[TABLE]

where $p_{1}$ is a function only of $r$ and $g(r)$ is an unknown function of integration.

In the pseudo-plug region from (25) and (33) we obtain the second invariant $\tau^{p,-1}$ of the stress tensor, which is equal to $B$ since the pseudo-plug region is just at the point of yielding:

[TABLE]

After minor calculation we have

[TABLE]

Substituting (38) into (26), we obtain the equation for $\frac{\partial u^{p,1}}{\partial z}$. Solving this equation, integrating with respect to $z$ and matching the first order velocities $u^{s,1}$ (35) and $u^{p,1}$ at $z=z_{0}$, we get

[TABLE]

Inserting (33), (38) into (25) we find that

[TABLE]

Substituting (40) into (32), integrating the resulting equation and enforcing continuity of the pressure $p^{s,1}$ (35) and $p^{p,1}$ at the pseudo-yield surfaces gives

[TABLE]

We substitute (40), (41) into (32), integrate with respect to $z$, using $\tau_{rz}^{p,1}\left.\right|_{z=0}=0$, and find

[TABLE]

Enforcing continuity of the first order shear stress $\tau_{rz}^{s,1}$ (35) and $\tau_{rz}^{p,1}$ (3.2) at the pseudo-yield surface leads to

[TABLE]

To find $p_{1}$, we insert (36), (39) into the flow rate constraint (17), integrate and find

[TABLE]

Using (44), (34) we get the expression for $p_{1}^{\prime}$:

[TABLE]

Smyrnaios and Tsamopoulos [9] showed that unyielded material could only exist around the two stagnation points of flow. The numerical modeling for Bingham fluid [9], [16], [23], experiments [10] and asymptotic solution [22] confirmed the results [9]. It is interesting to investigate the obtained asymptotic solution for Herschel-Bulkley fluid near the stagnation point. We will examine the second invariant of the stress $\tau$ at the plate at $z=1$.

[TABLE]

where $p_{1}^{\prime}(r)$ and $g(r)$ are determined by (43), (45). The leading order $\tau^{0}=\frac{B}{z_{0}(r)}$ exceeds the yield stress, because of $z_{0}(r)\leq 1$, hence, according to the zero-order solution, the Herschel-Bulkley fluid in the shear region is yielded. To analyze the behavior of the function $\tau_{rz}(r,1)$ near the axis of symmetry, we expand function $z_{0}(r)$ in Taylor series. The expressions for $\tau_{rz}^{0}(r,1)$ (22) and $\tau_{rz}^{1}(r,1)$ (35), (43), (45) takes the form:

[TABLE]

[TABLE]

The expression in the brackets on the right-hand side of (47) tends to infinity when $r\rightarrow 0$. In the limit $r\rightarrow 0$, we have, \tau_{rz}(r,1)=-\frac{B}{z_{0}(r)}+\varepsilon\Bigl{(}p_{1}^{\prime}(r)+g(r)\Bigr{)}\rightarrow\infty. Therefore, about $r=0$ there is a point $r=r_{0}$, at which $\tau_{rz}=-B$ and $\tau=B$. For $r\leq r_{0}$ we have $\tau<B$. So, the asymptotic analysis predicts a region of unyielded material near the axis of symmetry. Consequently, our asymptotic expansion breaks down.

Fig. 3 shows that the values of the second invariant $\tau$ of the stress tensor along the disk surface increase with the Bingham number $B$. The solid lines, corresponding to $\tau(r,1)=\tau^{0}(r,1)+\varepsilon\tau^{1}(r,1)$, is below the dashed line, corresponding to $\tau^{0}(r,1)$. The function $\tau(r)$ increases with $r$ and decreases sharply near the axis of symmetry $r=0$. The points on the graphs for which $\tau=B$ are marked diamonds on the solid lines.

3.3 Pressure distribution and squeeze force

In order to obtain the squeeze force we first calculate the normal radial stress $\sigma_{rr}$, neglecting the terms $O(\varepsilon^{2})$:

[TABLE]

We see that $\sigma_{rr}$ is dependent on $z$, which makes it impossible to satisfy the zero-stress condition at $r=1$ exactly. We impose the average boundary condition, suggested in [8]:

[TABLE]

We substitute (48) into the average condition (49) and obtain

[TABLE]

The pressure gradient zero-order is given by $p_{0}^{\prime}(r)=-\frac{B}{z_{0}}$. The integration, using the results from (23), (34):

[TABLE]

[TABLE]

produces the following pressure distribution ($p_{0}(1)=p_{R}$):

[TABLE]

The function $p_{1}(r)$ can be found from (45) numerically with $p_{1}(1)=0$.

The solid lines show computed final states, the dotted lines denote the leading-order result (8) and the dashed line shows the higher-order prediction (9) .

Fig. 4(a) shows the pressure distribution along the disk surface for different values of $B$ . The values of $\varepsilon p_{R}$ are very small for $B=0.01,0.1$, increase with the Bingham number. We see that adding $\varepsilon p^{s,1}$ reduces the value of pressure, so the solid line ($p^{s}$) is lower than the dotted line ($p^{s,0}$). Since for $B=0.01$, $0.1$ the value of $\varepsilon p^{s,1}$ is very small additive to the $p^{s,0}$, the plots of $p^{s}$ and $p^{s,0}$ coincide. We see in Fig. 4(b) that as power-law index $n$ decreases, the pressure becomes more significant in the center of the disk and smaller at the edge of the disc. With the increase in the number of Bingham, the graphs for different power-law indices $n$ less differ among themselves.

We calculate the squeeze force by integrating the normal axial stress $\sigma_{zz}^{s}=-p^{s}(r)+\varepsilon^{2}\tau_{zz}^{s,0}(r,z)$ over the plate surface. In the lubrication solution up to second order in $\varepsilon$ the component $\tau_{zz}$ is negligible. So

[TABLE]

We calculate $F^{0}$, using (22), (55), (51), (52),

[TABLE]

To calculate $F^{1}$, we substitute the expression for $p_{1}^{\prime}$ (45) into (55), integrate by parts , and obtain

[TABLE]

Fig. 4 suggests that there is very little difference between the zeroth order and first order results for the pressure over the plate. As a result, there will be very little difference between the zeroth and first order forces. In Fig. 5 we plot separate results for the total force $F$, $F^{0}$, $|\varepsilon F^{1}|$, $\varepsilon p_{R}$ as the functions of $B$. (Since $F^{1}$ is negative and we use a logarithmic scale, we plot $|\varepsilon F^{1}|$.) All these quantities increase with increasing B, the function $\varepsilon p_{R}(B)$ grows much faster than $|\varepsilon F^{1}(B)|$. For $B<1$ the graphs $F$ and $F^{0}$ coincide, for $B>1$ the graph of $F$ is higher than the graph of $F_{0}$. The total force $F$ is larger than $F^{0}$ due to $\varepsilon p_{R}$, and this difference increases with $B$.

Fig. 6 shows the dependence of the total force on the Bingham number and the power index. We observe an interesting effect: for small Bingham numbers, the total force greater for low power indices; however, but with increasing Bingham number, the total force becomes greater for higher power indices. Because the graphs in a logarithmic scale are the same for large values of B, in Figure 6(b), we show the plots for large Bingham number in a linear scale.