University of Birmingham Linear Rayleigh-Bénard stability of a transversely isotropic fluid

Suspended ﬁbres signiﬁcantly alter ﬂuid rheology, as exhibited in for example solutions of DNA, RNA and synthetic biological nanoﬁbres. It is of interest to determine how this altered rheology aﬀects ﬂow stability. Motivated by the fact thermal gradients may occur in biomolecular analytic devices, and recent stability results, we examine the problem of Rayleigh–B´enard convection of the transversely isotropic ﬂuid of Ericksen. A transversely isotropic ﬂuid treats these suspensions as a continuum with an evolving preferred direction, through a modiﬁed stress tensor incorporating four viscosity-like parameters. We consider the linear stability of a stationary, passive, transversely isotropic ﬂuid contained between two parallel boundaries, with the lower boundary at a higher temperature than the upper. To determine the marginal stability curves the Chebyshev collocation method is applied, and we consider a range of initially uniform preferred directions, from horizontal to vertical, and three orders of magnitude in the viscosity-like anisotropic parameters. Determining the critical wave and Rayleigh numbers, we ﬁnd that transversely isotropic eﬀects delay the onset of instability; this eﬀect is felt most strongly through the incorporation of the anisotropic shear viscosity, although the anisotropic extensional viscosity also contributes. Our analysis conﬁrms the importance of anisotropic rheology in the setting of convection.


Introduction
Suspended fibres significantly alter the rheology of the fluid, as exhibited in for example suspensions of DNA (Marrington et al., 2005), fibrous proteins of the cytoskeleton (Dafforn et al., 2004;Kruse et al., 2005), synthetic bio-nanofibres (McLachlan et al., 2013), extracellular matrix (Dyson et al., 2015) and plant cell walls (Dyson & Jensen, 2010).It is of interest to determine how this altered rheology affects flow stability; motivated by the impact of anisotropic effects on Taylor-Couette instability (Holloway et al., 2015), and the thermal gradients that may occur in devices which rely on nanofibre alignment for biomolecular analysis (Nordh et al., 1986), we examine the Rayleigh-Bénard instability of the transversely isotropic fluid of Ericksen.We consider the linear stability of a transversely isotropic fluid contained between two infinitely long horizontal boundaries of different temperatures (as shown in Figure 1), to , respectively.The leading order preferred direction is given by the angle θ (0) .a small arbitrary perturbation.Three different combinations of boundary types will be considered, (1) both boundaries are rigid, (2) both are stress free, and (3) the bottom boundary is rigid and the top is stress free.One application of our theory is to fibre-laden fluids, however it holds for any fluid which may be described as transversely isotropic.
In this paper, we adopt Ericksen's transversely isotropic fluid (Ericksen, 1960), which has been used to describe fibre-reinforced media (Cupples et al., 2017;Dyson et al., 2015;Green & Friedman, 2008;Holloway et al., 2015;Lee & Ockendon, 2005).Ericksen's model consists of mass and momentum conservation equations together with an evolution equation for the fibre director field.The stress tensor depends on the fibre orientation and linearly on the rate of strain; it takes the simplest form that satisfies the required invariances.Recently, Ericksen's model has been linked to suspensions of active particles (Holloway et al., 2018), such as self-propelling bacteria, algae and sperm (Saintillan & Shelley, 2013).Rayleigh (1916) was the first to form a mathematical model of the Rayleigh-Bénard system, using equations for the energy and state of an infinite layer of fluid, bounded by two stationary horizontal boundaries of different constant uniform temperatures.We work with the Boussinesq approximation that the flow is incompressible with non-constant density entering only through a buoyancy term.
In his original study, Rayleigh (1916) was able to find a closed-form solution in the case of both upper and lower boundaries being stress free, i.e. zero tangential stress; this setup has been simulated in experiments by replacing the bottom boundary with a layer of much less viscous fluid (leaving the top boundary stress free) (Goldstein & Graham, 1969).To determine the conditions where instability occurs for other combinations of boundary types, numerical techniques are required (Drazin, 2002).
We briefly discuss the equations and derive the steady state of the transversely isotropic model (Section 2), and then undertake a linear stability analysis, leading to an eigenvalue problem which is solved numerically (Sections 3 and 4).The effect of variations in viscositylike parameters and the steady state preferred direction on the marginal stability curves is considered (Section 5), then we conclude with a discussion of the results in Section 6.
Linear Rayleigh-Bénard stability of a TI fluid 3 and momentum balance leads to the generalised Navier-Stokes equations where ρ * 0 is the density at temperature T * 0 of the lower boundary, ρ * (x * , t * ) is the variable density of the fluid, t * is time, p * is the pressure, g * is acceleration due to gravity, ẑ is the unit vector in the z * -direction and τ * is the transversely isotropic stress tensor proposed by Ericksen (1960) (Dyson & Jensen, 2010;Holloway et al., 2018), μ * 1 implies the existence of a stress in the fluid even if it is instantaneously at rest, and can be interpreted as a tension or active behaviour in the fibre direction (Green & Friedman, 2008;Holloway et al., 2018), whilst the parameters μ * 2 and μ * 3 give the enhancement to the extensional and shear viscosities in the fibre direction, respectively, termed the anisotropic extensional and shear viscosities (Dyson & Jensen, 2010;Green & Friedman, 2008;Rogers, 1989).When the transversely isotropic model is interpreted as a suspension of fibres, the viscosity-like parameters may be calculated from the viscosity of the solvent fluid and the aspect ratio of the suspended fibres (Batchelor, 1970;Dyson & Jensen, 2010;Holloway et al., 2018).
We model the evolution of the preferred direction via the kinematic equation proposed by Green & Friedman (2008) which is a special case of the equation proposed by Ericksen (1960), appropriate for elongated particles with large aspect ratio.This equation takes into account the passive advection and reorientation of the fibres by the fluid flow.In the present study, we assume there is no active behaviour, i.e. μ * 1 = 0 (Holloway et al., 2018), therefore, the stress tensor is given by Temperature is governed by an advection-diffusion equation, where κ * is the coefficient of thermal conductivity (Chandrasekhar, 2013), and the constitutive relation for density is given as which is a linear function of temperature and independent of pressure (Drazin, 2002).
Here, α * is the coefficient of volume expansion, and we have assumed both quantities T * and ρ * are independent of the particles.We will consider two types of bounding surfaces; for both types of surface, we assume perfect conduction of heat and that the normal component of velocity is zero, i.e.
T * = T * 0 and w * = 0, at z * = 0, T * = T * 1 and w * = 0, at z * = d * . (2.8) The distinction between the types of bounding surfaces is then made through the final two boundary conditions.If the surface is rigid we impose no-slip boundary conditions, if the surface is stress free we impose zero-tangential stress, i.e.
u * = 0 on a rigid surface, τ * 12 = 0 on a stress-free surface. (2.9) Results will be presented from three groups of boundary conditions: both surfaces are rigid, both surfaces are stress free, and the bottom surface is rigid and the top surface is stress free.The evolution equation for the preferred direction (2.4) contains only a convective derivative for a and a term algebraic in a.Given that the velocity boundary conditions imply that there is zero inflow, it is sufficient to specify an initial condition for a only.

Steady state
Assuming that the parallel boundaries are infinitely long in the x-direction, a steady-state solution is given by where p 0 is some arbitrary pressure constant and the preferred fibre direction is described by the constant angle θ (0) to the x-axis (Figure 1).Convective flow of a suspension that has been uniformly aligned by a prior shear has previously been reported in the context of linear dichroism analysis of microtubules (Nordh et al., 1986).

Stability
We now examine the linear stability of the steady state described by equation (2.21), for the three different combinations of boundary types.We derive the first-order equations for an arbitrary perturbation, which are transformed into a generalised eigenvalue problem by assuming the solution takes the form of normal modes.

Linear stability analysis
We consider the stability of the steady-state solution to a perturbation where 0 < ε 1.As we have proposed a perturbation to the fibre orientation angle θ (0) , and not the alignment vector a directly, the form of a is given by Cupples et al. (2017) (3.5) Here, we have utilised the Taylor expansions for cos θ and sin θ.
Using the ansatz given in equations (3.1)-(3.5),we may state the following governing equations at first order.The incompressibility condition (2.11) becomes The first-order constitutive relations for stress (2.15) and fluid density (2.17) are given by τ (1) = P 2e (1) + μ 2 a (0) a (0) a (0) a (0) : e (1) + 2μ 3 a (0) a (0) • e (1) + e (1) • a (0) a (0) , (3.8) where e (1) = (∇ ∇ ∇u (1) + (∇ ∇ ∇u (1) ) T )/2 is the first-order rate-of-strain tensor.Notice that equations (3.6)-(3.9)are independent of the first-order alignment vector which is in turn governed by Finally, the equation governing temperature at next order is (3.12) The boundary conditions become homogeneous at first order, and are given by w (1) = u (1) = T (1) = 0 on a rigid surface, (3.13) w (1) = τ (1) 12 = T (1) = 0 on a stress-free surface. (3.14) After eliminating pressure and substituting for stress, the components of the momentum equation (3.7) are given by (3.16) Manipulating the components of the kinematic equation (2.12) yields an equation for the evolution of fibre direction (3.17) Notice equations (3.15)-(3.17)are decoupled, and so, we may solve the stability problem by considering only equations (3.12) and (3.16) with appropriate boundary conditions on w (1) and T (1) .The x-component of velocity and alignment angle may then be calculated from the solution for w (1) .We seek normal mode solutions to equations (3.12) and (3.16) of the form where k is the wave-number and s is the growth rate.Using this ansatz, equations (3.12) and (3.16) become ) and (3.20), i.e. for a given dimensionless wave-number k there will be non-trivial solutions (w , T ) to equations (3.19) and (3.20)only for certain values of s.We establish for each wave-number k the maximum Rayleigh number R l (k) such that the real part of all eigenvalues s are negative, i.e. the largest Rayleigh number such that the perturbation is stable and any disturbance decays to zero.The minimum of R l (k) is of particular interest, and is termed the critical Rayleigh number (R c ); it is used to determine the physical conditions under which instability first occurs (Acheson, 1990;Drazin, 2002;Koschmieder, 1993).If for a given experimental setup R < R c , then any perturbation decays exponentially to zero.The corresponding value of k at R c is also of interest; it describes the inverse wave-length of the convection currents and is termed the critical wave-number (k c ).

Numerical solution method
In order to determine the marginal stability curves R l (k), we must solve the eigenvalue problem (3.19) and (3.20) with boundary conditions given by (3.21).This is achieved using Chebyshev collocation, a spectral method that is capable of achieving high accuracy for low computational cost (Trefethen, 2000).
Using the Chebyshev differentiation matrix D (Trefethen, 2000) the linear operators on w and T in equations (3.19) and (3.20) may be approximated.This allows us to form the generalised matrix eigenvalue problem where s is the growth rate and eigenvalue of the problem, A and B are matrices, which are discrete representations of the linear operators which act on w and T , and the vector x contains the coefficients of the Lagrange polynomials, which approximate w and T at the Chebyshev points (equivalently the values of w and T at the Chebyshev points) (Trefethen, 2000).The matrices A and B may be constructed in MATLAB for each tuple of parameters θ (0) , μ 2 and μ 3 ; however, the matrices are not full rank as boundary conditions must be applied to close the problem.These constraints are applied using the method described by Hoepffner (2007); the solution space is reduced to consider only interpolants, which satisfy the boundary conditions.We may therefore compute the eigenvalue s for a range of parameters (θ (0) , μ 2 , μ 3 ) and Rayleigh number R using the inbuilt eigenvalue solver in MATLAB eig; this solver employs the QZ-algorithm for generalized eigenvalue problems.We then determine the Rayleigh number for which the eigenvalue is zero using the MATLAB function fzero, i.e. disturbances neither grow nor decay, and fminsearch to determine the critical wave and Rayleigh numbers.
As the Prandtl number (P) only appears in combination with the growth rate s, we do not consider variations in P as we are interested in the marginal stability curves where s = 0, i.e. the boundary between stability and instability.
To accommodate uncertainty in parameter values, we have performed an extensive parameter search for a wide range of steady state preferred directions (0 6 θ (0) 6 π/2) and viscosities (0 6 μ 2 , μ 3 6 1, 000).Notice that the solution is periodic in the steady state preferred direction with period π/2.
To validate our numerical procedure, we compared our results with those of Dominguez-Lerma et al. (1984) and Rayleigh (1916) for the Newtonian case, i.e. μ 2 = μ 3 = 0; we will denote the critical Rayleigh number for the Newtonian case R N .When both boundaries are stress free, our numerical approximation of the Rayleigh number is within 10 −12 of the known analytical result R N = 27π 4 /4.When both boundaries are rigid our numerical approximation of the Rayleigh number is within 10 −7 of the value of R N found by Dominguez-Lerma et al. (1984).

Results
In Section 5.1, we first determine the marginal stability curves R l (k); for any value of k, an experimental set-up satisfying R < R l (k) is stable for that wavelength, whereas if R lies above R l (k) the system is unstable.We calculate these curves for a range of https://www.cambridge.org/core/terms.https://doi.org/10.1017/S0956792518000359Downloaded from https://www.cambridge.org/core.University of Birmingham, on 06 Jul 2018 at 15:08:14, subject to the Cambridge Core terms of use, available at non-dimensional parameters representing the steady state preferred direction θ (0) , the anisotropic extensional viscosity μ 2 and the anisotropic shear viscosity μ 3 for different combinations of boundary conditions.We determine the critical wave and Rayleigh numbers for each tuple of non-dimensional parameters (θ (0) , μ 2 and μ 3 ) by finding the wave-number at which R l (k) is minimal.Provided that the Rayleigh number for a given experiment lies below this critical value, the system will be stable to small perturbations for all wavelengths and the fluid will be motionless.In Section 5.2, we will make an empirical approximation to the dependence of the critical wave and Rayleigh numbers on the problem parameters.

Critical wave and Rayleigh number
Figure 2 shows the critical wave-number (k c ) as a function of the steady state preferred direction (θ (0) ), for selected values of the anisotropic extensional (μ 2 ) and shear (μ 3 ) viscosities, with both boundaries rigid.The critical wave-number is related to the width of a convection cell; increases in k c reduce the width of the convection cell.Notice that Figures 2(a)-(d) are symmetric about θ (0) = π/4, where the maximum of k c is achieved.In Figure 2(a), we examine the effect of the anisotropic extensional viscosity with the anisotropic shear viscosity set to zero.The horizontal line corresponds to the Newtonian/isotropic case, and hence there is no dependence on the fibre direction θ (0) .As μ 2 is increased, the limiting form of the critical curve between π/8 6 θ (0) 6 3π/8 is quickly approached, with changes to μ 2 above 100 having only a small effect.In the ranges 0 6 θ (0) 6 π/8 and 3π/8 6 θ (0) 6 π/2, the changes to the critical wave-number occur much more slowly with respect to μ 2 , with a local minimum occurring for values of μ 2 above 250 around θ (0) = 0.1 and θ (0) = 1.5.The impact of changing the anisotropic extensional viscosity on the wave-number is therefore dependent on the steady state fibre direction.If the fibres are aligned near horizontal or vertical, the wave-number is decreased and the width of the convection cell increased; if the fibre direction is at π/4 to the horizontal at steady state, then the wave-number increases, and hence the width of the convection cell decreases.Observing how the critical curves change between Figures 2(a)-(d) allows us to identify the impact of the anisotropic shear viscosity μ 3 .As μ 3 is increased it dampens changes to the critical wave-number caused by changes in μ 2 , nearly removing the dependence on θ (0) completely in Figure 2(d), where μ 3 = 1, 000.
Figure 3 shows the critical Rayleigh number (R c ) as a function of θ (0) for changes in μ 2 and μ 3 with both boundaries rigid.Figure 3(a) shows the change in R c neglecting anisotropic shear viscosity, i.e. μ 3 = 0.The lowest horizontal line corresponds to the Newtonian case, and has no dependence on θ (0) as expected.For μ 2 .100 this horizontal line is simply translated to higher values of R c , with little to no dependence on θ (0) .As μ 2 is increased further the shape of the critical curves change dramatically.Global minima occur at θ (0) ≈ 0.3, 1.2, local maxima occur at θ (0) = 0, π/2, and the global maximum at θ (0) = π/4; the difference between the global minimum and maximum is approximately 6 × 10 4 for μ 2 = 1, 000.Therefore, when the anisotropic extensional viscosity is large and the anisotropic shear viscosity is negligible the steady state is most unstable for steady state fibre orientations close to π/16 of horizontal or vertical, however for smaller angles to the horizontal or vertical the stability sharply increases.when the steady state direction is at π/4 to the horizontal.Examining Figures 3(a)-(d) allows us to identify how the anisotropic shear viscosity affects the stability of the steady state.We observe that increasing μ 3 increases R c , hence making the steady state more stable.However, this relationship is not uniform for different values of θ (0) , as can be seen by noting that when μ 2 = 1, 000 and μ 3 = 0 the most stable value of θ (0) is π/4, but as μ 3 is increased to 1, 000 then θ (0) = π/4 becomes the most unstable value.Therefore, increasing the anisotropic shear viscosity has the most stabilising effect for values of the steady state fibre orientation near horizontal and vertical, and a slightly weaker effect when the steady state direction is π/4.However, increases in the anisotropic shear viscosity always stabilise the steady state for all choices of anisotropic extensional viscosity and steady state preferred directions.
(a) (b) Similar results are obtained for critical wave-number when both boundaries are stress free (Figure 4), but where the critical wave-number of a Newtonian fluid (k N ) is smaller.Figure 5 shows the dependence of R c on θ (0) for selected values of μ 2 and μ 3 with both boundaries stress free.In Figure 5(a), μ 3 = 0 and the Newtonian case is represented by the lowest horizontal line.As μ 2 is increased a global maximum occurs at θ (0) = 0 and π/2 and global minimum at θ (0) = π/4, where R c does not increase from the critical Rayleigh number for the Newtonian/isotropic case (R c ≈ R N ).Near horizontal or vertical fibre-orientation, increasing the anisotropic extensional viscosity increases the threshold at which instability occurs, but when the steady state preferred direction is π/4 there is little change to the stability threshold as anisotropic extensional viscosity is varied.Figures 5(a)-(d) show that as μ 3 is increased, R c increases regularly, smoothing out the points of inflection that occur for small values of μ 2 .Therefore, increasing μ 3 stabilises the steady state for all values of θ (0) and μ 2 .Changes in anisotropic shear viscosity affect the magnitude of the critical Rayleigh number much more than changes to the anisotropic extensional viscosity.
Figure 6 shows k c as a function of θ (0) for selected values of μ 2 and μ 3 when the lower boundary is rigid and the top stress free, and shows a more intricate dependence on the tuple of parameters (θ (0) , μ 2 , μ 3 ) than when upper and lower boundaries match.In Figure 6(a), the horizontal line corresponds to the Newtonian/isotropic case, and hence has no dependence upon θ (0) , as expected.As μ 2 is increased the critical curves become more complex, in the range 0 6 θ (0) 6 π/8 and 3π/8 6 θ (0) 6 π/2 similar behaviour is observed to when both boundaries are the same, with the appearance of a local maximum at  θ (0) = 0, π/2 and a global minimum for values of θ (0) = 0.1, 1.4.However, for θ (0) between π/8 and π/4 extra inflection points are introduced compared with the matching boundary cases, but this variation becomes small for values of μ 2 larger than 100.We again identify from Figures 6(a)-(d) that μ 3 dampens the change in the critical wave-number due to μ 2 , eventually removing the dependence on θ (0) (Figure 6(d)).
Figure 7 shows the dependence of R c on θ (0) for selected values of μ 2 and μ 3 , with the lower boundary rigid and the upper stress free.In Figure 7(a), we examine how μ 2 affects R c when μ 3 = 0; the Newtonian/isotropic case is shown by the lowest horizontal line.As μ 2 is increased, R c increases, however this increase is not uniform with respect to θ (0) .Global maxima occur at θ (0) = 0 and π/2, local maxima at θ (0) ≈ π/8 and 3π/8, local minima at θ (0) ≈ π/16 and 7π/16, and the global minimum at θ (0) = π/4.Therefore, increasing the anisotropic extensional viscosity increases the stability threshold most when, at steady state, the fibres are either horizontal or vertical, and least when they are directed at an angle π/4 radians.As μ 3 increases, R c is increased, with the additional local maxima and minima becoming less pronounced, and disappearing completely once  that when the top and bottom boundaries are the same, the curves for the critical wavenumber take the same form, but with lower critical wave-numbers for the stress freestress free boundaries than the rigid-rigid case.When the boundaries are mixed (Figure 6), and the anisotropic shear viscosity is negligible, additional inflection points occur between π/8 6 θ (0) 6 3π/8.However, variation between the critical curves is small for medium to large values of the anisotropic extensional viscosity, and all changes are dampened as the anisotropic shear viscosity is increased, similarly to the matching boundary case.Therefore, in all cases, the anisotropic extensional viscosity gives rise to variations in the critical wave-number with respect to the steady state preferred direction, which are dampened by increases in the anisotropic shear viscosity.
Comparing Figures 3, 5, 7 allows us to compare how the different boundary conditions affect the critical Rayleigh number.Similarly to the Newtonian/isotropic case the most stable pair of boundaries is rigid-rigid, with the most unstable being stress free-stress free.In all boundary pairs increasing either the anisotropic extensional or shear viscosities increases the critical Rayleigh number, however changes to the anisotropic shear viscosity affect the stability threshold much more than equivalent changes to the anisotropic extensional viscosity.
We notice in Figures 2-7 that the critical wave and Rayleigh numbers are the same for θ (0) and π/2 − θ (0) , i.e. the material has the same stability characteristics when the steady state preferred direction is horizontal or vertical, the dependence on this angle being symmetric about π/4.Examining equation (3.19), we see that the real part of the left-hand side depends on angle via terms that can be expressed as being proportional to cos 4θ (0) , which provides a mathematical explanation for this phenomenon.The physical consequence of this result is that a suspension of horizontal fibres has the same stability properties as a suspension of vertical fibres.

Empirical forms of critical curves
The critical wave and Rayleigh numbers depend on three parameters μ 2 , μ 3 , and θ (0) .Examining Figures 2-7, we notice that, for medium to large values of the anisotropic shear viscosity, the critical curves are continuous with no sharp extrema (i.e.we expect the rate of change of the critical values with θ (0) to be continuous also).We may therefore attempt to fit analytic functions to the numerical results, for k c and R c .Such fitted functions can be used to determine the critical values for a given set of parameter values without using costly simulations.The critical wave-number is fit by minimising the maximum absolute error between the function and the numerical results through the simplex search method of Lagarias et al. (1998) (fminsearch in MATLAB); the critical Rayleigh number R c via the trust region method of nonlinear least squares fitting (fit in MATLAB).
For μ 3 > 40, we fit the critical wave-number to the empirical form where k N is the critical wave-number for a Newtonian fluid and f 1 to f 5 are fitting parameters.The form of this function implies that for μ 3 μ 2 , and μ 3 1, the critical wave-number approximates its Newtonian value.For μ 2 = 0 the critical wave-number k c does not depend on θ (0) or μ 3 , and takes the same value as in the Newtonian case. Figure 8 shows a comparison between the sampled numerical results and fitted function, where the fitted parameters are given in Table 1; excellent qualitative and good quantitative agreement is found.
For μ 3 > 40 we fit the critical Rayleigh number to the empirical form The fit for (a) (μ3 = 50, θ (0) = 0) and (b) (μ3 = 1, 000, θ (0) = 0), (c) (μ3 = 50, θ (0) = π/4), (d) (μ3 = 1000, θ (0) = π/4) where black, blue and red correspond to rigid-rigid (RR), rigid-stress free (RF) and stress free-stress free (FF) boundary pairs and circles represent numerical results.We choose θ (0) = 0, π/4 so that we consider the maximum and minimum values of cos 4θ (0) .where g 1 to g 6 are fitting parameters.The form of this function implies that R c is dependent upon θ (0) for μ 2 μ 3 , however this dependence is dampened as μ 3 is increased.We also observe from the relative sizes of g 4 and g 5 , in Table 2, that the increase of R c with μ 3 is much greater than that with μ 2 .Figure 9 shows a comparison between the sampled numerical results and fitted function, with the fitted values found in Table 2; a reasonable quantitative agreement and a good qualitative agreement is found.Note that g 6 plays the role of R N , and numerically is extremely close to its known Newtonian values.

Conclusion
In this paper, we extended the work of Rayleigh (1916) to study the linear stability of a transversely isotropic viscous fluid, contained between two horizontal boundaries of different temperatures, which are either rigid or stress free.We used the stress tensor first proposed by Ericksen (1960), with μ 1 = 0 (equivalent to a passive fluid (Holloway et al., 2018)), and a kinematic equation for the fibre-director field to model a transversely isotropic fluid.Numerically, we presented results for a range of steady state, initially uniform, preferred fibre directions from horizontal to vertical; this is equivalent to the full range of directions as the governing equations have a period of π/2.
As found recently for the Taylor-Couette flow of a transversely isotropic fluid (Holloway et al., 2015), the anisotropic shear viscosity μ 3 , is much more important in determining the stability of the flow than the anisotropic extensional viscosity μ 2 .The influence of this pair of parameters upon the stability of the flow depends on the uniform steady state preferred direction θ (0) , as well as the boundary conditions.Similarly to a Newtonian fluid, the most stable pair of boundaries is rigid-rigid, for which the temperature difference between the two boundaries required to induce instability is the largest; the least stable boundary pair is stress free-stress free.
The rheological parameters μ 2 and μ 3 also have an impact on the critical wave-number k c , which describes the width of the convection cells.We find for steady state preferred directions near horizontal or vertical, the width of the convection cell increases with the anisotropic extensional viscosity, when compared to a Newtonian fluid, and decreases when the preferred direction makes an angle of π/4 with the horizontal.The anisotropic shear viscosity dampens any changes to the critical wave-number caused by increases in the anisotropic extensional viscosity.If μ 3 μ 2 and μ 3 1, then there is very little change to the critical wave-number, and hence convection cell size, with changes to the anisotropic extensional viscosity or steady state preferred direction.taken into account.More complex models of anisotropic rheology include nematic liquid crystals (Landau & Lifshitz, 1986), and Brownian suspension mechanics (Batchelor, 1970;Hinch & Leal, 1972).The transversely isotropic fluid is a rational limit of an aligned suspension (Holloway et al., 2018) and conclusions drawn from the simpler model are useful in guiding intuition regarding more complex and general constitutive laws.
Fluids which exhibit transversely isotropic rheology are commonly found in many industrial and biological applications, therefore it is necessary to gain a better understanding of the underlying mechanics governing the behaviour of these materials.As a classical fluid mechanics problem modified to incorporate anisotropic rheology, we hope the Rayleigh-Bénard stability analysis undertaken here will motivate research into this fascinating area.

Figure 1 .
Figure 1.A schematic diagram of the Rayleigh-Bénard setup.The lower and upper boundaries are located at z * = 0 and z * = d * at temperatures T * 0 and T * 1, respectively.The leading order preferred direction is given by the angle θ(0) .
.20) where we have adopted the convention D = d/ dz.Equations (3.19) and (3.20) form an eigenvalue problem, which must be solved subject to the boundary conditions (2.19) and (2.20) rewritten as w = Dw = T = 0 on a rigid surface, w = D 2 w = T = 0 on a stress-free surface, (3.21) (see Appendix A).The growth rate s represents an eigenvalue to equations (3.19
μ 3 & 100, as can be identified by comparing Figures7(a)-(d).Again, changes in μ 3 affect R c far more than similar changes to μ 2 .Therefore, increases in the anisotropic shear viscosity stabilise the steady state, with the most stabilisation occurring when the fibres are oriented horizontally or vertically at steady state, at least when θ (0) = π/4.Comparing Figures 2-7 allows us to examine the effect of the boundary conditions on the critical wave and critical Rayleigh numbers.Examining Figures2 and 4, we identify

Table 1 .
The parameter values from curve fitting for the critical wave-number given in equation (5.1) for the different combinations of boundary conditions Comparison of fitted curves with the numerical results for the critical wave-number kc.

Table 2 .
The parameter values from curve fitting for the critical Rayleigh number in equation (5.2) for the different combinations of boundary conditions