## Abstract

Fibre-laden fluids are found in a variety of situations, while Couette devices are used for flow spectroscopy of long biological molecules, such as DNA and proteins in suspension. The presence of these fibres can significantly alter the rheology of the fluid, and hence must be incorporated in any modelling undertaken. A transversely isotropic fluid treats these suspensions as a continuum with an evolving preferred direction, through a modified stress tensor incorporating four viscosity-like parameters. We consider the axisymmetric linear stability of a transversely isotropic viscous fluid, contained between two rotating co-axial cylinders, and determine the critical wave and Taylor numbers for varying gap width and inner cylinder velocity (assuming the outer cylinder is fixed). Through the inclusion of transversely isotropic effects, the onset of instability is delayed, increasing the range of stable operating regimes. This effect is felt most strongly through incorporation of the anisotropic shear viscosity, although the anisotropic extensional viscosity also contributes. The changes to the rheology induced by the presence of the fibres therefore significantly alter the dynamics of the system, and hence should not be neglected.

## 1. Introduction

Fibre-laden fluids are found in numerous industrial and biological contexts, including DNA polymers, fibrous proteins of the cytoskeleton and lipid environments [1], extracellular matrix [2] and plant cell walls [3]. The presence of fibres can significantly alter the rheology and thus their effects must be taken into account. We consider the linear stability of transversely isotropic fluids within a Couette device [4] (two rotating co-axial cylinders, as shown in figure 1), to an axisymmetric perturbation. This flow set-up is pertinent to flow linear dichroism (LD), a technique used for analysing long biological molecules in suspension; it is also a canonical example in traditional fluid mechanics for the stability analysis of viscous flows [5].

At low angular velocity, the flow of a Newtonian fluid is laminar and purely azimuthal. Taylor [6] examined the stability of this flow both experimentally and analytically, under the assumption of small gap width (relative to the outer cylinder's radius). As angular velocity is increased, the flow loses stability and gives way to axisymmetric Taylor vortices; these can be visualized as a stack of fluid doughnuts. This is the first in a sequence of instabilities leading to turbulence as angular velocity is increased further. Using linear stability analysis, Taylor established a condition on the angular velocities of the rotating cylinders and kinematic viscosity of the fluid. The assumption of small gap width has since been relaxed via numerical techniques (e.g. [7]).

LD spectroscopy is a technique for studying the structure of bio-molecules; this technique exploits the difference in absorption of horizontally and vertically polarized light by a horizontally aligned sample. Flow LD works with bio-molecules in suspension and produces alignment through an applied shear flow [8] such as Couette [9,10] or channel flow [1]. The method is particularly useful for particles with high aspect ratio, such as DNA [11], protein fibres [12–14] and photosynthetic pigment–protein complexes [15]; one current area of interest is exploiting flow LD for pathogen detection [16]. The requirement to maximize signal production leads to a trade-off however; large shear rates produce higher degrees of fibre alignment (and hence improved signal-to-noise ratio), but increasing the inner cylinder angular velocity too high leads to a fluid instability which destroys all signal.

Stability of Couette flow of a fibre-laden fluid is therefore of practical as well as theoretical interest. These suspensions can be modelled as a class of anisotropic fluids, for which there is a single preferred fibre direction, and all properties of the fluid are identical normal to this direction. Fluids which exhibit these properties are termed *transversely isotropic*.

To examine the behaviour of a transversely isotropic fluid, we must define an appropriate constitutive law. We use the constitutive equation derived by Ericksen [17], the simplest form of stress tensor that is linearly dependent upon strain rate and obeys the required invariances to have a single preferred direction. This model introduces four viscosity-like parameters describing the anisotropic behaviour of the material, along with a preferred direction that evolves according to a kinematic equation. This equation involves a parameter λ which can be related to the shape of the fibres. The four viscosity-like parameters intuitively correspond to an isotropic solvent viscosity, a tension due to the presence of fibres, and anisotropic extensional and shear viscosities. The range 0<λ<1 corresponds to finite-length particles which, as shown by Jeffery [18], exhibit time-periodic behaviour and hence no steady state is achievable when Brownian effects are neglected. Green & Friedman [2] derived a kinematic equation corresponding to λ=1; this limit of Jeffery's model for fibres which have a large aspect ratio [3] does exhibit steady-state alignment.

Leslie [19] extended Taylor's analysis to a transversely isotropic fluid, still assuming a small gap width, taking |λ|>1 and neglecting the tension parameter. Leslie's analysis focused upon the effect of λ, as opposed to the different viscosity-like parameters, on the onset of instability. Wan & Lin [20] expanded on Leslie's study to consider finite gap widths, and the effect of extensional anisotropic viscosity. However, they neglected the effects of anisotropy in the shear viscosity and assume that the fibres instantaneously align with the streamlines.

Motivated by the LD applications discussed above, in which we are considering very long and thin molecules within a Couette device with a non-small gap width (relative to the outer cylinder's radius), we will work with λ=1, finite gap width, and retain all four viscosity-like parameters. We consider the stability of an axisymmetric perturbation because, by analogy with the isotropic case, we suggest this will be the first instability to arise.

We briefly discuss the equations and derive the steady state of the transversely isotropic model (§2), and then undertake a linear stability analysis, leading to an eigenvalue problem which is solved numerically (§3). The effect of variations in viscosity-like parameters and gap width on the marginal stability curves is considered (§4), while we conclude with discussion of the results in §5.

## 2. Governing equations and steady-state solution

We adopt a cylindrical coordinate system with radial, azimuthal and axial components (*r**,*θ*,*z**), and velocity vector **u***, consisting of components (*u**,*v**,*w**). Here, and in the following, stars denotes dimensional variables and parameters. Mass conservation and momentum balance leads to the generalized Navier–Stokes equations
*ρ** is the constant density of the fluid, *t** is time and ** σ*** is a transversely isotropic stress tensor that must be specified.

### (a) Transversely isotropic stress tensor

Ericksen's stress tensor requires a single preferred direction of the fibres, given by the unit vector **a**, which may vary with both position and time [17]. This relation occurs naturally in statistically homogeneous suspensions of dilute rods [21] (see also [3] for further discussion) and is the most general symmetric form of the stress tensor which is linear in rate of strain and where the direction of the fibres is physically indistinguishable. By neglecting that fibres may rotate about the axis **a**, the stress tensor is given as
*p** is the pressure, **e***=(**∇*****u***+**∇*****u**^{*T})/2 is the rate-of-strain tensor, **I** is the identity tensor, and *μ**,*μ**_{2} and *μ**_{3} are viscosities, while *μ**_{1} is a stress.

By setting *μ**_{1}=*μ**_{2}=*μ**_{3}=0, the stress tensor for an incompressible isotropic fluid is recovered, with viscosity *μ**. We may therefore interpret *μ** as the viscosity for shear flow in the direction transverse to the fibres [22], which is equivalent to interpreting *μ** as the isotropic component of the solvent viscosity, modified by the volume fraction of fibres [3]. We observe that the third term in equation (2.3) is independent of the strain rate and only related to **a**. Therefore, *μ**_{2} and *μ**_{3} by considering two-dimensional deformations in the plane of the fibres. The extensional viscosity parallel to the fibre direction is given by *μ**_{2} contributes only to *μ**_{3} distinguishes *μ**_{2} is the anisotropic extensional viscosity, while *μ**_{3} is the anisotropic shear viscosity [2,3,22].

### (b) Kinematic equation for fibre evolution

A kinematic equation is required to explain how the preferred direction of the fibres evolves depending on their physical properties. Ericksen's evolution equation for **a** based on the appropriate invariances and linear dependence upon **∇*****u*** is
**w***=[**∇*****u***−(**∇*****u***)^{T}]/2, and λ is a constant which describes the properties of the fibres in suspension [17].

We consider a special case of this condition (corresponding to λ=1) that is applicable for fibres with large aspect ratio, a derivation of which can be found in Green & Friedman [2],
**a**|=1, the model encodes only the local alignment direction of fibres and not length. This result generalizes Jeffery's treatment of the motion of long ellipsoidal particles [18], which align asymptotically in the direction of the principal rate of strain [3].

### (c) Boundary conditions

We take the inner and outer cylinders to have radii *ω** (figure 1). In the remainder of this paper, we will look for an axisymmetric solution which automatically satisfies the required periodicity conditions.

### (d) Non-dimensionalization

We non-dimensionalize the model by scaling the independent and dependent variables via:

The incompressibility condition (2.1) and the kinematic equation (2.5) remain unchanged by this scaling,
**e**=(**∇****u**+**∇****u**^{T})/2. The momentum balance (2.2) becomes
*μ*_{1} is the ratio of the effects of tension in the fibres to the transverse shear viscosity, while *μ*_{2} and *μ*_{3} are the ratios of the extensional viscosity and shear viscosity in the fibre direction to the transverse shear viscosity, respectively [2].

Finally, the boundary conditions (2.6), in dimensionless form, are
**u**, *p*, and **a**, respectively, subject to constitutive law (2.11) and boundary conditions (2.13).

### (e) Steady state

Assuming axisymmetry and treating the cylinders as infinitely long, we look for a steady-state solution that only depends on *r*. Taking the fluid velocity and fibres to align purely in the azimuthal direction, we make the ansatz
*β*=2*η*^{2}*Ω*/(1−*η*^{2}), and *p*_{0} is some arbitrary pressure constant.

## 3. Stability

### (a) Linear stability analysis

We consider the stability of solution (2.15) to an axisymmetric perturbation, such that
*ϵ*≪1. As **a** is a unit vector

We propose the solution takes the form
*k* is the wavenumber and *s* is the growth rate. Using this ansatz, equations (3.10) and (3.11) are reduced to
*μ*_{1}) plays no role in determining the stability of the perturbation. As in general *s*≠0, equations (3.7)–(3.9) then become
*D*=d/d*r* and *μ*_{2}+4*μ*_{3} represents the enhancement to the anisotropic extensional viscosity versus that perpendicular to fibres and 1+*μ*_{3} represents the shear viscosity of the fluid in the fibre direction, enhanced by the presence of fibres. Equations (3.14) must be solved subject to the boundary conditions (2.13) rewritten as

Following the same convention as Dominguez-Lerma *et al.* [7], we present our results in terms of the Taylor number *T*, which may be related to *β* and the physical parameters via
*k*, there will be non-trivial solution *u*′, *v*′ to the problem (3.14) only for certain values of *T*. There will be some least eigenvalue *T*_{l}(*k*) corresponding to each particular wavenumber *k*, and it is the least of these values over all *k* that we seek. This minimum value is often termed the critical Taylor number (*T*_{c}) and is used to determine the physical conditions under which instability first occurs [5].

### (b) Numerical solution method

To determine the marginal stability curve, we set *s*=0 in equations (3.14) and calculate the minimum eigenvalue *β* of the resulting differential system, for a range of wavelengths *k*. We find the critical wavelength *k*_{c}, eigenvalue *β*_{c} and Taylor number *T*_{c} for variations in the non-dimensional parameters *η*, *μ*_{2} and *μ*_{3}. When *T* is greater than the critical value *T*_{c} (or equivalently *β* is greater than the critical value *β*_{c}), instability will occur. These critical values depend on the mechanical and geometric properties of the system, along with the angular velocity of the inner cylinder and determine whether a given experimental set-up will display instability. These instabilities occur in the form of fluid doughnuts, spaced with wavelength *k*_{c}.

We solve the boundary value problem (3.14) and (3.15) numerically using a shooting method [7,26,27]. Three independent solutions (*u*_{1},*v*_{1}), (*u*_{2},*v*_{2}) and (*u*_{3},*v*_{3}) to equations (3.14) are found, which satisfy the respective initial value problems (IVP1-3) at *r*=*η* (cf. [27]),
*r*=*η*, while the extended conditions in IVP1-3 are chosen to form a basis of solutions.

Hence, we may express
*C*_{1}, *C*_{2} and *C*_{3} are chosen to satisfy the boundary conditions (3.15) at *r*=1. The existence of non-zero solutions of the linear system (3.20) requires
*r*=1 [7,27]. To validate our numerical procedure, we compared our results with those of Dominguez *et al.* [7] for the Newtonian case, i.e. *μ*_{2}=*μ*_{3}=0 and found agreement of the critical values to within 0.1%.

To accommodate uncertainty in parameter values, we have performed an extensive parameter search for a wide range of gap ratios and viscosities from 0 to 1000 times the isotropic component of viscosity. Results were calculated in intervals of 10 for the viscosity parameters *μ*_{2} and *μ*_{3} (similar behaviour was observed for viscosities below 100 when calculated in intervals of 0.1, however the data are not presented here), and 0.025 for the gap ratio *η*; data are available from the University of Birmingham institutional repository (see ‘Data accessibility’ section).

## 4. Results

We first determine the marginal stability curves *T*_{l}(*k*); for any value of *k*, an experimental set-up satisfying *T*<*T*_{l}(*k*) is stable for that wavelength, whereas if *T* lies above *T*_{l}(*k*) the system is unstable. We calculate these curves for a range of non-dimensional parameters representing the gap ratio *η*, the anisotropic extensional viscosity *μ*_{2} and the anisotropic shear viscosity *μ*_{3}. The effect of tension in the fibre direction *μ*_{1} was not varied, because this parameter does not affect the stability of the perturbation at this order. We determine the critical wave and Taylor numbers for each tuple of non-dimensional parameters (*η*,*μ*_{2} and *μ*_{3}) by finding the wavenumber at which *T*_{l}(*k*) is minimal. Provided that the Taylor number for a given experiment lies below this critical value, the system will be stable for all wavelengths.

Increasing the anisotropic extensional viscosity (*μ*_{2}), while neglecting anisotropic shear viscosity (*μ*_{3}=0), leads to more concave marginal stability curves, along with higher Taylor numbers for the corresponding wavenumber (figure 2*a*), i.e. instability is suppressed. Raising the anisotropic shear viscosity also leads to more concave marginal stability curves and higher Taylor numbers for corresponding wavenumbers (figure 2*b*). Instability is also suppressed as *μ*_{3} increases, however for commensurate increases in the anisotropic extensional and shear viscosities, the marginal stability curves are substantially different. Increases in the shear viscosity increase convexity and Taylor numbers far more than increases in the extensional viscosity. Greater suppression of the instability is obtained for increases in *μ*_{3} compared with *μ*_{2}.

Figure 3 shows the critical wavelength *k*_{c} as a function of *μ*_{3} for selected values of *μ*_{2}, and as a function of *μ*_{2} for selected values of *μ*_{3}. The critical wavenumber *k*_{c} is a decreasing function of each parameter, i.e. the fluid doughnuts which appear have shorter wavelength, however increases in *μ*_{3} cause a larger reduction to *k*_{c} than corresponding increases in *μ*_{2}.

Similar results are plotted in figure 4 for *T*_{c} as a function of *μ*_{2} and *μ*_{3}. Increasing either *μ*_{2} or *μ*_{3} leads to a higher critical Taylor number (greater suppression of instability). However, again the dependence of the critical Taylor number on the anisotropic extensional viscosity *μ*_{2} is much weaker than that on the anisotropic shear viscosity *μ*_{3}.

Finally, we consider the effect of the gap width parameter *η*. Decreasing the gap between the cylinders (increasing *η*) leads to an increase in the critical wavenumber *k*_{c} (i.e. larger fluid doughnuts are formed). The change in the critical wavenumber due to *μ*_{2} and *μ*_{3} is greater when the gap between the cylinders is large (small *η*), than when the gap width is small (large *η*) as shown in figure 5*a*.

Rather than look at the effect of *η* directly on the critical Taylor number (which itself contains *η*), we instead consider the critical value of the angular velocity of the inner cylinder, *Ω*_{c}. Increasing *η* (decreasing the inter-cylinder gap) decreases the critical angular velocity *Ω*_{c}, until *η* approaches approximately 0.8, whereupon the critical angular velocity increases again figure 5*b*.

## 5. Discussion

In this paper, we extended the work of Taylor and others [6,7,26,27] to study the axisymmetric linear stability of a transversely isotropic viscous fluid, contained in a Couette device. Similarly to Leslie [19], we used the stress tensor and kinematic equation first proposed by Ericksen [17] to model a transversely isotropic fluid. However, we chose an alternative case of the kinematic equation which corresponded to fibres with a large aspect ratio, this was motivated by considering long slender fibres that align in shear flow. Numerically, we presented results for a range of gap ratios from 0.2 to 0.975; in typical applications, this ratio is approximately 0.82 [9].

The parameter *μ*_{3}, which describes the ratio of the shear viscosity in the fibre direction to the transverse shear viscosity, is much more important in determining the stability of the flow than *μ*_{2}, which describes the ratio of the extensional viscosity to the transverse shear viscosity. While our results agree with those of Wan & Lin [20] in the effect of *μ*_{2}, we find that *μ*_{3} has a more significant effect and therefore must be included. Although the tension in the fibres, *μ*_{1}, was retained within the calculations, it did not influence the final stability analysis.

For values of the gap ratio less than 0.8, i.e. for a large gap width (relative to the outer cylinder's radius) between the cylinders, the critical angular velocity of the inner cylinder is a decreasing function of the gap ratio. However, for larger values of the gap ratio, i.e. a narrow gap width, the critical angular velocity increases as viscous effects dominate over the inertial effects giving rise to instability; more energy is then needed for the perturbation to become unstable.

Understanding the stability of fibre-laden fluids in Couette devices is relevant to flow LD spectroscopy [9]. Maximizing the shear rate, while still remaining in the stable flow regime, produces the strongest fibre alignment and hence signal. By considering the rheological changes induced by the presence of fibres, we showed that the critical angular velocity of the inner cylinder at the onset of instability is much higher than previously thought. Exploiting this understanding may enable stronger signals to be recovered from the use of LD Couette devices across a range of fibre types.

It is instructive to consider how our results on critical angular velocity may relate to the strength of shear-induced alignment, and hence LD signal. To do this, we must consider the dependence of average shear rate on signal, and the production of signal as polarized light passes through the depth of fluid (figure 6). In figure 6*a*, we consider the average shear rate
*Ω*_{c}. Alignment is an increasing function of shear rate [1] and hence higher alignment should be obtained for narrower gap widths (relative to the cylinder ratio). This alignment function *b*, depending on the experimental set-up [1]. Conversely, a narrower gap width provides less sample to interact with incident light producing a weaker signal for a given alignment. Thus the interplay of these two effects must be taken into consideration, as shown in figure 6*c*,d. We see that in the saturating regime, an optimal gap width may be determined (observe the local maxima in figure 6d), this optimum is a feature of incorporating the mechanical properties of the fibres and the precise location will depend on the functional form of

In order to use this model in practice, it is necessary to estimate the rheological parameters of the fluid, the most crucial being *μ** and *et al.* [28]. Our linear stability analysis could then be used to determine the other relevant parameter * Ψ*(≪1) is the volume fraction of fibres and

In future work, we will consider the instability of a non-symmetric perturbation to the steady state, in addition to the dispersion of fibres about the average direction. We also aim to include the effect of entanglement between individual fibres and establish how this may change the stability of the flow.

Fibre-laden fluids are ubiquitous within industrial and biological contexts, therefore it is necessary to gain a better understanding of how these fibres modify the rheological properties and physical behaviours in a range of settings. As a canonical example from traditional fluid mechanics, we hope the Taylor–Couette stability analysis undertaken here will provide insight and spark further interest in this area of research.

## Data accessibility

Supporting data have been submitted to the University of Birmingham institutional repository epapers.bham.ac.uk/1955, item ID 1955.

## Authors' contributions

All authors designed the research and co-wrote the paper, C.R.H. carried out the mathematical analysis and carried out the numerical simulations. R.J.D. supervised the project.

## Conflict of interests

We have no competing interests.

## Funding statement

C.R.H. is supported by an Engineering and Physical Sciences Research Council (EPSRC) doctoral training award (EP/J500367/1). R.J.D. gratefully acknowledges the support of the University of Birmingham's System Science for Health initiative. Computations were performed using the University of Birmingham's BlueBEAR HPC service, which was purchased through HEFCE SRIF-3 funds. See http://www.bear.bham.ac.uk.

## Acknowledgements

The authors are grateful to Prof. Alison Rodger (University of Warwick) for inspiring us to study the problem.

- Received March 2, 2015.
- Accepted April 16, 2015.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.