1. Introduction
It is well known that a purely oscillatory boundary condition typically induces a velocity field with a nonzero time average, due to nonlinear effects in the underlying equations of motion. The timeaverage flow is referred to as the steady streaming, and thorough reviews of this phenomenon for the case of Newtonian fluids are provided by Riley (1967, 2001). A general methodology for studying the steady streaming in a viscoelastic fluid, which is of interest in the present paper, is described by Nikolakis (1992). 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 1992, 1994; Nickerson et al.
2008; Swindle, Hamilton & Ravi 2008; SharifKashani et al.
2011). A recent review of the state of the art in modelling the flow of vitreous humour is provided by Siggers & Ethier (2012).
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. (2012) used magnetic resonance imaging (MRI) to measure vitreous velocity fields under controlled sinusoidal rotations of the eye bulb. Rossi et al. (2012) 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 (2008) studied, using both analytical and experimental methods, the steady streaming component of the flow in a viscous fluid within a sphere that is performing smallamplitude, 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 (2011) studied the motion of a viscoelastic fluid in the same scenario, considering smallamplitude oscillations and finding the leadingorder 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 2006; Repetto, Siggers & Stocchino 2010; Meskauskas, Repetto & Siggers 2012), or by considering a realistic geometry either numerically (Balachandran & Barocas 2011; Abouali et al.
2012; Modareszadeh et al.
2012) or experimentally (Stocchino, Repetto & Cafferata 2007; Stocchino, Repetto & Siggers 2010; Bonfiglio et al.
2013).
In this study, we extend the work of Repetto et al. (2008) to the case of a viscoelastic fluid filling the chamber, closely following the approach of Böhme (1992), 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 (1990, 1992) 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 leadingorder flow and the secondorder 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 (1990, 1992) 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 (1974), 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 secondorder by
where
$\unicode[STIX]{x1D669}$
is the extrastress 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 secondorder relaxation function (Astarita & Marrucci 1974; Böhme 1992). 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
$ts$
, 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
Following Repetto et al. (2008), we model the vitreous chamber as a sphere of radius
$R$
performing smallamplitude, sinusoidal, torsional oscillations of prescribed angular displacement
${\it\beta}(t)$
, where
We expand the velocity field
$\boldsymbol{u}$
and the pressure field
$p$
in terms of
${\it\epsilon}$
up to second order:
Using this perturbation approach it can be shown (Böhme
1992) that the constitutive equation (
2.1) can be expressed only in terms of the history of Eulerian quantities.
The leadingorder components
$\boldsymbol{u}_{1}$
,
$p_{1}$
were derived analytically by Meskauskas et al. (2011), and, since they oscillate harmonically with frequency
${\it\omega}$
, we can write
where
$\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.
2008), the flow is composed of steady streaming components and components whose frequency is double that of the primary flow:
The aim of this paper is to calculate the steady streaming components
$\boldsymbol{u}_{2}^{st}(\boldsymbol{r})$
and
$p_{2}^{st}(\boldsymbol{r})$
. Böhme (
1992) derived the governing equations
(equations (40) and (41) in that paper), where overbars denote complex conjugates, and the symbols appearing in (
2.7a
) are explained in the following. The complex viscosity is given by
${\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
with
and
The rateofdeformation tensor is
and
$\unicode[STIX]{x1D659}_{1}$
is expanded as
Finally, the parameter
${\it\kappa}$
is given by
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 secondorder shear stress and
${\it\kappa}$
from the normal stress. We note that the righthand 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
where
$\text{}^{\ast }$
denotes dimensionless variables, and (2.7a
) and (2.7b
) become
subject to
$\boldsymbol{u}_{2}^{st\ast }=0$
at
$r^{\ast }=1$
. The problem is governed by four dimensionless parameters:
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 secondorder relaxation function.
3. Solution
Much of the work in this section has been derived in a previous work by Nikolakis (1990), 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 leadingorder solution, found by Meskauskas et al. (2011), is purely azimuthal:
and
$\tilde{p}_{1}$
is a constant.
We let the righthand 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
where
where
$\text{c.c.}$
denotes the complex conjugate.
Following Repetto et al. (2008), we expand the vectorvalued 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 (1995), Arfken & Weber (2001) and Meskauskas et al. (2011). 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}$
:
The functions
$\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 nonaxisymmetric 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.
2008), we find that all of these functions are zero, except for the following three:
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. (2008) ((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. (2008) (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 righthand 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
and
and
Thus the leadingorder part of the steady streaming flow is given by
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
with
Finally, we compute the leadingorder contribution to the nondimensional 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}$
):
where
$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 zeroshear viscosity
${\it\eta}_{0}$
in the expression for the shear stress (leading to the factor
$1/{\it\Gamma}$
upon nondimensionalisation). On the other hand the contributions to the normal stress, apart from the pressure, derive firstly from the elastic behaviour of the leadingorder flow
$\boldsymbol{u}_{1}$
, leading to a term proportional to
${\it\eta}^{\prime \prime }({\it\omega})$
(which, after nondimensionalisation, produces the first two terms of the second line of (
3.12b
)) and secondly from a term proportional to the secondorder 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
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. (2008) 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.
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.
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.
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 viscousdominated case, region 1.
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.
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.
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$
.
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}$
,
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
If the coefficient of
${\it\alpha}^{4}$
is not very small, then (since
${\it\alpha}$
is small)
${\it\psi}_{02}^{st}$
is singlesigned 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})$
$(2K1)/{\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(2K1)/{\it\alpha}^{2}2$
, which is in the range
$[0,1]$
if
$0<V<{\it\alpha}^{2}/(77(2K1))$
.
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.
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. (2008), in which the same setup was studied, but using a Newtonian fluid. A semianalytical 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 1990 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 secondorder 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 righthand side of (2.15a
).
The gradient of a vector
$\boldsymbol{u}$
and divergence of a tensor
$\unicode[STIX]{x1D64F}$
are given by
Hence the leadingorder velocity gradient
$\tilde{\unicode[STIX]{x1D647}}_{1}$
is given by
Using these, we can find:
Hence we can use (
A 2) to find
$\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})$
.
References
Abouali, O., Modareszadeh, A., Ghaffariyeh, A. & Tu, J.
2012
Numerical simulation of the fluid dynamics in vitreous cavity due to saccadic eye movement. Med. Engng Phys.
34 (6), 681–692.
Arfken, G. B. & Weber, H. J.
2001
Mathematical Methods for Physicists, 5th edn.
IAP Harcourt Academic Press.
Astarita, G. & Marrucci, G.
1974
Principles of NonNewtonian Fluid Mechanics. McGrawHill.
Balachandran, R. K. & Barocas, V. H.
2011
Contribution of saccadic motion to intravitreal drug transport: theoretical analysis. Pharmaceut. Res.
28 (5), 1049–1064.
Böhme, G.
1992
On steady streaming in viscoelastic liquids. J. NonNewtonian Fluid Mech.
44, 149–170.
Bonfiglio, A., Repetto, R., Siggers, J. H. & Stocchino, A.
2013
Investigation of the motion of a viscous fluid in the vitreous cavity induced by eye rotations and implications for drug delivery. Phys. Med. Biol.
58 (6), 1969–1982.
Lee, B., Litt, M. & Buchsbaum, G.
1992
Rheology of the vitreous body. Part 1: viscoelasticity of human vitreous. Biorheology
29, 521–533.
Lee, B., Litt, M. & Buchsbaum, G.
1994
Rheology of the vitreous body. Part 2: viscoelasticity of bovine and porcine vitreous. Biorheology
31 (4), 327–338; PMID: 7981433.
Meskauskas, J., Repetto, R. & Siggers, J. H.
2011
Oscillatory motion of a viscoelastic fluid within a spherical cavity. J. Fluid Mech.
685, 1–22.
Meskauskas, J., Repetto, R. & Siggers, J. H.
2012
Shape change of the vitreous chamber influences retinal detachment and reattachment processes: is mechanical stress during eye rotations a factor?
Invest. Ophthalmol. Vis. Sci.
53 (10), 6271–6281, PMID: 22899755.
Modareszadeh, A., Abouali, O., Ghaffarieh, A. & Ahmadi, G.
2012
Saccade movements effect on the intravitreal drug delivery in vitreous substitutes: a numerical study. Biomech. Model. Mechanobiol.
12 (2), 281–290.
Nickerson, C. S., Park, J., Kornfield, J. A. & Karageozian, H.
2008
Rheological properties of the vitreous and the role of hyaluronic acid. J. Biomech.
41 (9), 1840–1846.
Nikolakis, D.
1990
Strömung einer viskoelastischen Flüssigkeit in einem Kugelpendel. Z. Angew. Math. Mech.
5, T366–T368.
Nikolakis, D.
1992
Eine Theorie für Stationäre Drifterscheinungen Viskoelastischer Flüssigkeiten, Series 7: Strömungstechnik, vol. 209. VDI.
Piccirelli, M., Bergamin, O., Landau, K., Boesiger, P. & Luechinger, R.
2012
Vitreous deformation during eye movement. NMR Biomed.
25 (1), 59–66; PMID: 21567512.
Quartapelle, L. & Verri, M.
1995
On the spectral solution of the threedimensional Navier–Stokes equations in spherical and cylindrical regions. Comput. Phys. Commun.
90, 1–43.
Repetto, R.
2006
An analytical model of the dynamics of the liquefied vitreous induced by saccadic eye movements. Meccanica
41, 101–117.
Repetto, R., Siggers, J. H. & Stocchino, A.
2008
Steady streaming within a periodically rotating sphere. J. Fluid Mech.
608, 71–80.
Repetto, R., Siggers, J. H. & Stocchino, A.
2010
Mathematical model of flow in the vitreous humor induced by saccadic eye rotations: effect of geometry. Biomech. Model. Mechanobiol.
9 (1), 65–76.
Riley, N.
1967
Oscillatory viscous flows. Review and extension. J. Inst. Maths Applics.
3, 419–434.
Riley, N.
2001
Steady streaming. Annu. Rev. Fluid Mech.
33, 43–65.
Rossi, T., Querzoli, G., Pasqualitto, G., Iossa, M., Placentino, L., Repetto, R., Stocchino, A. & Ripandelli, G.
2012
Ultrasound imaging velocimetry of the human vitreous. Exp. Eye Res.
99 (1), 98–104, PMID: 22516112.
SharifKashani, P., Hubschman, J., Sassoon, D. & Kavehpour, H. P.
2011
Rheology of the vitreous gel: effects of macromolecule organization on the viscoelastic properties. J. Biomech.
44 (3), 419–423, PMID: 21040921.
Siggers, J. H. & Ethier, C. R.
2012
Fluid mechanics of the eye. Annu. Rev. Fluid Mech.
44 (1), 347–372.
Stocchino, A., Repetto, R. & Cafferata, C.
2007
Eye rotation induced dynamics of a Newtonian fluid within the vitreous cavity: the effect of the chamber shape. Phys. Med. Biol.
52, 2021–2034.
Stocchino, A., Repetto, R. & Siggers, J. H.
2010
Mixing processes in the vitreous chamber induced by eye rotations. Phys. Med. Biol.
55 (2), 453–467.
Swindle, K., Hamilton, P. & Ravi, N.
2008
In situ formation of hydrogels as vitreous substitutes: viscoelastic comparison to porcine vitreous. J. Biomed. Mater. Res. A
87A (3), 656–665.