3.1. Velocity and flow rate–pressure gradient relation in dimensional form

We consider low-Reynolds-number flow of an incompressible shear-thinning fluid in a narrow configuration where $h\ll r_{o}$ , which allows for the use of the lubrication approximation. Specifically, the lubrication approximation applies provided that the aspect ratio $\epsilon = h / r_o$ and the reduced Reynolds number $\epsilon {Re}$ are small, i.e. $\epsilon \ll 1$ and $\epsilon {Re} \ll 1$ (Leal Reference Leal2007).

The reduced Reynolds number $\epsilon {Re}$ is the ratio of fluid inertia, $\rho u_c^{2}/r_o$ , to viscous stress, $\eta (\dot {\gamma }) u_c/h^2$ , where $u_c=q/2 \pi h r_o$ is the characteristic velocity scale and $\rho$ is the density of the fluid. In contrast to a Newtonian fluid, where the reduced Reynolds number $\epsilon {Re}$ depends linearly on velocity (or shear rate $\dot {\gamma } \sim u_c/h$ ), the reduced Reynolds number for a shear-thinning fluid may exhibit either linear or nonlinear dependence on velocity, depending on the magnitude of the shear rate or shear stress. For low shear stresses (or shear rates), $\tau /\tau _{1/2}\ll 1$ , the viscosity of the Ellis fluid is $\eta (\dot {\gamma }) \approx \eta _0$ , so that the small reduced Reynolds number assumption requires $\epsilon {Re}_{\it small \; El}=\rho u_{c}h^2/\eta _{0}r_o \ll 1$ , similar to a Newtonian fluid. However, for high shear stresses (or shear rates), $\tau /\tau _{1/2}\gg 1$ , the viscosity of the Ellis fluid follows the power-law dependence (2.4) and is approximately $\eta (\dot {\gamma })\approx \eta _{0} (\eta _{0}\dot {\gamma }/\tau _{1/2})^{(1/n_e)-1}$ . Therefore, the small reduced Reynolds number assumption now requires that $\epsilon {Re}_{\it power\text{-}law} = (\rho u_c h^2 / \eta _0 r_o) (\eta _0 u_c / h \tau _{1/2})^{1 - (1/n_e)} \ll 1$ , which can also be expressed as $\epsilon {Re}_{\it power\text{-}law} = (\epsilon {Re}_{\it small \; El}) El^{1 - (1/n_e)} \ll 1$ , where $El = \eta _0 u_c / h \tau _{1/2}$ is the Ellis number defined in (3.10). As expected, in the power-law regime, the reduced Reynolds number increases nonlinearly with velocity (or the Ellis number) due to the shear-thinning behaviour of the fluid. Clearly, in the power-law regime where $El \gg 1$ , the condition $\epsilon {Re}_{\it power\text{-}law} = (\epsilon {Re}_{\it small \; El}) El^{1 - (1/n_e)} \ll 1$ is more restrictive than $\epsilon {Re}_{\it small \; El} \ll 1$ , since the latter does not necessarily imply that $\epsilon {Re}_{\it power\text{-}law} \ll 1$ . Nevertheless, for typical values of the physical parameters associated with the radial flow of an Ellis fluid, listed in table 4, both conditions, $\epsilon {Re}_{\it small \; El} \ll 1$ and $\epsilon {Re}_{\it power\text{-}law} \ll 1$ , are well satisfied, as shown in Appendix C.

In this lubrication limit, fluid inertia and longitudinal gradients are negligible, and we can assume a unidirectional radial flow, $\boldsymbol{u}=u_{r}(r,z)\boldsymbol{e}_{r}$ (see, e.g. Park Reference Park2020). Under this assumption, the momentum equations (2.1b ) in conjunction with (2.2) and (2.3) reduce to

(3.1a,b) \begin{align} \frac {\partial p}{\partial r}=\frac {\partial }{\partial z}\left (\eta (\tau )\frac {\partial u_{r}}{\partial z}\right)\!,\qquad \frac {\partial p}{\partial z}=0, \end{align}

where the magnitude of the shear stress $\tau =\sqrt {(\boldsymbol{\tau }:\boldsymbol{\tau })/2}$ is given by

(3.2) \begin{equation} \tau =|\tau _{rz}|=\eta (\tau )\dot {\gamma } \quad \textrm{with}\; \dot {\gamma }=\left |\frac {\partial u_{r}}{\partial z}\right |. \end{equation}

The governing equations (3.1) represent the lubrication equations that can be derived by introducing a small parameter $\epsilon =h/r_o\ll 1$ , and then formally expanding the velocity and pressure fields in powers of $\epsilon$ and considering the leading order terms (see, e.g. Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Leal Reference Leal2007). From (3.1b ), it follows that $p=p(r)$ , i.e. the pressure is constant across the gap $h$ , but varies along the $r$ -direction.

Integrating (3.1a ) with respect to $z$ and applying the symmetry condition at $z=0$ yields

(3.3) \begin{equation} \frac {\textrm{d}p}{\textrm{d}r}z=\eta (\tau )\frac {\partial u_{r}}{\partial z}. \end{equation}

Substituting (2.3) into (3.3) and rearranging, we obtain

(3.4) \begin{equation} \frac {\partial u_{r}}{\partial z}=\frac {\textrm{d}p}{\textrm{d}r}\frac {z}{\eta (\tau )}=\frac {1}{\eta _0}\frac {\textrm{d}p}{\textrm{d}r}z(1 + (\tau /\tau _{1/2})^{n_e-1}). \end{equation}

Noting that $\dot {\gamma }=-\partial u_{r}/\partial z$ for $0\leqslant z\leqslant h/2$ , and using (3.2) and (3.3), the magnitude of the shear stress $\tau$ can be expressed as

(3.5) \begin{equation} \tau =-\eta (\tau )\frac {\partial u_{r}}{\partial z}=-\frac {\textrm{d}p}{\textrm{d}r}z \quad \textrm{for}\quad 0 \leqslant z\leqslant h/2, \end{equation}

where we expect the pressure gradient to be negative, $\textrm{d}p/\textrm{d}r\lt 0$ .

Combining (3.4) and (3.5) provides the expression for $\partial u_{r}/\partial z$ . Integrating this result with respect to $z$ and applying the no-slip boundary conditions at the plate walls, we obtain an explicit expression for the radial velocity profile of an Ellis fluid,

(3.6) \begin{equation} u_r(r,z)=\frac {1}{\eta _0}\frac {\textrm{d}p}{\textrm{d}r}\left [\frac {1}{2}\left (z^2-\frac {h^2}{4}\right )+\frac {1}{n_e+1}\left (-\frac {\textrm{d}p}{\textrm{d}r} \frac {1}{\tau _{1/2}}\right )^{n_e-1}\left (z^{n_e+1}-\frac {h^{^{n_e+1}}}{2^{^{n_e+1}}}\right )\right ]\!. \end{equation}

Finally, substituting (3.6) into the integral constraint for the volumetric flux, $q=2\pi \int _{-h/2}^{h/2} u_r r \, \textrm{d}z$ provides a nonlinear equation for the pressure gradient $\textrm{d}p(r)/\textrm{d}r$ for a given flow rate $q$ ,

(3.7) \begin{equation} \frac {q}{4\pi r} = -\frac {1}{\eta _0}\frac {\textrm{d}p}{\textrm{d}r}\left [\frac {h^3}{24}+\frac {h^{n_e+2}}{2^{n_e+2}(n_e+2)}\left (-\frac {\textrm{d}p}{\textrm{d}r} \frac {1}{\tau _{1/2}}\right )^{n_e-1}\right ]\!. \end{equation}

In addition, the integral constraint for the volumetric flux, $q=2\pi \int _{-h/2}^{h/2} u_r r \, \textrm{d}z$ , implies that the radial velocity has the fully developed form of $ u_r(r,z)=f(z)/r$ (see, e.g. Park Reference Park2020). We will return to this point in § 4.

3.2. Non-dimensionalisation

To non-dimensionalise our hydrodynamic problem, we introduce dimensionless variables based on the lubrication theory,

(3.8) \begin{align} R=\frac {r}{r_o}, \quad Z=\frac {z}{h}, \quad U_r=\frac {u_r}{u_c}, \quad P=\frac {p}{\eta _0 u_c r_o /h^2}, \quad \Delta P=\frac {\Delta p}{\eta _0 u_c r_o /h^2}, \quad \mathcal{H}=\frac {\eta }{\eta _0}, \end{align}

where $u_c=q/2 \pi h r_o$ is the characteristic velocity scale. In addition, we introduce the inlet-to-outlet aspect ratio $\alpha =r_i/r_o$ , which is assumed to be small, $\alpha \ll 1$ , consistent with the experiments of Laurencena & Williams (Reference Laurencena and Williams1974).

Using the non-dimensionalisation (3.8), the governing equation for the pressure gradient (3.7) takes the form

(3.9) \begin{equation} \frac {1}{R} = -\frac {\textrm{d}P}{\textrm{d}R}\left [\frac {1}{12}+\frac {1}{2^{n_e+1}(n_e+2)}\left (-El\frac {\textrm{d}P}{\textrm{d}R} \right )^{n_e-1}\right ], \end{equation}

where $\textit{El}$ is the non-dimensional Ellis number defined as

(3.10) \begin{equation} El=\frac { \dot {\gamma }_{c}}{ \tau _{1/2}/\eta _{0}}=\frac {\eta _{0} u_{c}}{h \tau _{1/2}}=\frac {\eta _{0} q}{2 \pi h^2 r_o \tau _{1/2}}. \end{equation}

The Ellis number is the ratio of the characteristic shear rate in the flow, $\dot {\gamma }_{c}=u_{c}/h$ , to the cross-over shear rate of the fluid, $\tau _{1/2}/\eta _{0}$ , representing the relative importance of the shear thinning. For $El\ll 1$ , the onset of shear-thinning effect occurs only at sufficiently high shear rates, and the fluid is expected to behave as Newtonian. In contrast, when $El\gg 1$ , we expect the shear-thinning effect to become apparent at low shear rates, and the shear-thinning fluid follows the power-law model. We note that the definition of the Ellis number, (3.10), is similar to the definition of the Carreau number, provided $\tau _{1/2}=\eta _{0}/\lambda$ , where $\lambda$ is the inverse of a characteristic shear rate of the Carreau fluid (see, e.g. Boyko & Stone Reference Boyko and Stone2021b ; Chun et al. Reference Chun, Boyko, Christov and Feng2024). However, our definition of the Ellis number is inverse to the definition used in recent studies of Picchi et al. (Reference Picchi, Ullmann, Brauner and Poesio2021) and Chun et al. (Reference Chun, Ji, Yang, Malik and Feng2022).

Using the non-dimensionalisation (3.8), the Ellis and power-law models, (2.3) and (2.4), can be expressed in the dimensionless form in terms of the Ellis number,

(3.11) \begin{equation} \mathcal{H}_{ \it Ellis} = \frac {1}{1 + ( \mathcal{H}_{\it Ellis} El)^{n_e-1}} \quad \text{and} \quad \mathcal{H}_{ \it power\text{-}law}=El^{\frac {1}{n_e}-1},\quad \text{with} \; El=\frac { \dot {\gamma }}{ \tau _{1/2}/\eta _{0}}. \end{equation}

We present in figure 2 the dimensionless viscosity $\mathcal{H}=\eta /\eta _0$ as a function of the Ellis number $El= \eta _{0} \dot {\gamma }/ \tau _{1/2}$ for $n_e=1.5, 2.5$ and 5. As expected, for small values of $El$ , we recover the Newtonian plateau with $\mathcal{H}=1$ , whereas for large values of $El$ , we obtain the power-law behaviour with $\mathcal{H}=El^{({1}/{n_e})-1}$ , represented by red dashed lines. It is evident that as the index $n_e$ characterising the degree of shear thinning increases, the viscosity transitions from the small- $El$ (Newtonian) to the power-law behaviour at smaller values of the Ellis number.

Figure 2. Non-dimensional viscosity $\mathcal{H}=\eta /\eta _0$ as a function of the Ellis number $El= \eta _{0} \dot {\gamma }/ \tau _{1/2}$ for different values of the shear-thinning index $n_e$ . Red dashed lines represent the power-law asymptotic limit valid for large values of $El$ .

For arbitrary values of the Ellis number, the pressure gradient $\textrm{d}P/\textrm{d}R \equiv P'(R)$ can be calculated numerically from (3.9), as detailed in Appendix A. Once the pressure gradient $P'(R)$ is known, the pressure distribution $P(R)$ and pressure drop $\Delta P\equiv P(R=\alpha )-P(R=1)$ , expressed in terms of the integrals,

(3.12a,b) \begin{align} P(R) = -\int _R^1 P'(R) \,\textrm {d}R \quad \text{and} \quad \Delta P =P(\alpha )= -\int _\alpha ^1 P'(R) \textrm{d}R, \end{align}

can be obtained using numerical integration (see Appendix A). Note that, in (3.12), we have used the fact that $P(1)=0$ .

Furthermore, having obtained the pressure gradient $\textrm{d}P/\textrm{d}R$ , we can determine the non-dimensional radial velocity of an Ellis fluid using (3.6) as

(3.13) \begin{equation} U_r(R,Z)=\frac {\textrm{d}P}{\textrm{d}R}\left [\frac {1}{2}\left (Z^2-\frac {1}{4}\right )+\frac {1}{n_e+1}\left (-El\frac {\textrm{d}P}{\textrm{d}R}\right )^{n_e-1}\left (Z^{n_e+1}-\frac {1}{2^{^{n_e+1}}}\right )\right ]\!. \end{equation}

3.3. Asymptotic results for small and large values of Ellis number

Although no closed-form analytical solution is available for (3.9), in this subsection, we provide asymptotic solutions in the limit of small and large values of $El$ , thus allowing us to describe analytically almost the entire range of Ellis numbers.

For weak actuation, corresponding to the limit $El\ll 1$ , we seek the solution for the pressure gradient $\textrm{d}P/\textrm{d}R \equiv P'(R)$ in the form

(3.14) \begin{equation} P'(R)=P_0^{\prime}(R)+El^{n_e-1}P_1^{\prime}(R)+El^{2(n_e-1)}P_2^{\prime}(R)+O(El^{3(n_e-1)}). \end{equation}

Note that since $n_e\geqslant 1$ , we have $El^{n_e-1}\ll 1$ and $El^{2(n_e-1)}\ll 1$ in the small- $El$ limit. Substituting this expansion into (3.9) and solving order by order, we obtain

(3.15) \begin{equation} P_0^{\prime}(R)=-\frac {12}{R},\quad P_1^{\prime}(R)=\frac {6^{n_e +1}}{n_e+2}\left (\frac {1}{R}\right )^{n_e}, \quad P_2^{\prime}(R)=-\frac {3^{2n_e+1} 4^{n_e} n_e }{(n_e+2)^2} \left (\frac {1}{R}\right )^{2n_e-1}. \end{equation}

Integrating (3.15) with respect to $R$ and using the boundary condition $P(R=1)=0$ provides the expression for the pressure,

(3.16) \begin{align} P(R)&=\underset {P_{0}(R)}{\underbrace {-12\ln R}}+El^{n_{e}-1}\underset {P_{1}(R)}{\underbrace {\frac {6^{n_{e}+1}}{(n_{e}+2)(n_{e}-1)}(1-R^{1-n_{e}})}} \nonumber \\ &\quad + El^{2(n_{e}-1)}\underset {P_{2}(R)}{\underbrace {\frac {2^{2n_{e}-1}3^{2n_{e}+1}n_{e}}{(n_{e}+2)^{2}(n_{e}-1)}(R^{2-2n_{e}}-1)}}\quad \textrm{for}\; El\ll 1. \end{align}

Thus, the pressure drop $\Delta P$ in the small- $El$ limit is

(3.17) \begin{align} \Delta P &= \underset {\Delta P_{0}}{\underbrace {-12\ln \alpha }}+El^{n_{e}-1}\underset {\Delta P_{1}}{\underbrace {\frac {6^{n_{e}+1}}{(n_{e}+2)(n_{e}-1)}(1-\alpha ^{1-n_{e}})}} \nonumber \\ &\quad + El^{2(n_{e}-1)}\underset {\Delta P_{2}}{\underbrace {\frac {2^{2n_{e}-1}3^{2n_{e}+1}n_{e}}{(n_{e}+2)^{2}(n_{e}-1)}(\alpha ^{2-2n_{e}}-1)}}\quad \textrm{for}\; El\ll 1 .\end{align}

For strong actuation, corresponding to the limit of $El\gg 1$ , the Ellis model reduces to the power-law model and we expect to achieve a power-law regime. Specifically, the first term in the brackets on the right-hand side of (3.9) is negligible compared with the second term, so that the pressure gradient $ P'(R)$ can be expressed as

(3.18) \begin{equation} P'(R)=-2^{(n_{e}+1)/n_{e}}(n_{e}+2)^{1/n_{e}}\frac {El^{(1-n_{e})/n_{e}}}{R^{1/n_{e}}} \quad \textrm{for}\; El\gg 1. \end{equation}

Integrating (3.18) with respect to $R$ and using the boundary condition $P(R=1)=0$ yields the expression for the pressure in the large-Ellis-number limit

(3.19) \begin{equation} P(R)=-\frac {2^{(n_{e}+1)/n_{e}} n_e (n_e+2)^{1/n_e} }{ n_e-1}El^{(1-n_{e})/n_{e}} (R^{(n_e-1)/n_e}-1)\quad \textrm{for}\; El\gg 1, \end{equation}

so that the pressure drop in the large- $El$ limit is

(3.20) \begin{equation} \Delta P=-\frac {2^{(n_{e}+1)/n_{e}} n_e (n_e+2)^{1/n_e} }{ n_e-1}El^{(1-n_{e})/n_{e}} (\alpha ^{(n_e-1)/n_e}-1)\quad \textrm{for}\; El\gg 1. \end{equation}

We note that, as expected, (3.20), when written in a dimensional form, agrees with the expression for the pressure drop of the power-law fluid (see, e.g. Na & Hansen Reference Na and Hansen1967; Laurencena & Williams Reference Laurencena and Williams1974). Furthermore, for $n_e=1$ , (3.20) reduces to $\Delta P=-12\ln {\alpha }$ , consistent with the Newtonian limit.

4.1. Pressure drop and pressure distribution

First, in figure 3(a,b), we present the non-dimensional pressure drop $\Delta P=\Delta p/(\eta _{0}q/2 \pi h^{3})$ as a function of $El$ for the radial flow of an Ellis fluid for $n_e=1.5$ in panel (a) and $n_e=2.5$ in panel (b), with $\alpha = 0.05$ . Black triangles represent the finite-element simulation results obtained from calculating the pressure drop along the midplane $Z=0$ . Grey circles represent the theoretical results obtained by solving numerically (3.9) and (3.12b ). Cyan dotted and red dashed curves represent the asymptotic solutions (3.17) and (3.20) for small and large values of $El$ , respectively. Clearly, there is excellent agreement between the theoretical predictions based on our reduced-order model and finite-element simulation results. We observe that, as expected, the power-law asymptotic solution (3.20) captures fairly well the variation of the pressure drop with $El$ for large values of the Ellis number. However, the small- $El$ asymptotic solution (3.17) cannot accurately capture the pressure drop except for very small values of $El$ . Nevertheless, using small- $El$ asymptotic expressions for $\Delta P_{0}$ , $\Delta P_{1}$ and $\Delta P_{2}$ given in (3.17), we can improve the convergence of the asymptotic series by applying the diagonal Padé [1/1] approximation (Hinch Reference Hinch1991; Housiadas Reference Housiadas2017),

(4.1) \begin{equation} \Delta P_{ \it Pade}=\Delta P_{0}+El^{n_{e}-1} \frac {(\Delta P_{1})^2}{\Delta P_{1}-El^{n_{e}-1}\Delta P_{2}}\quad \textrm{for}\; El\ll 1, \end{equation}

represented in figure 3(a,b) by purple dash-dotted curves.

Figure 3. Non-dimensional pressure drop of the shear-thinning Ellis fluid in a radial flow between two disk-shaped plates. (a,b) Dimensionless pressure drop $\Delta P=\Delta p/(\eta _{0}q/2 \pi h^{3})$ as a function of $El=\eta _{0} q/(2 \pi h^2 r_o \tau _{1/2})$ for (a) $n_e=1.5$ and (b) $n_e=2.5$ . Black triangles represent the results of the finite-element simulation. Grey circles represent the theoretical results obtained by solving numerically (3.9) and (3.12b ). Cyan dotted curves represent the small- $El$ asymptotic solution (3.17). Purple dash-dotted curves represent the Padé approximation (4.1). Red dashed lines represent the power-law asymptotic solution (3.20). All calculations were performed using $\alpha = 0.05$ .

It is evident from figure 3(a,b) that the Padé approximation (4.1) significantly improves the agreement with the semi-analytical and finite-element simulation results for both $n_e=1.5$ and 2.5. For example, in the case of $n_e=1.5$ , we have a modest relative error of approximately 8 % for up to $El=0.04$ . As expected, when $El$ increases, the agreement between the Padé approximation (4.1) and simulations deteriorates, and the small- $El$ Padé approximation overpredicts the pressure drop, yielding a relative error of approximately 18 % for $El=0.115$ . Nevertheless, our theoretical predictions based on (3.9) and (3.12b ) are in excellent agreement with simulations even for order-one Ellis numbers, resulting in a relative error below 0.5 % for $El=0.87$ and $n_e=1.5$ . Therefore, our theory accurately describes the pressure drop behaviour and captures well the transition to the power-law dependence $\Delta P \sim El^{({1}/{n_e})-1}$ , given in (3.20), for larger values of $El$ . The dimensionless pressure drop, which physically represents the dimensionless flow resistance ( $\Delta p/q$ ) of the Ellis fluid, monotonically decreases with $El$ . Interestingly, for $El\gg 1$ , the scaling for the pressure $\Delta P \sim El^{({1}/{n_e})-1}$ is identical to the power-law dependence of the dimensionless viscosity $\mathcal{H}=El^{({1}/{n_e})-1}$ , thus highlighting that the pressure drop reduction for shear-thinning fluids is due to the decrease in viscosity $\eta$ with the applied shear rate $\dot {\gamma }$ or the shear stress $\tau$ .

Next, we present in figure 4(a,b) the non-dimensional pressure distribution $ P=2\pi p h^3 / (\eta _0 q)$ as a function of the radial coordinate $R=r/r_o$ for the radial flow of an Ellis fluid for $n_e=1.5$ in panel (a) and $n_e=2.5$ in panel (b), with $El = 0.05$ and $El = 25$ . Similar to the pressure drop, our theoretical predictions for the pressure (grey solid curves) accurately capture the finite-element simulation results (black dots) for all considered values of $El$ . Furthermore, we observe an excellent agreement between the small- and large- $El$ asymptotic solutions (purple dash-dotted and red dashed curves) and the finite-element simulation results. Note that the small- $El$ asymptotic solution is based on the application of the Padé approximation to the expressions in (3.16), similar to (4.1).

Figure 4. Non-dimensional pressure distribution of the shear-thinning Ellis fluid in a radial flow between two disk-shaped plates. (a,b) Dimensionless pressure $ P=2\pi p h^3 / (\eta _0 q)$ as a function of the radial coordinate $R=r/r_o$ for small ( $El=0.05$ ) and large ( $El=25$ ) Ellis numbers with (a) $n_e=1.5$ and (b) $n_e=2.5$ . Black dots represent the results of the finite-element simulation. Grey curves represent the theoretical results obtained by solving numerically (3.9) and (3.12a ). Purple dash-dotted curves represent the small- $El$ Padé approximation based on (3.16). Red dashed curves represent the power-law asymptotic solution (3.19). All calculations were performed using $\alpha = 0.05$ .

As expected, the results in figure 4(a,b) reveal that, for $El \ll 1$ , the pressure distribution shows a weak dependence on the shear-thinning index $n_e$ . However, for $El \gg 1$ , the pressure magnitude decreases as $n_e$ increases, consistent with the results shown in figure 2, indicating that the viscosity reduction for a given $El$ in the power-law regime becomes more pronounced with $n_e$ .

4.2. Comparison between the theoretical predictions and the finite-element simulations for the radial velocity

Our theoretical analysis presented in § 3 relies on the assumption of a unidirectional radial flow, $\boldsymbol{U}=U_{r}(R,Z)\boldsymbol{e}_{r}$ , which together with the integral constraint for the volumetric flux, $\int _{-1/2}^{1/2} U_r R \, \textrm{d}Z=1$ , implies that the radial velocity takes the fully developed form of $ U_r(R,Z)=\mathcal{F}(Z)/R$ for $\alpha \leqslant R\leqslant 1$ (see, e.g. Na & Hansen Reference Na and Hansen1967; Park Reference Park2020). Clearly, the assumption of unidirectional radial flow is not valid near the inlet, where the shear-thinning fluid enters from the tube into the gap between the plates. Indeed, even for a Newtonian fluid at low Reynolds numbers, there is a radial entrance length required to reach a fully developed velocity profile $U_{r}(R,Z)=-({3}/{2R})(4Z^{2}-1)$ (Chatterjee Reference Chatterjee1993, Reference Chatterjee2000). Therefore, it is of particular interest to compare our theoretical predictions for the radial velocity of the shear-thinning Ellis fluid with the results of the numerical simulations and elucidate over which length scale the radial velocity reaches its fully developed form, $ U_r(R,Z)=\mathcal{F}(Z)/R$ .

Figure 5. Non-dimensional radial velocity distribution of the shear-thinning Ellis fluid in a radial flow between two disk-shaped plates. (a–d) Contour plot of the radial velocity distribution, $U_r$ , as a function of the $(R,Z)$ coordinates for (a,b) $\alpha =0.1$ and (c,d) $\alpha =0.05$ , obtained from (a,c) our theory and (b,d) finite-element simulations. White dashed lines represent the midplane $Z=0$ . (e, f) Radial velocity $U_r$ multiplied by the radial coordinate $R$ , $U_r R$ , as a function of $Z$ for $R=\alpha , 2\alpha$ and $3\alpha$ , with (e) $\alpha =0.1$ and (f) $\alpha =0.05$ . Red solid lines represent the theoretical results obtained by solving (3.9) and (3.13). Dots, crosses and circles represent the results of the finite-element simulation. (g) Relative error between the theory and numerical simulations for the radial velocity along the midplane, $U_r(R,Z=0)$ , as a function of the radial coordinate $R$ for $\alpha =0.1$ (grey circles) and $\alpha =0.05$ (black dots). All calculations were performed using $El=0.1$ and $n_e=2.5$ .

In figure 5, we present the results of the dimensionless radial velocity of the shear-thinning Ellis fluid in a radial flow for $n_e=2.5$ and the Ellis number $El=0.1$ , corresponding to the transition from the small- $El$ limit to the power-law regime, as shown in figure 3(b). Figure 5(a–d) shows a comparison of theoretical predictions (figures 5 a and 5 c) and finite-element simulation results (figures 5 b and 5 d) for contours of the radial velocity, $U_{r}$ , as a function of the $(R,Z)$ coordinates for the Ellis fluid with $\alpha =0.1$ in panel (a,b) and $\alpha =0.05$ in panel (c,d). We observe a very good agreement between the theoretical and numerical results for the radial velocity $U_{z}$ for both $\alpha =0.1$ and $\alpha =0.05$ in the entire domain. However, as expected, near the inlet, the agreement deteriorates and for $\alpha =0.1$ , we see that the theoretically predicted radial velocity slightly overpredicts the magnitude of $U_r(R=\alpha ,Z=0)$ obtained from the simulation results.

To understand how fast the radial velocity attains its fully developed form, $ U_r(R,Z)=\mathcal{F}(Z)/R$ , we present in figure 5(e, f) the scaled radial velocity $U_r R$ as a function of $Z$ for $R=\alpha , 2\alpha$ and $3\alpha$ , with $\alpha =0.1$ in panel (e) and $\alpha =0.05$ in panel ( f). Red solid lines represent our theoretical predictions, and dots, crosses and circles represent the finite-element simulation results. It is evident that the radial velocity at $R=\alpha$ (grey circles) is not symmetric with respect to the midplane $Z=0$ and deviates from the theoretically predicted fully developed profile due to the entering fluid flow from the tube into the gap between the plates. However, already at $R=2\alpha$ , we observe that the radial velocity becomes symmetric and our theory captures well the velocity profile for both $\alpha =0.1$ and $\alpha =0.05$ . For further clarification, figure 5(g) presents the relative error between the theory and numerical simulations for the radial velocity along the midplane, $U_r(R,Z=0)$ , as a function of the radial coordinate $R$ for $\alpha =0.1$ (grey circles) and $\alpha =0.05$ (black dots). It follows from figure 5(g) that the relative error, defined as $|(U_{r}^{\it theory}-U_{r}^{ \it sim})/U_{r}^{ \it sim}|_{Z=0}$ , decays very fast and, for both $\alpha =0.1$ and $\alpha =0.05$ , it is below 0.5 % when $R\gtrsim 0.105$ .

So far, in figure 5, we have presented results for the radial velocity of a shear-thinning Ellis fluid with $El=0.1$ and $n_e=2.5$ . To elucidate the influence of the Ellis number and the shear-thinning index on the development of the radial velocity, we present in figure 6 the radial entrance length $R_{\it entry}$ (scaled by $\alpha =0.05$ ), defined as the radial position $R$ where the relative error $|(U_{r}^{\it theory}-U_{r}^{ \it sim})/U_{r}^{ \it sim}|_{Z=0}$ falls below 1 %, as a function of $El$ for three different values of the shear-thinning indices $n_e=1.5, 2.5$ and 5. We observe that for $El\ll 1$ , the entrance length is independent of the shear-thinning index $n_e$ and has a constant value of $R_{\it entry}/\alpha =1.05$ . However, as the Ellis number increases, the entrance length $R_{\it entry}$ exhibits a strong dependence on $n_e$ . While for $n_e=1.5$ , the entrance length remains nearly constant throughout the investigated range of Ellis numbers, for higher shear-thinning indices, the entrance length $R_{\it entry}$ starts to increase with $El$ above a certain value of the Ellis number, as observed for $n_e=2.5$ and $n_e=5$ .

We note that our small- $El$ results for the shear-thinning Ellis fluid are consistent with the previous estimate for the radial entrance length $R_{\it entry}$ of a Newtonian fluid at low Reynolds numbers. Specifically, Chatterjee (Reference Chatterjee2000) studied the radial entrance flow of a Newtonian fluid, considering a semi-infinite configuration in the radial direction, and defined the entrance length $R_{\it entry}$ as the value of $R$ at which the midplane radial velocity $U_r(R,Z=0)$ is at 1 % from its fully developed value. In particular, Chatterjee (Reference Chatterjee2000) found that $R_{\it entry}=1.05\alpha$ for $h/r_i=0.25$ in the absence of inertia. In our finite-element simulations, we have $\epsilon =0.01$ and for $\alpha =0.05$ , we find that $R_{\it entry}\approx 1.05\alpha =0.0525$ for $\epsilon /\alpha =h/r_i=0.2$ and $El\leqslant 10^{-2}$ , in close agreement with the Newtonian value of Chatterjee (Reference Chatterjee2000).

Figure 6. Radial entrance length $R_{\it entry}$ , defined as the radial position $R$ where the relative error between the theory and numerical simulations for the radial velocity along the midplane, $|(U_{r}^{\it theory}-U_{r}^{ \it sim})/U_{r}^{ \it sim}|_{Z=0}$ , falls below 1 %, as a function of $El$ for three values of the shear-thinning indices $n_e=1.5, 2.5$ and 5. All calculations were performed using $\alpha =0.05$ .

To better understand the variation of the entrance length $R_{\it entry}$ with the Ellis number and the shear-thinning index, we present in figure 7(a) the scaled radial velocity $U_r R$ , evaluated at $(R,Z)=(0.1,0)$ , as a function of $El$ for three different values of the shear-thinning indices $n_e=1.5, 2.5$ and 5. Solid lines represent the theoretical results obtained by solving (3.9) and (3.13). Dots, crosses and circles represent the results of the finite-element simulation. The cyan dotted line represents the small- $El$ (Newtonian) asymptotic solution,

(4.2) \begin{equation} U_r(R,Z)=-\frac {6}{R}\left (Z^{2}-\frac {1}{4}\right )\quad \textrm{for}\; El\ll 1, \end{equation}

and red dashed lines represent the power-law asymptotic solution,

(4.3) \begin{equation} U_r(R,Z)=-\frac {2^{n_e+1}(n_e+2)}{(n_e+1)R}\left (Z^{n_e+1}-\frac {1}{2^{^{n_e+1}}}\right )\quad \textrm{for}\; El\gg 1, \end{equation}

obtained from substituting (3.18) into (3.13) while neglecting the first term in the brackets of (3.13).

It is evident from figure 7(a) that the radial velocity shows a weak dependence on $El$ and $n_e$ for small Ellis numbers, consistent with the small- $El$ asymptotic solution (4.2). Furthermore, the radial velocity exhibits a weak dependence on $El$ for large Ellis numbers, yet retains a dependence on $n_e$ , in agreement with the power-law asymptotic solution (4.3). However, in the transitional range of Ellis numbers, from the small- $El$ (Newtonian) regime to the power-law behaviour, the radial velocity depends on both the Ellis number and the shear-thinning index.

The weak dependence of $U_r$ on $El$ and $n_e$ for small Ellis numbers rationalises the weak dependence of the entrance length $R_{\it entry}$ on the fluid rheology for $El \ll 1$ , as shown in figure 6. This is also highlighted in figure 7(b), which shows the relative error between the theoretical prediction and numerical simulations for the radial velocity at $(R, Z) = (0.1, 0)$ as a function of $El$ for three different values of the shear-thinning indices $n_e=1.5, 2.5$ and 5. Clearly, for $n_e = 1.5$ and $n_e = 2.5$ , the relative error at $(R, Z) = (0.1, 0)$ remains below 1 % when $El \leqslant 0.1$ , supporting the conclusion that the entrance length $R_{\it entry}$ is smaller than 0.1 in this regime. However, for $n_e = 5$ , the relative error exceeds 1 % when $El \geqslant 0.1$ , indicating that the entrance length is greater than 0.1 in this case.

Figure 7. (a) Radial velocity $U_r$ multiplied by the radial coordinate $R$ , evaluated at $(R,Z)=(0.1,0)$ , as a function of $El$ for three values of the shear-thinning indices $n_e=1.5, 2.5$ , and 5. Solid lines represent the theoretical results obtained by solving (3.9) and (3.13). Dots, crosses and circles represent the results of the finite-element simulation. The cyan dotted line represents the small- $El$ asymptotic solution (4.2). Red dashed lines represent the power-law asymptotic solution (4.3). (b) Relative error between the theory and numerical simulations for the radial velocity at $(R,Z)=(0.1,0)$ as a function of $El$ for three values of the shear-thinning indices $n_e=1.5, 2.5$ and 5. All calculations were performed using $\alpha =0.05$ .

5.1. Fit of viscosity data and rheological parameters of the Ellis model

Laurencena & Williams (Reference Laurencena and Williams1974) provided shear-rate-dependent viscosity data for the Natrosol solution measured with a Weissenberg rheogoniometer (see their figure 3). In addition, they used a falling ball viscometer to measure the zero-shear-rate viscosity $\eta _0$ . We reproduce in figure 8 the viscosity data from Laurencena & Williams (Reference Laurencena and Williams1974) for the aqueous Natrosol solution, obtained using a Weissenberg rheogoniometer. The black curve represents the fit of the complete set of the rheological data to the Ellis model (2.3) obtained using MATLAB’s nonlinear least-squares routine lsqcurvefit. The corresponding rheological parameters are provided in the upper row of table 2.

Figure 8. Experimental data from Laurencena & Williams (Reference Laurencena and Williams1974) and the fitting curves for viscosity as a function of shear rate for the Natrosol solution. The black curve represents the fit of the complete set of the rheological data () to the Ellis model (2.3). The red dashed curve represents the fit of the rheological data to the power-law model (2.4). The rheological parameters obtained from the fitting are summarised in table 2.

Table 2. Rheological parameters obtained from fitting the viscosity dependence on the shear rate of the Natrosol solution to the Ellis and power-law models based on the complete data set (upper row) and based solely on the data from the power-law regime (middle row). The lower row includes the fitting parameters reported by Laurencena & Williams (Reference Laurencena and Williams1974), who also determined $\eta _0$ using a falling ball viscometer at 21.5 $^\circ$ C. Certainty bounds are calculated based on the method presented by Ashkenazi & Solav (Reference Ashkenazi and Solav2025).

We observe that the leftmost data point in figure 8 corresponds to the transition from the low-shear-rate viscosity limit to the power-law regime. The red dashed curve represents the fit of the rheological data without this data point to the power-law model (2.4). The resulting rheological parameters, provided in the middle row of table 2, are in close agreement with the parameters from Laurencena & Williams (Reference Laurencena and Williams1974) for the power-law regime, as given in the lower row. It follows from table 2 that $n\approx 1/n_e$ , consistent with (2.4). While Laurencena & Williams (Reference Laurencena and Williams1974) used a falling ball viscometer to obtain the zero-shear-rate viscosity $\eta _0=38.0$ Pa s (see the lower row in table 2), unfortunately, they did not provide the error bars nor shear rate conditions for this measurement. Therefore, in our calculations, we use the value of $\eta _0$ obtained from the fitting to the Ellis model, which is a bit larger.

5.2. A quantitative comparison between theory and experiments

First, we summarise in table 3 the values of the physical parameters of the experimental system of Laurencena & Williams (Reference Laurencena and Williams1974) used for the measurements of the pressure distribution of the shear-thinning Natrosol solution as a function of the radial position for two different flow rates and heights. We observe that experimental values satisfy the lubrication assumptions of our theoretical analysis, $\epsilon =h/r_o\ll 1$ .

Table 3. Values of the physical parameters for the experimental system of Laurencena & Williams (Reference Laurencena and Williams1974), which measured the pressure distribution of the shear-thinning Natrosol solution as a function of the radial position for two different flow rates and heights. The values of $El$ are calculated using $\eta _{0}$ and $\tau _{1/2}$ from the upper row of table 2, corresponding to fitting the viscosity dependence on the shear rate to the Ellis model.

Next, we present in figure 9(a,b) a comparison of our theoretical, asymptotic and finite-element simulation results with the experimental measurements of Laurencena & Williams (Reference Laurencena and Williams1974) for the pressure distribution of the shear-thinning Natrosol solution. Figure 9(a) represents the experiment 1 (the upper row in table 3), corresponding to $El=0.233$ , and figure 9(b) represents the experiment 2 (the lower row in table 3), corresponding to $El=26.34$ . Looking first at figure 9(b), we note the corresponding Ellis number $El=26.34$ is within the power-law regime. In this case, we observe excellent agreement between our theoretical (grey curve) and finite-element simulation (black dots) predictions and the experimental results (grey circles). Furthermore, the power-law asymptotic solution (3.19) (red dashed curve) is indistinguishable from the full theoretical solution.

Figure 9(a) presents the pressure distribution as a function of the radial position for $El=0.233$ , corresponding to a transitional Ellis number from the small- $El$ (Newtonian) to the power-law behaviour, consistent with the results shown in figures 2 and 3(b) for $n_e=2.5$ . First, there is good agreement between our theoretical predictions based on solving (3.9) and (3.12a ) and the finite-element simulation results. However, neither the small- $El$ Padé approximation nor the power-law asymptotic solution (3.19) can accurately capture the numerically predicted pressure variation at this transitional Ellis number, consistent with the results shown in figure 3(b).

Figure 9. Comparison between our theory and the experimental data from Laurencena & Williams (Reference Laurencena and Williams1974) for the pressure distribution of the shear-thinning Natrosol solution in a radial flow between two disk-shaped plates. (a,b) Pressure $p$ as a function of the radial coordinate $R=r/r_o$ for (a) experiment 1 and (b) experiment 2, whose details are summarised in table 3. Grey circles represent the experimental data, black dots represent the finite-element simulation results and grey curves represent the theoretical predictions. The purple dash-dotted curve represents the small- $El$ Padé approximation based on (3.16) and red dashed curves represent the power-law asymptotic solution (3.19). The shaded regions indicate a 20 % uncertainty in the value of the gap $h$ in the theoretical calculations due to the difficulties in measuring $h$ in the experiments of Laurencena & Williams (Reference Laurencena and Williams1974).

Second, and more importantly, the experimental pressure distribution of the Natrosol solution, shown in figure 9(a), systematically overpredicts our theoretical and numerical results for $El = 0.233$ , which are based on the inelastic shear-thinning Ellis model. We attribute this discrepancy to several factors. The first reason is the complex rheology of the non-Newtonian Natrosol solution, which is not a purely shear-thinning fluid but exhibits some elasticity, as reported by Laurencena & Williams (Reference Laurencena and Williams1974). While shear thinning tends to decrease the pressure, fluid viscoelasticity has the opposite effect, increasing it. Therefore, the observed overprediction in figure 9(a) may be attributed to fluid viscoelasticity, which tends to increase the pressure. As noted by Co & Bird (Reference Co and Bird1977), the net effect on pressure generally depends on the relative magnitudes of shear-thinning and viscoelasticity effects, provided that fluid inertia is negligible. Therefore, given the excellent agreement between theory and experiments in the power-law regime, as shown in figure 9(b), we expect the shear-thinning effect to dominate, while the effect of fluid elasticity is expected to be weak in this case. The second and major reason for the discrepancy is the uncertainty in the gap height $h$ . Indeed, Laurencena & Williams (Reference Laurencena and Williams1974) reported the difficulties in accurately measuring $h$ in the experiments due to a lack of rigidity in the support bearings. In experiment 1, shown in figure 9(a), the values of the flow rate $q$ and height $h$ are much smaller than in experiment 2, which we believe may increase uncertainty in the measured value of $h$ . Therefore, we added a shaded region, based on taking a 20 % range from the reported experimental value of $h$ . Clearly, we obtain a much better agreement between our theoretical and numerical predictions and the experimental results. We note that Laurencena & Williams (Reference Laurencena and Williams1974) did not provide the error bars for the pressure measurements, thus making it difficult to elucidate the error between the theory and experiments.

Finally, we note that, in addition to the pressure distribution measurements, Laurencena & Williams (Reference Laurencena and Williams1974) performed the flow rate–pressure drop measurements for the Natrosol solution. However, unfortunately, they did not provide the values of $h$ , which varied for each flow rate in the experiments. Thus, it is difficult to directly compare their experimental pressure drop data with our theoretical predictions.