1. Introduction
It is well known that a purely oscillatory boundary condition typically induces a velocity field with a non-zero time average, due to nonlinear effects in the underlying equations of motion. The time-average flow is referred to as the steady streaming, and thorough reviews of this phenomenon for the case of Newtonian fluids are provided by Riley (Reference Riley1967, Reference Riley2001). A general methodology for studying the steady streaming in a viscoelastic fluid, which is of interest in the present paper, is described by Nikolakis (Reference Nikolakis1992). Steady streaming can often be important, because it can play a major role in mass transport processes; indeed this is the main reason why steady streaming flows have received so much attention in the literature.
The present investigation is motivated by the need to understand the dynamics of the vitreous humour in the human eye induced by eye rotations. The vitreous chamber is the largest chamber of the eye, delimited anteriorly by the lens and posteriorly by the retina. It is filled with vitreous humour, which in young healthy subjects can be modelled as a viscoelastic fluid (Lee, Litt & Buchsbaum Reference Lee, Litt and Buchsbaum1992, Reference Lee, Litt and Buchsbaum1994; Nickerson et al. Reference Nickerson, Park, Kornfield and Karageozian2008; Swindle, Hamilton & Ravi Reference Swindle, Hamilton and Ravi2008; Sharif-Kashani et al. Reference Sharif-Kashani, Hubschman, Sassoon and Kavehpour2011). A recent review of the state of the art in modelling the flow of vitreous humour is provided by Siggers & Ethier (Reference Siggers and Ethier2012).
The dynamics of the vitreous humour induced by eye rotations have important implications for two main reasons. Firstly, there are numerous indications that vitreous stresses on the retina during eye rotations play a role in the mechanisms leading to retinal detachment. Secondly, mass transport processes within the vitreous chamber are significantly affected by the motion of the vitreous humour, with important implications for drug delivery procedures.
There are few in vivo experimental studies on the motion of the vitreous body induced by eye rotations. Piccirelli et al. (Reference Piccirelli, Bergamin, Landau, Boesiger and Luechinger2012) used magnetic resonance imaging (MRI) to measure vitreous velocity fields under controlled sinusoidal rotations of the eye bulb. Rossi et al. (Reference Rossi, Querzoli, Pasqualitto, Iossa, Placentino, Repetto, Stocchino and Ripandelli2012) employed the particle image velocimetry (PIV) technique to measure the vitreous velocity field on planes orthogonal to the axis of rotation.
On the modelling side, Repetto, Siggers & Stocchino (Reference Repetto, Siggers and Stocchino2008) studied, using both analytical and experimental methods, the steady streaming component of the flow in a viscous fluid within a sphere that is performing small-amplitude, sinusoidal, torsional oscillations, and they found excellent agreement between their theoretical predictions and experimental measurements. The overall flow is axisymmetric, and its steady streaming component consists of two circulation cells, one in each hemisphere. Particles close to the axis of rotation move along the axis to the poles, then around the wall to the equator of the sphere, finally moving radially across the equatorial plane to the centre of the sphere before moving along the axis again.
Meskauskas, Repetto & Siggers (Reference Meskauskas, Repetto and Siggers2011) studied the motion of a viscoelastic fluid in the same scenario, considering small-amplitude oscillations and finding the leading-order oscillatory flow. They showed that, using the rheological properties derived by various experimental studies in the literature, it is possible for resonant excitation of the motion of the vitreous humour to occur for frequencies typical of real eye rotations. Under conditions of resonance, there may be especially large velocities within the domain and therefore correspondingly large shear stresses on the wall.
In reality the vitreous chamber is not perfectly spherical, and the details of its geometry have a significant effect on the steady streaming component of the flow. Various authors have accounted for the lack of sphericity, by considering a weak departure from the spherical shape analytically (Repetto Reference Repetto2006; Repetto, Siggers & Stocchino Reference Repetto, Siggers and Stocchino2010; Meskauskas, Repetto & Siggers Reference Meskauskas, Repetto and Siggers2012), or by considering a realistic geometry either numerically (Balachandran & Barocas Reference Balachandran and Barocas2011; Abouali et al. Reference Abouali, Modareszadeh, Ghaffariyeh and Tu2012; Modareszadeh et al. Reference Modareszadeh, Abouali, Ghaffarieh and Ahmadi2012) or experimentally (Stocchino, Repetto & Cafferata Reference Stocchino, Repetto and Cafferata2007; Stocchino, Repetto & Siggers Reference Stocchino, Repetto and Siggers2010; Bonfiglio et al. Reference Bonfiglio, Repetto, Siggers and Stocchino2013).
In this study, we extend the work of Repetto et al. (Reference Repetto, Siggers and Stocchino2008) to the case of a viscoelastic fluid filling the chamber, closely following the approach of Böhme (Reference Böhme1992), who derived in detail a method for studying the steady streaming in a viscoelastic fluid. In spite of the likely importance of the shape of the vitreous chamber in the dynamics of the vitreous humour, as a first step and for the sake of simplicity, in this paper we will consider a domain with spherical shape and we will focus on the effect of the viscoelastic properties of the fluid on the characteristics and intensity of the steady streaming.
A similar problem has been considered by Nikolakis (Reference Nikolakis1990, Reference Nikolakis1992) who studied the motion of a viscoelastic liquid in a rigid spherical cavity attached to the end of a pendulum. He showed that this flow is equivalent to the flow of the same liquid in a torsionally oscillating sphere, and, assuming the amplitude of the oscillations to be small, found the leading-order flow and the second-order correction. He found that there were three qualitatively different flow patterns, including complete flow reversal, due to the competition between normal stress and inertial effects, and he found the regions in the parameter space wherein these three flow patterns exist.
Our paper extends the work of Nikolakis (Reference Nikolakis1990, Reference Nikolakis1992) in that we find two further qualitatively different types of flow, and we extend the parameter space study accordingly. We also study the dependence of the streaming intensity on the controlling parameters and evaluate the stresses on the sphere wall, which might have some importance in the application that motivated the present work. We finally derive analytically the solution for the steady streaming flow in the limit of small values of the Womersley number.
Owing to the lack of reliable rheological data, it is hard to infer information directly applicable to the realistic situation in the eye from the present model, especially in ageing subjects who have inhomogeneous properties in their vitreous humour. However, this work is intended as a useful conceptual step towards the understanding of transport processes in young healthy eyes.
2. Mathematical formulation
We develop a mathematical model of the dynamics of vitreous humour in an eye performing repeated saccadic movements. Following Astarita & Marrucci (Reference Astarita and Marrucci1974), we make the following assumptions about the material properties of the viscoelastic fluid representing the vitreous humour: the stress at a given point is determined by the history of deformation in the neighbourhood of that material point; the fluid has no preferred configuration (that is all possible configurations of material points are equivalent in the sense that differences in stress are due only to differences in the history of deformation); and the fluid has a fading memory. In addition, we assume that the history of deformation is entirely described by the deformation gradient,
$\unicode[STIX]{x1D641}$
, and that the fluid density is constant.
Under these conditions, and further assuming that the material deformations are small (which will be true for the small oscillatory rotations that are considered herein), it can be shown that the constitutive equation may be approximated at second-order by
$$\begin{eqnarray}\unicode[STIX]{x1D669}(\boldsymbol{r},t)=-\int _{0}^{\infty }G(s)\frac{\partial \unicode[STIX]{x1D63E}(\boldsymbol{r},s,t)}{\partial s}\text{d}s+\int _{0}^{\infty }\int _{0}^{\infty }m(s,s^{\prime })\frac{\partial \unicode[STIX]{x1D63E}(\boldsymbol{r},s,t)}{\partial s}\frac{\partial \unicode[STIX]{x1D63E}(\boldsymbol{r},s^{\prime },t)}{\partial s^{\prime }}\text{d}s\text{d}s^{\prime },\end{eqnarray}$$
where
$\unicode[STIX]{x1D669}$
is the extra-stress tensor (meaning that the total stress tensor equals
$-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D669}$
, where
$\unicode[STIX]{x1D644}$
is the identity matrix),
$\boldsymbol{r}$
identifies a position in space,
$t$
is time,
$\unicode[STIX]{x1D63E}$
is the right Cauchy–Green tensor,
$G(s)$
is the linear relaxation function and
$m(s,s^{\prime })$
is the second-order relaxation function (Astarita & Marrucci Reference Astarita and Marrucci1974; Böhme Reference Böhme1992). The functions
$G$
and
$m$
provide a Lagrangian measure of the memory of the fluid, and they both need to tend to zero sufficiently rapidly when
$s,s^{\prime }\rightarrow \infty$
for the assumption of fading memory to hold. The right Cauchy–Green tensor that appears in the above expression can be computed from the deformation gradient tensor
$\unicode[STIX]{x1D641}$
as
$\unicode[STIX]{x1D63E}=\unicode[STIX]{x1D641}^{\text{T}}\unicode[STIX]{x1D641}$
, and
$\unicode[STIX]{x1D641}$
is defined as
$\unicode[STIX]{x1D641}=\partial \hat{\boldsymbol{r}}/\boldsymbol{r}$
, where
$\hat{\boldsymbol{r}}(\boldsymbol{r},s,t)$
denotes the position of a material particle at time
$t-s$
, which at time
$t$
is at
$\boldsymbol{r}$
. The Eulerian velocity, which is the primitive kinematic variable used in the following, is related to
$\hat{\boldsymbol{r}}$
by
$$\begin{eqnarray}\boldsymbol{u}(\hat{\boldsymbol{r}},t-s)=-\frac{\partial \hat{\boldsymbol{r}}(\boldsymbol{r},s,t)}{\partial s}.\end{eqnarray}$$
Following Repetto et al. (Reference Repetto, Siggers and Stocchino2008), we model the vitreous chamber as a sphere of radius
$R$
performing small-amplitude, sinusoidal, torsional oscillations of prescribed angular displacement
${\it\beta}(t)$
, where
$$\begin{eqnarray}{\it\beta}(t)={\it\epsilon}\cos {\it\omega}t.\end{eqnarray}$$
We expand the velocity field
$\boldsymbol{u}$
and the pressure field
$p$
in terms of
${\it\epsilon}$
up to second order:
$$\begin{eqnarray}\displaystyle \boldsymbol{u}(\boldsymbol{r},t) & = & \displaystyle {\it\epsilon}\boldsymbol{u}_{1}(\boldsymbol{r},t)+{\it\epsilon}^{2}\boldsymbol{u}_{2}(\boldsymbol{r},t)+O({\it\epsilon}^{3}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p(\boldsymbol{r},t) & = & \displaystyle {\it\epsilon}p_{1}(\boldsymbol{r},t)+{\it\epsilon}^{2}p_{2}(\boldsymbol{r},t)+O({\it\epsilon}^{3}).\end{eqnarray}$$
The leading-order components
$\boldsymbol{u}_{1}$
,
$p_{1}$
were derived analytically by Meskauskas et al. (Reference Meskauskas, Repetto and Siggers2011), and, since they oscillate harmonically with frequency
${\it\omega}$
, we can write
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{1}(\boldsymbol{r},t) & = & \displaystyle \text{Re}\{\tilde{\boldsymbol{u}}_{1}(\boldsymbol{r})\text{e}^{\text{i}{\it\omega}t}\},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{1}(\boldsymbol{r},t) & = & \displaystyle \text{Re}\{\tilde{p}_{1}(\boldsymbol{r})\text{e}^{\text{i}{\it\omega}t}\},\end{eqnarray}$$
$\text{Re}$
denotes the real part and the dimensionless expressions of the functions
$\tilde{\boldsymbol{u}}_{1}$
and
$\tilde{p}_{1}$
are given in (3.1).
At second order, and, as in the viscous case (Repetto et al. Reference Repetto, Siggers and Stocchino2008), the flow is composed of steady streaming components and components whose frequency is double that of the primary flow:
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{2}(\boldsymbol{r},t) & = & \displaystyle \boldsymbol{u}_{2}^{st}(\boldsymbol{r})+\text{Re}\{\tilde{\boldsymbol{u}}_{2}(\boldsymbol{r})\text{e}^{2\text{i}{\it\omega}t}\},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{2}(\boldsymbol{r},t) & = & \displaystyle p_{2}^{st}(\boldsymbol{r})+\text{Re}\{\tilde{p}_{2}(\boldsymbol{r})\text{e}^{2\text{i}{\it\omega}t}\}.\end{eqnarray}$$
$\boldsymbol{u}_{2}^{st}(\boldsymbol{r})$
and
$p_{2}^{st}(\boldsymbol{r})$
. Böhme (Reference Böhme1992) derived the governing equations 
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\eta}_{0}{\rm\nabla}^{2}\boldsymbol{u}_{2}^{st}-\boldsymbol{{\rm\nabla}}p_{2}^{st}=\frac{{\it\rho}}{2}\text{Re}\left\{\tilde{\unicode[STIX]{x1D647}}_{1}\bar{\tilde{\boldsymbol{ u}}}_{1}\right\}-\text{Re}\left\{\frac{{\it\eta}^{\ast }({\it\omega})-{\it\eta}_{0}}{2\text{i}{\it\omega}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D63D}+4{\it\kappa}({\it\omega})\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\tilde{\unicode[STIX]{x1D659}}_{1}\bar{\tilde{\unicode[STIX]{x1D659}}}_{1}\right)\!\right\}\!,\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{2}^{st}=0\qquad & \displaystyle\end{eqnarray}$$
${\it\eta}^{\ast }({\it\omega})={\it\eta}^{\prime }({\it\omega})-\text{i}{\it\eta}^{\prime \prime }({\it\omega})=\int _{0}^{\infty }G(s)\text{e}^{-\text{i}{\it\omega}s}\text{d}s$
, with
${\it\eta}_{0}={\it\eta}^{\ast }(0)$
. The velocity gradient is 
$$\begin{eqnarray}\unicode[STIX]{x1D647}(\boldsymbol{r},t)=\frac{\partial }{\partial \boldsymbol{r}}\boldsymbol{u}(\boldsymbol{r},t)={\it\epsilon}\unicode[STIX]{x1D647}_{1}(\boldsymbol{r},t)+{\it\epsilon}^{2}\unicode[STIX]{x1D647}_{2}(\boldsymbol{r},t)+O({\it\epsilon}^{3}),\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D647}_{1}(\boldsymbol{r},t)=\text{Re}\{\tilde{\unicode[STIX]{x1D647}}_{1}(\boldsymbol{r})\text{e}^{\text{i}{\it\omega}t}\},\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D63D}(\boldsymbol{r})=\frac{\partial }{\partial \boldsymbol{r}}\left(\!\tilde{\unicode[STIX]{x1D647}}_{1}\bar{\tilde{\boldsymbol{u}}}_{1}\right)+\left[\frac{\partial }{\partial \boldsymbol{r}}(\tilde{\unicode[STIX]{x1D647}}_{1}\bar{\tilde{\boldsymbol{u}}}_{1})\right]^{\text{T}}+\bar{\tilde{\unicode[STIX]{x1D647}}}_{1}^{\text{T}}\tilde{\unicode[STIX]{x1D647}}_{1}+\tilde{\unicode[STIX]{x1D647}}_{1}^{\text{T}}\bar{\tilde{\unicode[STIX]{x1D647}}}_{1}.\end{eqnarray}$$
The rate-of-deformation tensor is
$$\begin{eqnarray}\unicode[STIX]{x1D659}(\boldsymbol{r},t)={\textstyle \frac{1}{2}}\left[\unicode[STIX]{x1D647}(\boldsymbol{r},t)+\unicode[STIX]{x1D647}^{\text{T}}(\boldsymbol{r},t)\right]={\it\epsilon}\unicode[STIX]{x1D659}_{1}(\boldsymbol{r},t)+{\it\epsilon}^{2}\unicode[STIX]{x1D659}_{2}(\boldsymbol{r},t)+O({\it\epsilon}^{3})\end{eqnarray}$$
and
$\unicode[STIX]{x1D659}_{1}$
is expanded as
$$\begin{eqnarray}\unicode[STIX]{x1D659}_{1}(\boldsymbol{r},t)=\text{Re}\left\{\tilde{\unicode[STIX]{x1D659}}_{1}(\boldsymbol{r})\text{e}^{\text{i}{\it\omega}t}\right\}.\end{eqnarray}$$
Finally, the parameter
${\it\kappa}$
is given by
$$\begin{eqnarray}{\it\kappa}({\it\omega})=\frac{1}{2}\int _{0}^{\infty }\int _{0}^{\infty }m(s,s^{\prime })\cos {\it\omega}(s-s^{\prime })\text{d}s\text{d}s^{\prime }.\end{eqnarray}$$
Note that
${\it\kappa}$
is the only parameter that is required to be derived from the function
$m(s,s^{\prime })$
. As commented by Böhme, the parameters
${\it\eta}^{\prime }$
and
${\it\eta}^{\prime \prime }$
can be derived from the amplitude and phase of the oscillating shear stress, whilst
${\it\eta}_{0}$
can be found from the mean second-order shear stress and
${\it\kappa}$
from the normal stress. We note that the right-hand side of (2.7a
) can be found from the solutions
$\boldsymbol{u}_{1}$
and
$p_{1}$
, and the details of the calculation are provided in appendix A.
We now make the problem dimensionless by setting
$$\begin{eqnarray}\boldsymbol{u}^{\ast }=\frac{\boldsymbol{u}}{{\it\omega}R},\quad p^{\ast }=\frac{p}{{\it\eta}^{\prime }{\it\omega}},\quad r^{\ast }=\frac{r}{R},\quad t^{\ast }=t{\it\omega},\end{eqnarray}$$
where
$\text{}^{\ast }$
denotes dimensionless variables, and (2.7a
) and (2.7b
) become
$$\begin{eqnarray}\displaystyle {\rm\nabla}^{\ast 2}\boldsymbol{u}_{2}^{st\ast }-{\it\Gamma}\boldsymbol{{\rm\nabla}}^{\ast }p_{2}^{st\ast } & = & \displaystyle \frac{1}{2}{\it\alpha}^{2}{\it\Gamma}\text{Re}\{\tilde{\unicode[STIX]{x1D647}}_{1}^{\ast }\bar{\tilde{\boldsymbol{u}}}_{1}^{\ast }\}\nonumber\\ \displaystyle & & \displaystyle -\,\text{Re}\left\{\frac{{\it\Gamma}(1-\text{i}V)-1}{2\text{i}}\boldsymbol{{\rm\nabla}}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D63D}^{\ast }+4K{\it\Gamma}V\boldsymbol{{\rm\nabla}}^{\ast }\boldsymbol{\cdot }(\tilde{\unicode[STIX]{x1D659}}_{1}^{\ast }\bar{\tilde{\unicode[STIX]{x1D659}}}_{1}^{\ast })\right\},\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}^{\ast }\boldsymbol{\cdot }\boldsymbol{u}_{2}^{st\ast } & = & \displaystyle 0,\end{eqnarray}$$
$\boldsymbol{u}_{2}^{st\ast }=0$
at
$r^{\ast }=1$
. The problem is governed by four dimensionless parameters: 
$$\begin{eqnarray}{\it\alpha}=\sqrt{\frac{{\it\rho}{\it\omega}R^{2}}{{\it\eta}^{\prime }}},\quad V=\frac{{\it\eta}^{\prime \prime }}{{\it\eta}^{\prime }},\quad {\it\Gamma}=\frac{{\it\eta}^{\prime }}{{\it\eta}_{0}},\quad K=\frac{{\it\kappa}{\it\omega}}{{\it\eta}^{\prime \prime }}.\end{eqnarray}$$
The parameter
${\it\alpha}$
is the Womersley number, which is the only dimensionless parameter appearing in the analogous analysis in the case of a purely viscous fluid. The three additional dimensionless parameters are
$V$
, the ratio between the elastic and viscous components of the fluid, and can be interpreted as the inverse of the loss factor,
${\it\Gamma}$
, the ratio between the viscosity at frequency
${\it\omega}$
and the viscosity at zero frequency, and
$K$
, a parameter derived from the second-order relaxation function.
3. Solution
Much of the work in this section has been derived in a previous work by Nikolakis (Reference Nikolakis1990), but we present the complete derivation here so that the reader can follow the notation more easily. Hereinafter we work in terms of dimensionless variables but, for the sake of readability, we drop the stars on the variables. We work in spherical polar coordinates
$(r,{\it\theta},{\it\phi})$
in the radial, zenithal and azimuthal directions, respectively, with corresponding unit vectors
$\hat{\boldsymbol{r}}$
,
$\hat{{\bf\theta}}$
and
$\hat{{\bf\phi}}$
and velocity components
$(u_{r},u_{{\it\theta}},u_{{\it\phi}})$
. The leading-order solution, found by Meskauskas et al. (Reference Meskauskas, Repetto and Siggers2011), is purely azimuthal:
$$\begin{eqnarray}\tilde{\boldsymbol{u}}_{1}(\boldsymbol{r})={\tilde{u}}_{1{\it\phi}}({\it\theta},r)\hat{{\bf\phi}},\quad {\tilde{u}}_{1{\it\phi}}=g(r)\sin {\it\theta},\ g(r)=-\frac{\text{iJ}_{3/2}(ar)}{\sqrt{r}\text{J}_{3/2}(a)},\ a=\frac{{\it\alpha}\text{e}^{-\text{i}{\rm\pi}/4}}{\sqrt{1-\text{i}V}},\end{eqnarray}$$
and
$\tilde{p}_{1}$
is a constant.
We let the right-hand side of (2.15a
) be denoted
$\boldsymbol{{\mathcal{F}}}=\left({\mathcal{F}}_{r},{\mathcal{F}}_{{\it\theta}},{\mathcal{F}}_{{\it\phi}}\right)^{\text{T}}$
, and from (3.1) and (A 3)–(A 10) we obtain
$$\begin{eqnarray}{\mathcal{F}}_{r}=f_{r}(r)\sin ^{2}{\it\theta},\quad {\mathcal{F}}_{{\it\theta}}=f_{{\it\theta}}(r)\sin {\it\theta}\cos {\it\theta},\quad {\mathcal{F}}_{{\it\phi}}=0,\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle f_{r} & = & \displaystyle -\frac{{\it\alpha}^{2}{\it\Gamma}g{\bar{g}}}{2r}+{\it\Gamma}V\left(g^{\prime \prime }\left({\bar{g}}^{\prime }-\frac{{\bar{g}}}{r}\right)+\text{c.c.}\right)\nonumber\\ \displaystyle & & \displaystyle -\,K{\it\Gamma}V\left(\left\{g^{\prime \prime }\left({\bar{g}}^{\prime }-\frac{{\bar{g}}}{r}\right)+\frac{g^{\prime }{\bar{g}}}{r^{2}}+\text{c.c.}\right\}-\frac{g^{\prime }{\bar{g}}^{\prime }}{r}-\frac{g{\bar{g}}}{r^{3}}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle f_{{\it\theta}} & = & \displaystyle -\frac{{\it\alpha}^{2}{\it\Gamma}}{2}\frac{g{\bar{g}}}{r}-K{\it\Gamma}V\left(\left\{\frac{g^{\prime }{\bar{g}}}{r^{2}}+\text{c.c.}\right\}-\frac{g^{\prime }{\bar{g}}^{\prime }}{r}-\frac{g{\bar{g}}}{r^{3}}\right),\end{eqnarray}$$
$\text{c.c.}$
denotes the complex conjugate.
Following Repetto et al. (Reference Repetto, Siggers and Stocchino2008), we expand the vector-valued functions
$\boldsymbol{{\mathcal{F}}}$
and
$\boldsymbol{u}_{2}^{st}$
in terms of vector spherical harmonics
$(\mathbb{P}_{mn},\mathbb{B}_{mn},\mathbb{C}_{mn})$
and
$p_{2}^{st}$
in terms of spherical harmonics
$Y_{mn}$
, which are defined and explained in detail by Quartapelle & Verri (Reference Quartapelle and Verri1995), Arfken & Weber (Reference Arfken and Weber2001) and Meskauskas et al. (Reference Meskauskas, Repetto and Siggers2011). Since the solution is axisymmetric, we only require the set of axisymmetric harmonics
$\mathbb{P}_{0n}$
,
$\mathbb{B}_{0n}$
,
$\mathbb{C}_{0n}$
and
$Y_{0n}$
:
$$\begin{eqnarray}\displaystyle \boldsymbol{{\mathcal{F}}}(r,{\it\theta}) & = & \displaystyle \mathop{\sum }_{n=0}^{\infty }{\mathcal{F}}_{P0n}(r)\mathbb{P}_{0n}({\it\theta},{\it\phi})+\mathop{\sum }_{n=1}^{\infty }{\mathcal{F}}_{B0n}(r)\mathbb{B}_{0n}({\it\theta},{\it\phi})+\mathop{\sum }_{n=1}^{\infty }{\mathcal{F}}_{C0n}(r)\mathbb{C}_{0n}({\it\theta},{\it\phi}),\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{2}^{st} & = & \displaystyle \mathop{\sum }_{n=0}^{\infty }u_{2,0n}^{st}(r)\mathbb{P}_{0n}({\it\theta},{\it\phi})+\mathop{\sum }_{n=1}^{\infty }v_{2,0n}^{st}(r)\mathbb{B}_{0n}({\it\theta},{\it\phi})+\mathop{\sum }_{n=1}^{\infty }w_{2,0n}^{st}(r)\mathbb{C}_{0n}({\it\theta},{\it\phi}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{2}^{st} & = & \displaystyle \mathop{\sum }_{n=0}^{\infty }p_{2,0n}^{st}(r)Y_{0n}.\end{eqnarray}$$
$\mathbb{P}_{0n}$
,
$\mathbb{B}_{0n}$
and
$\mathbb{C}_{0n}$
are everywhere unit vectors pointing in the radial, zenithal and azimuthal directions, respectively (note that this property does not hold for non-axisymmetric vector spherical harmonics).
We first expand the function
$\boldsymbol{{\mathcal{F}}}(r,{\it\theta})$
to find the functions
${\mathcal{F}}_{P0n}(r)\ (n\geqslant 0)$
,
${\mathcal{F}}_{B0n}(r)\ (n\geqslant 1)$
and
${\mathcal{F}}_{C0n}(r)\ (n\geqslant 1)$
, and, as in the case of a purely viscous fluid (Repetto et al.
Reference Repetto, Siggers and Stocchino2008), we find that all of these functions are zero, except for the following three:
$$\begin{eqnarray}{\mathcal{F}}_{P00}(r)=\frac{4}{3}\sqrt{{\rm\pi}}f_{r}(r),\quad {\mathcal{F}}_{P02}(r)=-\frac{4}{3}\sqrt{\frac{{\rm\pi}}{5}}f_{r}(r),\quad {\mathcal{F}}_{B02}(r)=-4\sqrt{\frac{{\rm\pi}}{30}}f_{{\it\theta}}(r).\end{eqnarray}$$
The irrational factors that appear in the above expressions are due to normalising coefficients in the definitions of the vector spherical harmonics. Note that, in the case of a purely viscous fluid (
$V=0$
,
${\it\Gamma}=1$
and
$K=0$
), we recover the expressions found by Repetto et al. (Reference Repetto, Siggers and Stocchino2008) ((2.11) in that paper).
Hence all the functions
$u_{2,0n}^{st}$
,
$v_{2,0n}^{st}$
,
$w_{2,0n}^{st}$
and
$p_{2,0n}^{st}$
are zero, except for possibly
$u_{2,00}^{st}$
,
$u_{2,02}^{st}$
,
$v_{2,02}^{st}$
,
$p_{2,00}^{st}$
and
$p_{2,02}^{st}$
. We obtain the same governing equations that were found by Repetto et al. (Reference Repetto, Siggers and Stocchino2008) (printed as (2.13a–e) in that paper, but replacing the variables
$u_{20,0}$
,
$u_{20,2}$
,
$v_{20,2}$
,
$p_{20,0}$
and
$p_{20,2}$
by
$u_{2,00}^{st}$
,
$u_{2,02}^{st}$
,
$v_{2,02}^{st}$
,
$p_{2,00}^{st}$
and
$p_{2,02}^{st}$
, respectively, and replacing the right-hand sides by
${\it\alpha}^{2}{\mathcal{F}}_{P0}$
,
${\it\alpha}^{2}{\mathcal{F}}_{P2}$
,
${\it\alpha}^{2}{\mathcal{F}}_{B2}$
by
${\mathcal{F}}_{P00}$
,
${\mathcal{F}}_{P02}$
,
${\mathcal{F}}_{B02}$
, respectively).
Solving these equations, we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle u_{2,00}^{st}(r)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{2,02}^{st}(r)=\sqrt{\frac{{\rm\pi}}{5}}{\tilde{u}}_{2,02}^{st} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{where} & \displaystyle \nonumber\\ \displaystyle & \displaystyle {\tilde{u}}_{2,02}^{st}=c_{1}r+c_{2}r^{3}+rI_{1}(r)-\frac{I_{2}(r)}{r^{2}}-r^{3}I_{3}(r)+\frac{I_{4}(r)}{r^{4}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v_{2,02}^{st}(r)=\sqrt{\frac{{\rm\pi}}{30}}\tilde{v}_{2,02}^{st} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{where} & \displaystyle \nonumber\\ \displaystyle & \displaystyle \tilde{v}_{2,02}^{st}=3c_{1}r+5c_{2}r^{3}+3rI_{1}(r)-5r^{3}I_{3}(r)-\frac{2I_{4}(r)}{r^{4}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}c_{1}=-I_{1}(1)+{\textstyle \frac{5}{2}}I_{2}(1)-{\textstyle \frac{7}{2}}I_{4}(1),\quad c_{2}=-{\textstyle \frac{3}{2}}I_{2}(1)+I_{3}(1)+{\textstyle \frac{5}{2}}I_{4}(1),\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle I_{1}(r) & = & \displaystyle -\frac{4}{15}\int _{0}^{r}f_{r}(\tilde{r})\text{d}\tilde{r},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle I_{2}(r) & = & \displaystyle -\frac{2}{15}\int _{0}^{r}\tilde{r}^{3}\left(2f_{r}(\tilde{r})+3f_{{\it\theta}}(\tilde{r})\right)\text{d}\tilde{r},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle I_{3}(r) & = & \displaystyle -\frac{4}{35}\int _{0}^{r}\frac{1}{\tilde{r}^{2}}\left(f_{r}(\tilde{r})-f_{{\it\theta}}(\tilde{r})\right)\text{d}\tilde{r},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle I_{4}(r) & = & \displaystyle -\frac{2}{35}\int _{0}^{r}\tilde{r}^{5}\left(2f_{r}(\tilde{r})+5f_{{\it\theta}}(\tilde{r})\right)\text{d}\tilde{r}.\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{u}_{2}^{st}=({\textstyle \frac{1}{4}}{\tilde{u}}_{2,02}^{st}(3\cos ^{2}{\it\theta}-1),-{\textstyle \frac{1}{4}}\tilde{v}_{2,02}^{st}\cos {\it\theta}\sin {\it\theta},0),\end{eqnarray}$$
with
${\tilde{u}}_{2,02}^{st}$
and
$\tilde{v}_{2,02}^{st}$
given by (3.6c
) and (3.6e
), respectively. We can also define an associated streamfunction,
${\it\psi}_{02}^{st}$
, given by
$$\begin{eqnarray}{\it\psi}_{02}^{st}={\textstyle \frac{1}{4}}r^{2}{\tilde{u}}_{2,02}^{st}(r)\sin ^{2}{\it\theta}\cos {\it\theta},\end{eqnarray}$$
with
$$\begin{eqnarray}\boldsymbol{u}_{2}^{st}=\left(\frac{1}{r^{2}\sin {\it\theta}}\frac{\partial {\it\psi}_{02}^{st}}{\partial {\it\theta}},-\frac{1}{r\sin {\it\theta}}\frac{\partial {\it\psi}_{02}^{st}}{\partial r},0\right).\end{eqnarray}$$
Finally, we compute the leading-order contribution to the non-dimensional shear (
${\it\tau}_{2\Vert }^{st}$
) and normal (
${\it\tau}_{2\bot }^{st}$
) stress at the wall induced by the steady streaming flow (normalised with the same scale as the pressure,
${\it\eta}^{\prime }{\it\omega}$
):
$$\begin{eqnarray}\displaystyle {\it\tau}_{2\Vert }^{st} & = & \displaystyle {\displaystyle \frac{5}{4{\it\Gamma}}}\left(-3I_{2}(1)+7I_{4}(1)\right)\cos {\it\theta}\sin {\it\theta},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\tau}_{2\bot }^{st} & = & \displaystyle p_{0}+\frac{5}{8}\left(-5I_{2}(1)+7I_{4}(1)\right)(3\cos ^{2}{\it\theta}-1)\nonumber\\ \displaystyle & & \displaystyle +\,V\left.\left\{g^{\prime }\overline{g}^{\prime }-\left(\frac{g\overline{g}}{r}\right)^{\prime }-Kr^{2}\left(\frac{g}{r}\right)^{\prime }\left(\frac{\overline{g}}{r}\right)^{\prime }\right\}\right|_{r=1}\sin ^{2}{\it\theta},\end{eqnarray}$$
$p_{0}$
is an arbitrary constant. The maximum of the wall shear stress is thus predicted to be located invariably at
${\it\theta}={\rm\pi}/4$
, whereas the normal stress attains extreme values at
${\it\theta}=0$
and
${\it\theta}={\rm\pi}/2$
. Note that the fact that there is a factor
${\it\Gamma}^{-1}$
in
${\it\tau}_{2\Vert }^{st}$
and not in
${\it\tau}_{2\bot }^{st}$
might appear odd, and so we comment upon it here. It results from the fact that the shear stress is proportional to components of the steady part of
$\boldsymbol{u}_{2}$
and is therefore multiplied by the zero-shear viscosity
${\it\eta}_{0}$
in the expression for the shear stress (leading to the factor
$1/{\it\Gamma}$
upon non-dimensionalisation). On the other hand the contributions to the normal stress, apart from the pressure, derive firstly from the elastic behaviour of the leading-order flow
$\boldsymbol{u}_{1}$
, leading to a term proportional to
${\it\eta}^{\prime \prime }({\it\omega})$
(which, after non-dimensionalisation, produces the first two terms of the second line of (3.12b
)) and secondly from a term proportional to the second-order relaxation function
$m(s,s^{\prime })$
(which becomes the last term in (3.12b
)).
4. Results
In this section we investigate the dependence of the streaming flow (3.4b
) upon the values of the governing parameters
$V$
,
$K$
,
${\it\alpha}$
and
${\it\Gamma}$
. Since
${\it\Gamma}$
acts only as a multiplicative parameter, we set
${\it\Gamma}=1$
throughout. Furthermore, we focus attention on relatively small values of the Womersley number
${\it\alpha}$
. This is because the constitutive equation is based on the assumption of small strains, meaning that the model is unsuitable to describe the flow for large
${\it\alpha}$
, when a boundary layer develops at the wall, with consequent high strains there. Since we also assume the flow is both symmetric in the equatorial plane (
${\it\theta}={\rm\pi}/2$
) and axisymmetric (independent of
${\it\phi}$
), in the plots we always show the flow in the quadrant
${\it\phi}=0$
and in the range
$0\leqslant {\it\theta}\leqslant {\rm\pi}/2$
.
The flow in this quadrant consists of one or more circulations. The centres of circulation are given by maxima and minima of
${\it\psi}_{02}^{st}$
. Using (3.10), we have
$$\begin{eqnarray}\displaystyle \frac{\partial {\it\psi}_{02}^{st}}{\partial {\it\theta}} & = & \displaystyle \frac{1}{4}r^{2}{\tilde{u}}_{2,02}^{st}(r)\sin {\it\theta}(3\cos ^{2}{\it\theta}-1),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\partial {\it\psi}_{02}^{st}}{\partial r} & = & \displaystyle \frac{1}{4}r\tilde{v}_{2,02}^{st}(r)\sin ^{2}{\it\theta}\cos {\it\theta}.\end{eqnarray}$$
Thus centres of circulations are given by points
$(r^{\ast },{\it\theta}^{\ast })$
where
$v_{2,02}^{st}(r^{\ast })=0$
and
${\it\theta}^{\ast }=\arccos (\frac{1}{3})$
. The circulations occupy the ranges of
$r$
between the zeros of the function
$u_{2,02}^{st}$
; that is, there are
$n$
circulations if
$u_{2,02}^{st}$
has
$n$
zeros (counting
$r=1$
as a zero, but not
$r=0$
). If
$u_{2,02}^{st}$
has a single maximum or minimum between two neighbouring zeros then this corresponds to a simple circulation, whilst if there are multiple extrema, the circulation has multiple centres with saddle points between them (an example of this is shown in figure 5, which will be explained later).
We recall from Repetto et al. (Reference Repetto, Siggers and Stocchino2008) that, in the case of a viscous fluid, the steady streaming consists of a single toroidal circulation cell in the upper hemisphere and one in the lower hemisphere; their sense of rotation is clockwise (in
${\it\phi}=0$
and
$0\leqslant {\it\theta}\leqslant {\rm\pi}/2$
), meaning that particles close to the axis of rotation drift towards the poles. With the present model, and using parameter values representative of a viscous fluid, which are
$V=0$
,
$K=0$
, we recover these results for the corresponding value of
${\it\alpha}$
.
In the case of a viscoelastic fluid the situation is much more complex. In fact, for each of the values of
$K$
that we investigated, we found that the
${\it\alpha}{-}V$
plane is divided into five regions, in each of which the steady streaming has qualitatively different characteristics. We show the case
$K=1$
in figure 1.
Figure 1. Curves in the plane
${\it\alpha}{-}V$
separating regions in which the steady streaming structure has different characteristics. Region 1: one clockwise rotating cell. Regions 2 and
$2^{\prime }$
: two circulation cells with opposite sense of rotation (inner counter-clockwise and outer clockwise). Region 3: one counter-clockwise rotating cell. Region 4: two counter-clockwise rotating cells separated by a saddle point. Region 5: three circulation cells.
$K=1$
.
With
$K=1$
, and for sufficiently low values of the parameter
$V$
(the limit
$V\rightarrow 0$
corresponds to a purely viscous fluid), viscous effects play a dominant role as expected, and the steady streaming flow is qualitatively similar to that in a purely viscous fluid (we denote the region of parameter space in which the flow behaves qualitatively in this way as region 1). A typical example of the velocity field and streamlines in region 1 are shown in figure 2.
Figure 2. Example of the steady streaming flow in region 1 of parameter space (using
${\it\alpha}=5$
,
$V=0.1$
,
$K=1$
). (a) Velocity vectors and contour plot of the velocity magnitude; (b) streamlines.
For all values of
${\it\alpha}$
, as
$V$
increases, which corresponds to increasing the elastic effects, a transition occurs from region 1 to a new region, region 2, in which there is an additional anticlockwise circulation cell close to centre of the sphere, whilst the original cell is squeezed towards the wall of the domain. An example of the flow in region 2 is shown in figure 3.
Figure 3. Example of the steady streaming flow in region 2 (using
${\it\alpha}=5$
,
$V=0.25$
,
$K=1$
). (a) Velocity vectors and contour plot of the velocity magnitude, (b) streamlines. Here and in the following figures, in (b) light shading indicates anticlockwise rotation and dark shading indicates clockwise rotation.
Further increase of
$V$
leads to annihilation of the first circulation cell, leaving only the second cell. We call this region 3, and an example of the flow here is shown in figure 4. Thus the flow pattern in region 3 is, qualitatively, a reversal of the flow in the viscous-dominated case, region 1.
Figure 4. Example of the steady streaming flow in region 3 (using
${\it\alpha}=5$
,
$V=1$
,
$K=1$
). (a) Velocity vectors and contour plot of the velocity magnitude, (b) streamlines.
If
${\it\alpha}$
is sufficiently large, a further increase of
$V$
leads to the appearance of other complex structures in the steady streaming flow. We denote these regions 4, 5 and
$2^{\prime }$
, and their arrangement in parameter space can be seen in figure 1. In region 4, there are two anticlockwise circulation cells, which are separated by a saddle point, see the example in figure 5. The steady streaming flow is dominated by the flow near to the centre of the sphere and the axis of rotation, and the velocities in the outer circulation cell are very small. Region 5 has a third circulation cell, and the flow is shown in figure 6. As in region 4, the circulation cell closer to the axis of rotation is characterised by much higher velocities than the other two. Finally, in region
$2^{\prime }$
, the flow has two circulation cells, as in region 2; however, here the velocity in the inner cell is much higher than that in the outer one, see figure 7. We note that, in the cases we investigated, the maximum streaming velocity is always located along the axis of rotation.
Figure 5. Example of the steady streaming flow in region 4 (using
${\it\alpha}=12$
,
$V=1.4$
,
$K=1$
). (a) Velocity vectors and contour plot of the velocity magnitude, (b) streamlines.
Figure 6. Example of the steady streaming flow in region 5 (using
${\it\alpha}=12$
,
$V=1.6$
,
$K=1$
). (a) Velocity vectors and contour plot of the velocity magnitude, (b) streamlines.
Figure 7. Example of the steady streaming flow in region
$2^{\prime }$
(using
${\it\alpha}=12$
,
$V=2$
,
$K=1$
). (a) Velocity vectors and contour plot of the velocity magnitude, (b) streamlines.
The boundaries of the regions in the
${\it\alpha}{-}V$
plane show a strong dependence on the value of
$K$
. In figure 8 we show the cases
$K=0.8$
and
$K=1.2$
(previously we considered
$K=1$
), and it can be seen that the topological arrangement of the regions in
${\it\alpha}{-}V$
space remains the same, but their boundaries are moved significantly.
The intensity of the steady streaming, as represented by its maximum velocity, also varies significantly with the controlling parameters, as shown in figure 9. In this figure we plot the streaming intensity versus the Womersley number for different values of (a)
$V$
and (b)
$K$
. The most obvious feature appears to be a general increase of this intensity with increasing values of
${\it\alpha}$
and
$V$
. We can also see some changes corresponding to boundaries between neighbouring regions being crossed.
Figure 9. Maximum streaming velocity versus the Womersley number
${\it\alpha}$
for different values of
$V$
(a) and
$K$
(b). (a)
$K=1$
, (b)
$V=0.5$
.
There are several indications in various biological contexts that cells are sensitive to the temporally averaged shear stress, i.e. in this case, the stress associated with the steady streaming flow, for which we have derived an expression in (3.12a
). In figure 10 we plot contour lines of the maximum wall shear stress in the plane
${\it\alpha}{-}V$
, for
$K=1$
.
Figure 10. Contour lines of the wall shear stress for different values of the parameters
${\it\alpha}$
and
$V$
. The shear stress is in the meridional (
${\it\theta}$
) direction, and it is plotted on
${\it\theta}={\rm\pi}/4$
, because this is where its absolute value is maximised; positive values indicate it is directed towards the equator and negative values that it is towards the poles. The thick line corresponds to the value 0.
In the case of small values of
${\it\alpha}$
we may calculate the solution analytically. Using (3.1), we note that, in the limit of small
${\it\alpha}$
,
$$\begin{eqnarray}\displaystyle g & = & \displaystyle -\text{i}r\left(1+\frac{a^{2}}{10}(1-r^{2})+\frac{a^{4}}{1400}(1-r^{2})(9-5r^{2})\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{a^{6}}{378\,000}\left(133-243r^{2}+135r^{4}-25r^{4}\right)+O({\it\alpha}^{8})\right).\end{eqnarray}$$
Using (3.3a
) and (3.3b
), we can find corresponding series expansions of the forcing terms
$f^{r}$
and
$f^{{\it\theta}}$
. Hence we can find the coefficients of the vector spherical harmonics from (3.5), perform the integrals (3.8), and find the components of the velocity or streamfunction from (3.6c
)–(3.6e
), (3.10). After some algebra, we obtain
$$\begin{eqnarray}\displaystyle {\it\psi}_{02}^{st} & = & \displaystyle \frac{{\it\Gamma}}{2\,079\,000(1+V^{2})^{2}}r^{3}(1-r^{2})^{2}\left\{1155V(1+V^{2})(1-2K){\it\alpha}^{4}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\left[5(2+r^{2})+V^{2}(1-2K)(131-50r^{2})\right]{\it\alpha}^{6}+O({\it\alpha}^{8})\right\}\sin ^{2}{\it\theta}\cos {\it\theta}.\end{eqnarray}$$
If the coefficient of
${\it\alpha}^{4}$
is not very small, then (since
${\it\alpha}$
is small)
${\it\psi}_{02}^{st}$
is single-signed in
$0<{\it\theta}<{\rm\pi}/2$
. Small values of this coefficient occur if either
$V\ll 1$
or
$|K-\frac{1}{2}|\ll 1$
. If neither of these conditions hold then:
-
(i) if
$K<\frac{1}{2}$
there is only a clockwise rotation (along axis to poles, around the surface to the equator and radially inward through the midplane); -
(ii) if
$K>\frac{1}{2}$
there is only an anticlockwise rotation.If the coefficient of
${\it\alpha}^{4}$
is small: -
(iii) in the region
$|K-\frac{1}{2}|\ll 1$
, we have a zero of
${\it\psi}_{02}^{st}$
at
$r^{\ast 2}=231V(1+V^{2})$
$(2K-1)/{\it\alpha}^{2}-2$
, which is in the range
$[0,1]$
if
$\frac{1}{2}<K<\frac{1}{2}+{\it\alpha}^{2}/(154V(1+V^{2}))$
; -
(iv) in the region
$V\ll 1$
, we have a zero of
${\it\psi}_{02}^{st}$
at
$r^{\ast 2}=231V(2K-1)/{\it\alpha}^{2}-2$
, which is in the range
$[0,1]$
if
$0<V<{\it\alpha}^{2}/(77(2K-1))$
.
In the limit of small
${\it\alpha}$
, these results predict the boundary between region 1 and region 2 to be at
$K=\frac{1}{2}$
and that between region 2 and region 3 at
$K=\frac{1}{2}+{\it\alpha}^{2}/(154V(1+V^{2}))$
. We computed these boundaries numerically for the case
$V=1$
, which is shown in figure 11. Comparing the numerical and analytical results, we found excellent agreement for the location of the boundary between regions 2 and 3; however, the agreement was less satisfactory for the boundary between regions 1 and 2, because the velocities in the limit
${\it\alpha}\rightarrow 0$
and
$K\approx \frac{1}{2}$
are comparable to the numerical errors.
Figure 11. Curves in the
${\it\alpha}{-}K$
plane showing regions in which the steady streaming flow has different characteristics, for
$V=1$
. The labels correspond to the regions whose flow characteristics are described in § 4. The analytical results for small
${\it\alpha}$
predict that these curves lie at
$K=\frac{1}{2}$
and
$K=\frac{1}{2}+{\it\alpha}^{2}/(154V(1+V^{2}))$
.
5. Discussion and conclusions
We have studied the steady streaming flow of a viscoelastic fluid contained in a rigid sphere performing periodic torsional oscillations about an axis passing through its centre. The present work extends that by Repetto et al. (Reference Repetto, Siggers and Stocchino2008), in which the same set-up was studied, but using a Newtonian fluid. A semi-analytical solution of the problem has been found by expanding all variables in terms of powers of the amplitude of oscillations
${\it\epsilon}$
, assumed small. Moreover, the velocity and pressure fields have been expanded in terms of vector and scalar spherical harmonics, respectively. A fully analytical solution of the problem was found in the limit of small values of the Womersley number
${\it\alpha}$
.
We have shown that a steady streaming flow is generated, and that its structure can be significantly more complicated than in the case of a purely viscous fluid. In fact, we identified five different regions in parameter space in which the steady streaming flow assumes qualitatively different patterns (the three regions identified by Nikolakis Reference Nikolakis1990 and two additional ones). Large streaming velocities are found for large values of
$V$
. Moreover, we showed analytically that, in the limit of small
${\it\alpha}$
, the streaming intensity is proportional to
${\it\alpha}^{4}$
, whereas it is proportional to
${\it\alpha}^{6}$
in the case of a purely viscous fluid. Therefore, in viscoelastic fluids the steady streaming intensity is expected to be larger than in purely viscous fluids for sufficiently small frequencies.
The work was motivated by the need to understand mass transport processes in the vitreous chamber of the eye induced by eye rotations. The clinical motivation arises because drug delivery to the retina through direct injection into the vitreous body is a frequently used treatment, which it is often more effective than topical, oral or intravenous delivery. However, the effectiveness of the treatment depends on the amount of drug reaching the retina and its spatial and temporal dependence, which in turn depends on the convection induced by the steady streaming flow.
There are several studies of the rheology of the vitreous body in the literature; however, we still do not have a reliable characterisation, due to difficulties inherent in measuring the mechanical properties. In particular, as soon as the vitreous body is extracted from the eye it rapidly degrades, leading to changes in the mechanical properties. Moreover, standard rheological tests are complicated by the tendency of the vitreous humour to form a lubrication layer on solid surfaces. As a result of the above difficulties, the existing measurements of vitreous humour properties are very sparse and highly dependent on the measurement technique adopted. At present, none of the existing measurements in the literature allow us to estimate the second-order relaxation function
$m$
appearing in the constitutive equation (2.1), and as we have shown in this paper, this has a strong quantitative influence on the results through the dimensionless parameter
$K$
. We note, however, that, as discussed in § 2, it would be possible to design experiments to obtain an estimate of the parameter
$K$
.
For these reasons, the present model cannot yet be used to predict the behaviour of the vitreous humour in the vitreous chamber. However, the model is conceptually important, since it allows us to identify all possible steady streaming patterns, and predict the dependence of the flow on the values of the governing dimensionless parameters.
Various authors have shown that, for a purely viscous fluid filling the vitreous chamber of the eye, the Péclet number associated with diffusive and advective mass transport is much greater than one, implying that mass transport is primarily driven by advection rather than diffusion. The present results suggest that this is also the case for a viscoelastic fluid, which is a better representation of the rheology of the healthy vitreous humour. A comprehensive theoretical study of mass transport phenomena in the vitreous cavity for the case of viscoelastic fluids has not been done, but will be usefully informed by the results in the present paper. We note, however, that other mechanisms which have not been considered here might play a very important role, such as Taylor dispersion in particular, which is expected to produce a large effect on the solute transport in this problem.
A limitation of our approach to understand the steady streaming in the vitreous cavity is related to the assumption that the rotations of the sphere are modelled by a single frequency. Several Fourier modes would be needed for describing real eye movements, and, since the analysis in the present paper is nonlinear, interactions between different modes would significantly complicate the picture. Nevertheless, this work has a significant conceptual relevance, showing that large streaming velocities can be generated, and represents a building block for future work that could address the above important issue.
Finally, we note that one of the research frontiers in ophthalmology is the development of artificial viscoelastic vitreous substitutes. The present results can be used to inform the optimal rheology of such a fluid.
Appendix A. Mathematical expressions
In this section, we report the mathematical expressions in terms of spherical coordinates
$(r,{\it\theta},{\it\phi})$
of the quantities appearing on the right-hand side of (2.15a
).
The gradient of a vector
$\boldsymbol{u}$
and divergence of a tensor
$\unicode[STIX]{x1D64F}$
are given by
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{u}=\left(\begin{array}{@{}ccc@{}}\displaystyle \frac{\partial u_{r}}{\partial r} & \displaystyle \frac{1}{r}\frac{\partial u_{r}}{\partial {\it\theta}}-\frac{1}{r}u_{{\it\theta}} & \displaystyle \frac{1}{r\sin {\it\theta}}\frac{\partial u_{r}}{\partial {\it\varphi}}-\frac{1}{r}u_{{\it\varphi}}\\ \displaystyle \frac{\partial u_{{\it\theta}}}{\partial r} & \displaystyle \frac{1}{r}\frac{\partial u_{{\it\theta}}}{\partial {\it\theta}}+\frac{1}{r}u_{r} & \displaystyle \frac{1}{r\sin {\it\theta}}\frac{\partial u_{{\it\theta}}}{\partial {\it\varphi}}-\frac{\cot {\it\theta}}{r}u_{{\it\varphi}}\\ \displaystyle \frac{\partial u_{{\it\varphi}}}{\partial r} & \displaystyle \frac{1}{r}\frac{\partial u_{{\it\varphi}}}{\partial {\it\theta}} & \displaystyle \frac{1}{r\sin {\it\theta}}\frac{\partial u_{{\it\varphi}}}{\partial {\it\varphi}}+\frac{1}{r}u_{r}+\frac{\cot {\it\theta}}{r}u_{{\it\theta}}\end{array}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}=\left(\begin{array}{@{}c@{}}\displaystyle \frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}T_{rr}\right)+\frac{1}{r\sin {\it\theta}}\frac{\partial }{\partial {\it\theta}}\left(\sin {\it\theta}T_{{\it\theta}r}\right)+\frac{1}{r\sin {\it\theta}}\frac{\partial T_{{\it\phi}r}}{\partial {\it\phi}}-\frac{T_{{\it\theta}{\it\theta}}+T_{{\it\phi}{\it\phi}}}{r}\\ \displaystyle \frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}T_{r{\it\theta}}\right)+\frac{1}{r\sin {\it\theta}}\frac{\partial }{\partial {\it\theta}}\left(\sin {\it\theta}T_{{\it\theta}{\it\theta}}\right)+\frac{1}{r\sin {\it\theta}}\frac{\partial T_{{\it\phi}{\it\theta}}}{\partial {\it\phi}}+\frac{T_{{\it\theta}r}}{r}-\frac{\cot {\it\theta}T_{{\it\phi}{\it\phi}}}{r}\\ \displaystyle \frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}T_{r{\it\phi}}\right)+\frac{1}{r\sin {\it\theta}}\frac{\partial }{\partial {\it\theta}}\left(\sin {\it\theta}T_{{\it\theta}{\it\phi}}\right)+\frac{1}{r\sin {\it\theta}}\frac{\partial T_{{\it\phi}{\it\phi}}}{\partial {\it\phi}}+\frac{T_{{\it\phi}r}}{r}+\frac{\cot {\it\theta}T_{{\it\phi}{\it\theta}}}{r}\end{array}\right). & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Hence the leading-order velocity gradient
$\tilde{\unicode[STIX]{x1D647}}_{1}$
is given by
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D647}}_{1}(r,{\it\theta})=\left(\begin{array}{@{}ccc@{}}0 & 0 & \displaystyle -\frac{g\sin {\it\theta}}{r}\\ 0 & 0 & \displaystyle -\frac{g\cos {\it\theta}}{r}\\ g^{\prime }\sin {\it\theta} & \displaystyle \frac{g\cos {\it\theta}}{r} & 0\end{array}\right).\end{eqnarray}$$
Using these, we can find:
$$\begin{eqnarray}\displaystyle & \tilde{\unicode[STIX]{x1D647}}_{1}\bar{\tilde{\boldsymbol{u}}}_{1}=\left(\begin{array}{@{}c@{}}\displaystyle -\frac{g{\bar{g}}\sin ^{2}{\it\theta}}{r}\\ \displaystyle -\frac{g{\bar{g}}\sin {\it\theta}\cos {\it\theta}}{r}\\ 0\end{array}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}(\tilde{\unicode[STIX]{x1D647}}_{1}\bar{\tilde{\boldsymbol{u}}}_{1})=\left(\begin{array}{@{}ccc@{}}\displaystyle -\left(\frac{g{\bar{g}}}{r}\right)^{\prime }\sin ^{2}{\it\theta} & \displaystyle -\frac{g{\bar{g}}}{r^{2}}\sin {\it\theta}\cos {\it\theta} & 0\\ \displaystyle -\left(\frac{g{\bar{g}}}{r}\right)^{\prime }\sin {\it\theta}\cos {\it\theta} & \displaystyle -\frac{g{\bar{g}}}{r^{2}}\cos ^{2}{\it\theta} & 0\\ 0 & 0 & \displaystyle -\frac{g{\bar{g}}}{r^{2}}\end{array}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \bar{\tilde{\unicode[STIX]{x1D647}}}_{1}^{\text{T}}\tilde{\unicode[STIX]{x1D647}}_{1}=\left(\begin{array}{@{}ccc@{}}g^{\prime }{\bar{g}}^{\prime }\sin ^{2}{\it\theta} & \displaystyle \frac{g{\bar{g}}^{\prime }}{r}\sin {\it\theta}\cos {\it\theta} & 0\\ \displaystyle \frac{g^{\prime }{\bar{g}}}{r}\sin {\it\theta}\cos {\it\theta} & \displaystyle \frac{g{\bar{g}}}{r^{2}}\cos ^{2}{\it\theta} & 0\\ 0 & 0 & \displaystyle \frac{g{\bar{g}}}{r^{2}}\end{array}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \tilde{\unicode[STIX]{x1D647}}_{1}^{\text{T}}\bar{\tilde{\unicode[STIX]{x1D647}}}_{1}=\overline{\bar{\tilde{\unicode[STIX]{x1D647}}}_{1}^{\text{T}}\tilde{\unicode[STIX]{x1D647}}_{1}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D63D}=\left(\begin{array}{@{}ccc@{}}\displaystyle 2\left(g^{\prime }{\bar{g}}^{\prime }-\left(\frac{g{\bar{g}}}{r}\right)^{\prime }\right)\sin ^{2}{\it\theta} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0\end{array}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \tilde{\unicode[STIX]{x1D659}}_{1}=\left(\begin{array}{@{}ccc@{}}0 & 0 & \displaystyle \frac{r}{2}\left(\frac{g}{r}\right)^{\prime }\sin {\it\theta}\\ 0 & 0 & 0\\ \displaystyle \frac{r}{2}\left(\frac{g}{r}\right)^{\prime }\sin {\it\theta} & 0 & 0\end{array}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \tilde{\unicode[STIX]{x1D659}}_{1}\bar{\tilde{\unicode[STIX]{x1D659}}}_{1}=\left(\begin{array}{@{}ccc@{}}\displaystyle \frac{r^{2}}{4}\left(\frac{g}{r}\right)^{\prime }\left(\frac{{\bar{g}}}{r}\right)^{\prime }\sin ^{2}{\it\theta} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & \displaystyle \frac{r^{2}}{4}\left(\frac{g}{r}\right)^{\prime }\left(\frac{{\bar{g}}}{r}\right)^{\prime }\sin ^{2}{\it\theta}\end{array}\right). & \displaystyle\end{eqnarray}$$
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D63D}$
and
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\tilde{\unicode[STIX]{x1D659}}_{1}\bar{\tilde{\unicode[STIX]{x1D659}}}_{1})$
.
Open access
