Linear Rayleigh-Bénard stability of a transversely-isotropic fluid
UUnder consideration for publication in Euro. Jnl of Applied Mathematics Linear Rayleigh-B´enard stability of atransversely-isotropic fluid
C. R. HOLLOWAY , D. J. SMITH , and R. J. DYSON ∗ School of Mathematics, University of Birmingham, B15 2TT, U.K. Institute for Metabolism and Systems Research, University of Birmingham, B15 2TT, U.K. ∗ email : [email protected] ( Received 4 August 2017 ) Suspended fibres significantly alter fluid rheology, as exhibited in for example solutionsof DNA, RNA and synthetic biological nanofibres. It is of interest to determine how thisaltered rheology affects flow stability. Motivated by the fact thermal gradients may occurin biomolecular analytic devices, and recent stability results, we examine the problem ofRayleigh-B´enard convection of the transversely-isotropic fluid of Ericksen. 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 considerthe linear stability of a stationary, passive, transversely-isotropic fluid contained betweentwo 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. Determiningthe critical wave and Rayleigh numbers we find that transversely-isotropic effects delaythe onset of instability; this effect is felt most strongly through the incorporation of theanisotropic shear viscosity, although the anisotropic extensional viscosity also contributes.Our analysis confirms the importance of anisotropic rheology in the setting of convection.
Key Words:
Suspended fibres significantly alter the rheology of the fluid, as exhibited in for examplesuspensions of DNA (Marrington et al., 2005), fibrous proteins of the cytoskeleton (Daf-forn 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 isof interest to determine how this altered rheology affects flow stability; motivated by theimpact of anisotropic effects on Taylor-Couette instability (Holloway et al., 2015), andthe thermal gradients that may occur in devices which rely on nanofibre alignment forbiomolecular analysis (Nordh et al., 1986), we examine the Rayleigh-B´enard instabilityof the transversely-isotropic fluid of Ericksen.We consider the linear stability of a transversely-isotropic fluid contained between two a r X i v : . [ phy s i c s . f l u - dyn ] A ug Holloway, Smith and
Dyson T ∗ = T ∗ z ∗ = d ∗ T ∗ = T ∗ z ∗ = 0 θ (0) Figure 1. A schematic diagram of the Rayleigh-B´enard setup. The lower and upperboundaries are located at z ∗ = 0 and z ∗ = d ∗ at temperatures T ∗ and T ∗ . The leadingorder preferred direction is given by the angle θ (0) .infinitely-long horizontal boundaries of different temperatures (as shown in Figure 1), toa small arbitrary perturbation. Three different combinations of boundary types will beconsidered, (1) both boundaries are rigid, (2) both are free, and (3) the bottom boundaryis rigid and the top is free. One application of our theory is to fibre-laden fluids, howeverit holds for any fluid which may be described as transversely-isotropic.In this paper we adopt Ericksen’s transversely-isotropic fluid (Ericksen, 1960), whichhas 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’smodel consists of mass and momentum conservation equations together with an evolutionequation for the fibre director field. The stress tensor depends on the fibre orientationand linearly on the rate of strain; it takes the simplest form that satisfies the requiredinvariances. Recently Ericksen’s model has been linked to suspensions of active particles(Holloway et al., 2017), 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´enard sys-tem, using equations for the energy and state of an infinite layer of fluid, bounded by twostationary horizontal boundaries of different constant uniform temperatures. We workwith the Boussinesq approximation that the flow is incompressible with non-constantdensity entering only through a buoyancy term. Given we consider infinitesimal motionof a liquid the Boussinesq approximation is equally valid as for a Newtonian fluid.In his original study Rayleigh (1916) was able to find a closed-form solution in thecase of both upper and lower boundaries being free, i.e. zero tangential stress; this setuphas been simulated in experiments by replacing the bottom boundary with a layer ofmuch less viscous fluid (leaving the top boundary free) (Goldstein & Graham, 1969). Todetermine the conditions where instability occurs for other combinations of boundarytypes, numerical techniques are required (Drazin, 2002).We briefly discuss the equations and derive the steady state of the transversely-isotropicmodel (section 2), and then undertake a linear stability analysis, leading to an eigenvalueproblem which is solved numerically (sections 3-4). The effect of variations in viscosity-like parameters and the steady state preferred direction on the marginal stability curvesis considered (section 5), then we conclude with a discussion of the results in section 6. inear Rayleigh-B´enard stability of a TI fluid We adopt a two dimensional Cartesian coordinate system ( x ∗ , z ∗ ), and velocity vec-tor u ∗ = ( u ∗ , w ∗ ); stars denote dimensional variables and parameters. In formulatingour governing equations we make use of the Boussinesq approximation (Chandrasekhar,2013), treating the density as constant in all terms except bouyancy. Mass conservationand momentum balance leads to the generalised Navier-Stokes equations ∇ ∗ · u ∗ = , (2.1) ρ ∗ (cid:18) ∂ u ∗ ∂t ∗ + ( u ∗ · ∇ ∗ ) u ∗ (cid:19) = − ∇ ∗ p ∗ + ∇ ∗ · τ ∗ − ρ ∗ g ∗ ˆ z , (2.2)where ρ ∗ is the density at temperature T ∗ of the lower boundary, ρ ∗ ( x ∗ , t ∗ ) is the variabledensity of the fluid, t ∗ is time, p ∗ is the pressure, g ∗ is acceleration due to gravity, ˆ z is theunit vector in the z ∗ -direction, and τ ∗ is the transversely-isotropic stress tensor proposedby Ericksen (1960), τ ∗ = 2 µ ∗ e ∗ + µ ∗ a a + µ ∗ a a a a : e ∗ + 2 µ ∗ ( a a · e ∗ + e ∗ · a a ) . (2.3)Ericksen’s stress tensor incorporates the single preferred direction a ( x ∗ , t ∗ ), the rate-of-strain tensor e ∗ = ( ∇ ∗ u ∗ + ∇ ∗ u ∗ T ) and viscosity-like parameters µ ∗ , µ ∗ , µ ∗ and µ ∗ . The parameter µ ∗ is the isotropic components of the solvent viscosity, modified bythe volume fraction of the fibres (Dyson & Jensen, 2010; Holloway et al., 2017), µ ∗ implies the existence of a stress in the fluid even if it is instantaneously at rest, andcan be interpreted as a tension in the fibre direction (Green & Friedman, 2008), whilstthe parameters µ ∗ and µ ∗ may be interpreted as the anisotropic extensional and shearviscosities respectively (Dyson & Jensen, 2010; Green & Friedman, 2008; Rogers, 1989).We model the evolution of the fibre direction via the kinematic equation proposed byGreen & Friedman (2008) ∂ a ∂t ∗ + ( u ∗ · ∇ ∗ ) a + a a a : ∇ ∗ u ∗ = ( a · ∇ ∗ ) u ∗ , (2.4)which is a special case of the equation proposed by Ericksen (1960), appropriate for fibreswith large aspect ratio. In the present study, we assume there is no active behaviour, i.e. µ ∗ = 0 (Holloway et al., 2017), therefore the stress tensor is given by τ ∗ = 2 µ ∗ e ∗ + µ ∗ a a a a : e ∗ + 2 µ ∗ ( a a · e ∗ + e ∗ · a a ) . (2.5)Temperature is governed by an advection-diffusion equation, ∂T ∗ ∂t ∗ + ( u ∗ · ∇ ∗ ) T ∗ = κ ∗ ∇ ∗ T ∗ , (2.6)where κ ∗ is the coefficient of thermal conductivity (Chandrasekhar, 2013), and the con-stitutive relation for density is given as ρ ∗ = ρ ∗ (cid:0) − α ∗ ( T ∗ − T ∗ ) (cid:1) , (2.7)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 fibres. Holloway, Smith and
Dyson
We will consider two types of bounding surfaces; for both types of surface we assumeperfect conduction of heat and that the normal component of velocity is zero, i.e. T ∗ = T ∗ and w ∗ = 0 , at z ∗ = 0 ,T ∗ = T ∗ and w ∗ = 0 , at z ∗ = d ∗ . (2.8)The distinction between the types of bounding surfaces is then made through the finaltwo boundary conditions. If the surface is rigid we impose no-slip boundary conditions,if the surface is free we impose zero-tangential stress, i.e. u ∗ = 0 on a rigid surface ,∂u ∗ ∂x ∗ = 0 on a free surface . (2.9)Results will be presented from three groups of boundary conditions: both surfaces arerigid, both surfaces are free, and the bottom surface is rigid and the top surface is free. The model is non-dimensionalised by scaling the independent and dependent variablesvia: x ∗ = d ∗ x , t ∗ = d ∗ κ ∗ t, u ∗ = κ ∗ d ∗ u ,T ∗ = β ∗ d ∗ T, ( p ∗ , τ ∗ ) = ρ ∗ κ ∗ d ∗ ( p, τ ) , ρ ∗ = ρ ∗ ρ, (2.10)where variables without asterisks denote dimensionless quantities, and β ∗ is the verticaltemperature gradient, as chosen in Drazin (2002), i.e. β ∗ = ( T ∗ − T ∗ ) /d ∗ . The incom-pressibility condition (2.1) and the kinematic equation (2.4) remain unchanged by thisscaling, ∇ · u = 0 , (2.11) ∂ a ∂t + ( u · ∇ ) a + a a a : ∇ u = ( a · ∇ ) u . (2.12)The momentum balance (2.2) becomes ∂ u ∂t + ( u · ∇ ) u = − ∇ p + ∇ · τ − RPB ρ ˆ z , (2.13)where we have introduced the following dimensionless parameters B = α ∗ β ∗ d ∗ , R = α ∗ β ∗ d ∗ g ∗ ρ ∗ κµ ∗ , P = µ ∗ ρ ∗ κ ∗ . (2.14)The Rayleigh number R is a dimensionless parameter relating the stabilising effects ofmolecular diffusion of momentum to the destabilising effects of buoyancy (Drazin, 2002;Koschmieder, 1993; Sutton, 1950), and the Prandtl number P relates the diffusion ofmomentum to diffusion of thermal energy (Chandrasekhar, 2013). Non-dimensionalisingthe stress tensor (2.5) yields τ = P (cid:0) e + µ a a a a : e + 2 µ ( a a · e + e · a a ) (cid:1) , (2.15) inear Rayleigh-B´enard stability of a TI fluid e = ( ∇ u + ∇ u T ) /
2, and non-dimensional parameters µ = µ ∗ µ ∗ , µ = µ ∗ µ ∗ , (2.16)have been introduced. Here µ and µ are the ratios of the extensional viscosity andshear viscosity in the fibre direction to the transverse shear viscosity, respectively (Green& Friedman, 2008; Holloway et al., 2015).The constitutive equation (2.7) for variable density is non-dimensionalised to give ρ = 1 + B ( T − T ) , (2.17)where T = T ∗ /β ∗ d ∗ and equation (2.6), which governs the temperature distribution,becomes ∂ T∂t + ( u · ∇ ) T = ∇ T. (2.18)Finally, the boundary conditions (2.8) and (2.9), in dimensionless form, are T = T , and w = 0 , at z = 0 ,T = T , and w = 0 , at z = 1 , (2.19)where T = T ∗ /β ∗ d ∗ . The distinction between the type of surface remains unchanged, u = 0 on a rigid surface ,∂ u∂x = 0 on a free surface . (2.20)The model consists of four governing equations (2.11), (2.12), (2.13), (2.18) for u , a , p and T , respectively, subject to constitutive laws (2.15) and (2.17) with boundaryconditions (2.19) and (2.20). Assuming that the parallel boundaries are infinitely long in the x -direction, a steadystate solution is given by u (0) = 0 , p (0) = p − RPB (cid:32) z + B z (cid:33) ,T (0) = T − z, ρ (0) = 1 + B z,θ (0) = constant , a (0) = (cos θ (0) , sin θ (0) ) , (2.21)where p is some arbitrary pressure constant and the preferred fibre direction is describedby the constant angle θ (0) to the x -axis (Figure 1). We now examine the linear stability of the steady state described by equations (2.21), forthe three different combinations of boundary types. We derive the first-order equations for
Holloway, Smith and
Dyson an arbitrary perturbation, which are transformed into a generalised eigenvalue problemby assuming the solution takes the form of normal modes.
We consider the stability of the steady state solution to a perturbation, u ( x, z, t ) = ε u (1) ( x, z, t ) + O (cid:16) ε (cid:17) , (3.1) p ( r, z, t ) = p (0) + εp (1) ( x, z, t ) + O (cid:16) ε (cid:17) , (3.2) T ( x, z, t ) = T (0) + εT (1) ( x, z, t ) + O (cid:16) ε (cid:17) , (3.3) θ ( x, z, t ) = θ (0) + εθ (1) ( x, z, t ) + O (cid:16) ε (cid:17) . (3.4)where 0 < ε (cid:28)
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) a = (cid:16) cos θ (0) − εθ (1) sin θ (0) , sin θ (0) + εθ (1) cos θ (0) (cid:17) + O (cid:16) ε (cid:17) . (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 governingequations at first order. The incompressibility condition (2.11) becomes ∇ · u (1) = 0 , (3.6)with conservation of momentum (2.13) given by ∂ u (1) ∂t = − ∇ p (1) + ∇ · τ (1) − RPB ρ (1) ˆ z . (3.7)The first order constitutive relations for stress (2.15) and fluid density (2.17) are givenby τ (1) = P (cid:18) e (1) + µ a (0) a (0) a (0) a (0) : e (1) + 2 µ (cid:16) a (0) a (0) · e (1) + e (1) · a (0) a (0) (cid:17)(cid:19) , (3.8) ρ (1) = −B T (1) , (3.9)where e (1) = ( ∇ u (1) + ( ∇ u (1) ) T ) / a (1) = (cid:16) − θ (1) sin θ (0) , θ (1) cos θ (0) (cid:17) , (3.10)which is in turn governed by ∂ a (1) ∂t + a (0) a (0) a (0) : ∇ u (1) = (cid:16) a (0) · ∇ (cid:17) u (1) . (3.11)Finally, the equation governing temperature at next order is ∂ T (1) ∂t − w (1) = ∇ T (1) . (3.12) inear Rayleigh-B´enard stability of a TI fluid w (1) = u (1) = T (1) = 0 on a rigid surface , (3.13) w (1) = ∂ u (1) ∂x = T (1) = 0 on a free surface . (3.14)After eliminating pressure and substituting for stress, the components of the momentumequation (3.7) are given by1 P ∇ (cid:32) ∂ u (1) ∂t (cid:33) = (cid:32) µ sin θ (0) µ (cid:33) ∇ u (1) + µ sin 4 θ (0) (cid:32) ∂ u (1) ∂x∂z − ∂ u (1) ∂x ∂z (cid:33) + cos 4 θ (0) ∂ u (1) ∂x ∂z , (3.15)1 P ∇ (cid:32) ∂ w (1) ∂t (cid:33) = (cid:32) µ sin θ (0) µ (cid:33) ∇ w (1) + R ∂ T (1) ∂x + µ sin 4 θ (0) (cid:32) ∂ w (1) ∂x∂z − ∂ w (1) ∂x ∂z (cid:33) + cos 4 θ (0) ∂ w (1) ∂x ∂z . (3.16)Manipulating the components of the kinematic equation (2.12) yields an equation forthe evolution of fibre direction, ∂ θ (1) ∂t = cos θ (0) ∂ w (1) ∂x − sin θ (0) ∂ u (1) ∂z − sin 2 θ (0) ∂ u (1) ∂x , (3.17)Notice equations (3.15), (3.16) and (3.17) are decoupled, and so we may solve the stabilityproblem by considering only equations (3.12) and (3.16) with appropriate boundaryconditions on w (1) and T (1) . The x -component of velocity and alignment angle may thenbe calculated from the solution for w (1) .We propose the solution to equations (3.12) and (3.16) takes the form w (1) = w (cid:48) ( z ) e st + ikx , T (1) = T (cid:48) ( z ) e st + ikx , (3.18)where k is the wave-number and s is the growth rate. Using this ansatz, equations (3.12)and (3.16) become (cid:34) (cid:32) µ sin θ (0) µ (cid:33) (cid:16) D − k (cid:17) + µ (cid:32) i sin 4 θ (0) (cid:16) kD − k D (cid:17) − cos 4 θ (0) k D (cid:33)(cid:35) w (cid:48) − R k T (cid:48) = s P (cid:16) D − k (cid:17) w (cid:48) , (3.19) w (cid:48) + (cid:104) D − k (cid:105) T (cid:48) = sT (cid:48) , (3.20)where we have adopted the convention D = d / d z . Equations (3.19) and (3.20) form aneigenvalue problem which must be solved subject to the boundary conditions (2.19) and Holloway, Smith and
Dyson (2.20) rewritten as w (cid:48) = Dw (cid:48) = T (cid:48) = 0 on a rigid surface ,w (cid:48) = D w (cid:48) = T (cid:48) = 0 on a free surface . (3.21)The growth rate s represents an eigenvalue to equations (3.19) and (3.20), i.e. for agiven dimensionless wave-number k there will be non-trivial solutions ( w (cid:48) , T (cid:48) ) to equa-tions (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 arenegative, i.e. the largest Rayleigh number such that the perturbation is stable and anydisturbance decays to zero. The minimum of R l ( k ) is of particular interest, and is termedthe critical Rayleigh number ( R c ); it is used to determine the physical conditions underwhich instability first occurs (Acheson, 1990; Drazin, 2002; Koschmieder, 1993). If for agiven 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-lengthof the convection currents and is termed the critical wave-number ( k c ). In order to determine the marginal stability curves R l ( k ) we must solve the eigenvalueproblem (3.19) and (3.20) with boundary conditions given by (3.21). This is achievedusing Chebyshev collocation, a spectral method that is capable of achieving high accuracyfor low computational cost (Trefethen, 2000).Using the Chebyshev differentiation matrix D (Trefethen, 2000) the linear operatorson w (cid:48) and T (cid:48) in equations (3.19) and (3.20) may be approximated. This allows us to formthe generalised matrix eigenvalue problem A x = s B x , (4.1)where s is the growth rate and eigenvalue of the problem, A and B are matrices whichare discrete representations of the linear operators which act on w (cid:48) and T (cid:48) , and the vector x contains the coefficients of the Lagrange polynomials which approximate w (cid:48) and T (cid:48) at the Chebyshev points (equivalently the values of w (cid:48) and T (cid:48) at the Chebyshev points)(Trefethen, 2000). The matrices A and B may be constructed in MATLAB for eachtuple of parameters θ (0) , µ and µ ; however, the matrices are not full rank as boundaryconditions must be applied to close the problem. These constraints are applied usingthe method described by Hoepffner (2007); the solution space is reduced to consideronly interpolants which satisfy the boundary conditions. We may therefore computethe eigenvalue s for a range of parameters ( θ (0) , µ , µ ) and Rayleigh number R usingthe inbuilt eigenvalue solver in MATLAB eig ; this solver employs the QZ -algorithm forgeneralized eigenvalue problems. We then determine the Rayleigh number for which theeigenvalue is zero using the MATLAB function fzero , i.e. disturbances neither grow nordecay, 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 donot 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 inear Rayleigh-B´enard stability of a TI fluid (cid:54) θ (0) (cid:54) π/ (cid:54) µ , µ (cid:54) π/ i.e. µ = µ = 0; we willdenote the critical Rayleigh number for the Newtonian case R N . When both boundariesare free, our numerical approximation of the Rayleigh number is within 10 − of theknown analytical result R N = 27 π /
4. When both boundaries are rigid our numericalapproximation of the Rayleigh number is within 10 − of the value of R N found byDominguez-Lerma et al. (1984). 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, whereasif R lies above R l ( k ) the system is unstable. We calculate these curves for a range ofnon-dimensional parameters representing the steady state preferred direction θ (0) , theanisotropic extensional viscosity µ and the anisotropic shear viscosity µ for differ-ent combinations of boundary conditions. We determine the critical wave and Rayleighnumbers for each tuple of non-dimensional parameters ( θ (0) , µ and µ ) by finding thewave-number at which R l ( k ) is minimal. Provided that the Rayleigh number for a givenexperiment lies below this critical value, the system will be stable to small perturbationsfor all wavelengths and the fluid will be motionless. In section 5.2 we will make an em-pirical approximation to the dependence of the critical wave and Rayleigh numbers onthe problem parameters. Figure 2 shows the critical wave-number ( k c ) as a function of the steady state preferreddirection ( θ (0) ), for selected values of the anisotropic extensional ( µ ) and shear ( µ ) vis-cosities, with both boundaries rigid. The critical wave-number is related to the width of aconvection cell; increases in k c reduce the width of the convection cell. Notice that Figures2(a) – (d) are symmetric about θ (0) = π/
4, where the maximum of k c is achieved. In Fig-ure 2(a) we examine the effect of the anisotropic extensional viscosity with the anisotropicshear viscosity set to zero. The horizontal line corresponds to the Newtonian/isotropiccase, and hence there is no dependence on the fibre direction θ (0) . As µ is increased,the limiting form of the critical curve between π/ (cid:54) θ (0) (cid:54) π/ µ above 100 having only a small effect. In the ranges 0 (cid:54) θ (0) (cid:54) π/ π/ (cid:54) θ (0) (cid:54) π/
2, the changes to the critical wave-number occur much more slowlywith respect to µ , with a local minimum occurring for values of µ above 250 around θ (0) = 0 . θ (0) = 1 .
5. The impact of changing the anisotropic extensional viscosity onthe wave-number is therefore dependent on the steady state fibre direction. If the fibresare aligned near horizontal or vertical, the wave-number is decreased and the width of theconvection cell increased; if the fibre direction is at π/ Holloway, Smith and
Dyson (a) (b) θ (0) k c µ = 00 π/ π/ θ (0) k c µ = 100 π/ π/ (c) (d) θ (0) k c µ = 1000 π/ π/ θ (0) k c µ = 10000 π/ π/ Figure 2. Critical wave-number ( k c ) for changes in the anisotropic extensional viscosity( µ ), the anisotropic shear viscosity ( µ ), and the preferred direction in the fluid at steadystate ( θ (0) ) with both boundaries rigid. In each subfigure the arrows indicate increasing µ ( µ = 0 , , , , , µ = 0, (b) µ = 10, (c) µ = 100, and (d) µ = 1000.Observing how the critical curves change between Figures 2(a) – (d) allows us to identifythe impact of the anisotropic shear viscosity µ . As µ is increased it dampens changesto the critical wave-number caused by changes in µ , nearly removing the dependence on θ (0) completely in Figure 2(d) where µ = 1000. Similar results are obtained when bothboundaries are free (Figure 3), but where the critical wave-number of a Newtonian fluid( k N ) is smaller.Figure 4 shows k c as a function of θ (0) for selected values of µ and µ when the lowerboundary is rigid and the top free, and shows a more intricate dependence on the tuple ofparameters ( θ (0) , µ , µ ) than when upper and lower boundaries match. In Figure 4(a) thehorizontal line corresponds to the Newtonian/isotropic case, and hence has no dependenceupon θ (0) , as expected. As µ is increased the critical curves become more complex, inthe range 0 (cid:54) θ (0) (cid:54) π/ π/ (cid:54) θ (0) (cid:54) π/ θ (0) = 0 , π/ θ (0) = 0 . , .
4. However, for θ (0) between π/ π/ inear Rayleigh-B´enard stability of a TI fluid (a) (b) θ (0) k c µ = 00 π/ π/ θ (0) k c µ = 100 π/ π/ (c) (d) θ (0) k c µ = 1000 π/ π/ θ (0) k c µ = 10000 π/ π/ Figure 3. Critical wave-number ( k c ) for changes in the anisotropic extensional viscosity( µ ), the anisotropic shear viscosity ( µ ), and the preferred direction in the fluid at steadystate ( θ (0) ) with both boundaries free. In each subfigure the arrows indicate increasing µ ( µ = 0 , , , , , µ = 0, (b) µ = 10, (c) µ = 100, and (d) µ = 1000.variation becomes small for values of µ larger than 100. We again identify from Figures4(a) – (d) that µ dampens the change in the critical wave-number due to µ , eventuallyremoving the dependence on θ (0) (Figure 4(d)).Figure 5 shows the critical Rayleigh number ( R c ) as a function of θ (0) for changes in µ and µ wtih both boundaries rigid. Figure 5(a) shows the change in R c neglectinganisotropic shear viscosity, i.e. µ = 0. The lowest horizontal line corresponds to theNewtonian case, and has no dependence on θ (0) as expected. For µ (cid:46)
100 this horizontalline is simply translated to higher values of R c , with little to no dependence on θ (0) . As µ is increased further the shape of the critical curves change dramatically. Global minimaoccur at θ (0) ≈ . , .
2, local maxima occur at θ (0) = 0 , π/
2, and the global maximum at θ (0) = π/
4; the difference between the global minimum and maximum is approximately6 × for µ = 1000. Therefore when the anisotropic extensional viscosity is large andthe anisotropic shear viscosity is negligible the steady state is most unstable for steadystate fibre orientations close to π/
16 of horizontal or vertical, however for smaller anglesto the horizontal or vertical the stability sharply increases. The most stable case when2
Holloway, Smith and
Dyson (a) (b) θ (0) k c µ = 00 π/ π/ θ (0) k c µ = 100 π/ π/ (c) (d) θ (0) k c µ = 1000 π/ π/ θ (0) k c µ = 10000 π/ π/ Figure 4. Critical wave-number ( k c ) for changes in the anisotropic extensional viscosity( µ ), the anisotropic shear viscosity ( µ ), and the preferred direction in the fluid at steadystate ( θ (0) ) with the bottom boundary rigid and the top boundary free. In each subfigurethe arrows indicate increasing µ ( µ = 0 , , , , , µ = 0, (b) µ = 10, (c) µ = 100, and (d) µ = 1000.the steady state direction is at π/ µ 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 bynoting that when µ = 1000 and µ = 0 the most stable value of θ (0) is π/
4, but as µ isincreased to 1000 then θ (0) = π/ π/
4. However, increases in the anisotropic shear viscosity alwaysstabilise the steady state for all choices of anisotropic extensional viscosity and steadystate preferred directions.Figure 6 shows the dependence of R c on θ (0) for selected values of µ and µ withboth boundaries free. In Figure 6(a) µ = 0 and the Newtonian case is represented bythe lowest horizontal line. As µ is increased a global maximum occurs at θ (0) = 0 and π/ θ (0) = π/
4, where R c does not increase from the critical Rayleigh inear Rayleigh-B´enard stability of a TI fluid (a) (b) θ (0) R c µ = 00 π/ π/ . . . . . × θ (0) R c µ = 100 π/ π/ . . . . . × (c) (d) θ (0) R c µ = 1000 π/ π/ . . . . . . × θ (0) R c µ = 10000 π/ π/ × Figure 5. Critical Rayleigh number ( R c ) for changes in the anisotropic extensional vis-cosity ( µ ), the anisotropic shear viscosity ( µ ), and the preferred direction in the fluidat steady state ( θ (0) ) with both boundaries rigid. In each subfigure the arrows indicateincreasing µ ( µ = 0 , , , , , µ = 0, (b) µ = 10, (c) µ = 100,and (d) µ = 1000.number for the Newtonian/isotropic case ( R c ≈ R N ). Near horizontal or vertical fibre-orientation, increasing the anisotropic extensional viscosity increases the threshold atwhich instability occurs, but when the steady state preferred direction is π/ µ is increased, R c increases regularly, smoothing out the pointsof inflection that occur for small values of µ . Therefore, increasing µ stabilises thesteady state for all values of θ (0) and µ . Changes in anisotropic shear viscosity affect themagnitude of the critical Rayleigh number much more than changes to the anisotropicextensional viscosity.Figure 7 shows the dependence of R c on θ (0) for selected values of µ and µ , with thelower boundary rigid and the upper free. In Figure 7(a) we examine how µ affects R c when µ = 0; the Newtonian/isotropic case is shown by the lowest horizontal line. As µ 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, localminima at θ (0) ≈ π/
16 and 7 π/
16, and the global minimum at θ (0) = π/
4. Therefore4
Holloway, Smith and
Dyson (a) (b) θ (0) R c µ = 00 π/ π/ . . . . × θ (0) R c µ = 100 π/ π/ . . . . × (c) (d) θ (0) R c µ = 1000 π/ π/ . . . . . × θ (0) R c µ = 10000 π/ π/ . . . . . × Figure 6. Critical Rayleigh number ( R c ) for changes in the anisotropic extensional vis-cosity ( µ ), the anisotropic shear viscosity ( µ ), and the preferred direction in the fluidat steady state ( θ (0) ) with both boundaries free. In each subfigure the arrows indicateincreasing µ ( µ = 0 , , , , , µ = 0, (b) µ = 10, (c) µ = 100,and (d) µ = 1000.increasing the anisotropic extensional viscosity increases the stability threshold mostwhen, at steady state, the fibres are either horizontal or vertical, and least when they aredirected at an angle π/ µ increases, R c is increased, with the additionallocal maxima and minima becoming less pronounced, and disappearing completely once µ (cid:38) µ affect R c far more than similar changes to µ . Therefore increases in the anisotropic shearviscosity stabilise the steady state, with the most stabilisation occurring when the fibresare oriented horizontally or vertically at steady state, at least least when θ (0) = π/ π/ (cid:54) θ (0) (cid:54) π/ inear Rayleigh-B´enard stability of a TI fluid (a) (b) θ (0) R c µ = 00 π/ π/ . . . . × θ (0) R c µ = 100 π/ π/ . . . . × (c) (d) θ (0) R c µ = 1000 π/ π/ . . . . . . × θ (0) R c µ = 10000 π/ π/ . . . . . . × Figure 7. Critical Rayleigh number ( R c ) for changes in the anisotropic extensional vis-cosity ( µ ), the anisotropic shear viscosity ( µ ), and the preferred direction in the fluidat steady state ( θ (0) ) with the bottom boundary rigid and the top boundary free. Ineach subfigure the arrows indicate increasing µ ( µ = 0 , , , , , µ = 0, (b) µ = 10, (c) µ = 100, and (d) µ = 1000.anisotropic extensional viscosity, and all changes are dampened as the anisotropic shearviscosity is increased, similarly to the matching boundary case. Therefore, in all cases,the anisotropic extensional viscosity gives rise to variations in the critical wave-numberwith respect to the steady state preferred direction, which are dampened by increases inthe anisotropic shear viscosity.Comparing Figures 5 – 7 allows us to compare how the different boundary conditionsaffect the critical Rayleigh number. Similarly to the Newtonian/isotropic case the moststable pair of boundaries is rigid-rigid, with the most unstable being free-free. In allboundary pairs increasing either the anisotropic extensional or shear viscosities increasesthe critical Rayleigh number, however changes to the anisotropic shear viscosity affectthe stability threshold much more than equivalent changes to the anisotropic extensionalviscosity.We notice in Figures 5 – 7 that the critical wave and Rayleigh numbers are the samefor θ (0) and π/ − θ (0) , i.e. the material has the same stability characteristics when thesteady state preferred direction is horizontal or vertical, the dependance on this angle6 Holloway, Smith and
Dyson
Table 1. The parameter values from curve fitting for the critical wave-number given inequation (5.1) for the different combinations of boundary conditions.
Boundary Type f f f f f rigid-rigid 175 . . . × − . . × − rigid-free 91 . . . × − . . × − free-free 284 . . . × − . . × − being symmetric about π/
4. When the transversely-isotropic material is interpreted as asuspension of elongated particles, this result may be explained by noting that when theparticles are either horizontal or vertical, only translational, non-rotating, motion will beinduced. The stability characteristics of these states should therefore be similar.
Examining Figures 2 – 7 we notice that, for medium to large values of the anisotropicshear viscosity, the critical curves are continuous with no sharp extrema ( i.e. we expectthe rate of change of the critical values with θ (0) to be continuous also). We may thereforeattempt to fit analytic functions to the numerical results, for k c and R c . The critical wave-number is fit by minimising the maximum absolute error between the function and thenumerical 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 nonlinearleast squares fitting ( fit in MATLAB).For µ >
40 we fit the critical wave-number to the empirical form k c ( µ , µ , θ (0) ) ≈ f (cid:0) exp ( − f µ ) − (cid:1) f + µ cos(4 θ (0) ) − f µ f + µ + k N , (5.1)where k N is the critical wave-number for a Newtonian fluid and f to f are fittingparameters. The form of this function implies that for µ (cid:29) µ , and µ (cid:29)
1, the criticalwave-number approximates its Newtonian value. For µ = 0 the critical wave-number k c does not depend on θ (0) or µ , and takes the same value as in the Newtonian case. Figure8 shows a comparison between the sampled numerical results and fitted function, wherethe fitted parameters are given in table 1; excellent qualitative and good quantitativeagreement is found.For µ >
40 we fit the critical Rayleigh number to the empirical form R c ( µ , µ , θ (0) ) ≈ (cid:18) − g µ + 1 + g (cid:19) µ cos(4 θ (0) ) + (cid:32) − g µ / + 1 + g (cid:33) µ + ( g µ + g ) , (5.2)where g to g are fitting parameters. The form of this function implies that R c is depen-dent upon θ (0) for µ (cid:29) µ , however this dependence is dampened as µ is increased. Wealso observe from the relative sizes of g and g , in table 2, that the increase of R c with µ is much greater than that with µ . Figure 9 shows a comparison between the samplednumerical results and fitted function, with the fitted values found in table 2; a reasonable inear Rayleigh-B´enard stability of a TI fluid ( a ) ( b ) µ k c . . . . . µ k c . . . . . ( c ) ( d ) µ k c . . . . . µ k c . . . . . . Figure 8. Comparison of fitted curves with the numerical results for the critical wave-number k c . The fit for (a) ( µ = 50, θ (0) = 0) and (b) ( µ = 1000, θ (0) = 0), (c)( µ = 50, θ (0) = π/ µ = 1000, θ (0) = π/
4) where black, blue and red correspondto rigid-rigid ( RR ), rigid-free ( RF ) and free-free ( F F ) boundary pairs and circles repre-sent numerical results. We choose θ (0) = 0 , π/ θ (0) .Table 2. The parameter values from curve fitting for the critical Rayleigh number inequation (5.2) for the different combinations of boundary conditions. Boundary Type g g g g g g rigid-rigid 448 . . . . . . . . . . . . . . . . . . quantitative agreement and a good qualitative agreement is found. Note that g playsthe role of R N , and numerically is extremely close to its known Newtonian values.8 Holloway, Smith and
Dyson ( a ) ( b ) µ R c . . . × RRRFFF µ R c . . × RR RFFF ( c ) ( d ) µ R c . . × RRRFFF µ R c . . × RRRFFF
Figure 9. Comparison of fitted curves with the numerical results for the critical Rayleighnumber R c . The fit for (a) ( µ = 50, θ (0) = 0) and (b) ( µ = 1000, θ (0) ) (c) ( µ = 50, θ (0) = π/
4) (d) ( µ = 1000, θ (0) = π/ RR ), rigid-free ( RF ) and free-free ( F F ) boundary pairs and circles representthe numerical results. We choose θ (0) = 0 , π/ θ (0) . In this paper, we extended the work of Rayleigh to study the linear stability of atransversely-isotropic viscous fluid, contained between two horizontal boundaries, whichare either rigid or free, of different temperatures. We used the stress tensor first proposedby Ericksen (1960), with µ = 0 (equivalent to a passive fluid (Holloway et al., 2017)),and a kinematic equation for the fibre-director field to model a transversely-isotropicfluid. 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 rangeof directions as the governing equations have a period of π/ µ , is much more important in determiningthe stability of the flow than the anisotropic extensional viscosity µ . The influence of thispair of parameters upon the stability of the flow depends on the uniform steady state inear Rayleigh-B´enard stability of a TI fluid θ (0) , as well as the boundary conditions. Similarly to a Newtonianfluid, the most stable pair of boundaries is rigid-rigid, for which the temperature differencebetween the two boundaries required to induce instability is the largest; the least stableboundary pair is free-free.The rheological parameters µ and µ also have an impact on the critical wave-number k c , which describes the width of the convection cells. We find for steady state preferreddirections near horizontal or vertical, the width of the convection cell increases with theanisotropic extensional viscosity, when compared to a Newtonian fluid, and decreaseswhen the preferred direction makes an angle of π/ µ (cid:29) µ and µ (cid:29) Acknowledgment
CRH is supported by an Engineering and Physical Sciences Research Council (EPSRC)doctoral training award (EP/J500367/1) and RJD the support of the EPSRC grant(EP/M00015X/1). The authors thank Gemma Cupples for valuable discussions.
References
D. J. Acheson (1990).
Elementary fluid dynamics . Oxford University Press.H. B´enard (1901). ‘Les tourbillons cellulaires dans une nappe liquide.-M´ethodes optiquesd’observation et d’enregistrement’.
J. Phys.–Paris (1):254–266.S. Chandrasekhar (2013). Hydrodynamic and Hydromagnetic Stability . Courier DoverPublications.G. Cupples, et al. (2017). ‘Viscous propulsion in active transversely isotropic media’.
J.Fluid Mech. :501–524.0
Holloway, Smith and
Dyson
T. R. Dafforn, et al. (2004). ‘Protein fiber linear dichroism for structure determinationand kinetics in a low-volume, low-wavelength Couette flow cell’.
Biophys. J. (1):404–410.M. A. Dominguez-Lerma, et al. (1984). ‘Marginal stability curve and linear growthrate for rotating Couette–Taylor flow and Rayleigh–B´enard convection’. Phys. Fluids (4):856–860.P. G. Drazin (2002). Introduction to hydrodynamic stability . Cambridge University Press.R. J. Dyson, et al. (2015). ‘An investigation of the influence of extracellular matrixanisotropy and cell–matrix interactions on tissue architecture’.
J. Math. Biol. :1775–1809.R. J. Dyson & O. E. Jensen (2010). ‘A fibre-reinforced fluid model of anisotropic plantcell growth’. J. Fluid Mech. :472–503.J. L. Ericksen (1960). ‘Transversely isotropic fluids’.
Colloid. Polym. Sci. (2):117–122.R. J. Goldstein & D. J. Graham (1969). ‘Stability of a horizontal fluid layer with zeroshear boundaries’.
Phys. Fluids (6):1133–1137.J. E. F. Green & A. Friedman (2008). ‘The extensional flow of a thin sheet of incom-pressible, transversely isotropic fluid’. Eur. J. Appl. Math. (3):225–258.J. Hoepffner (2007). ‘Implementation of boundary conditions’. . [Online; accessed 03-Aug-2017].C. R. Holloway, et al. (2017). ‘Influences of transversely-isotropic rheology and transla-tional diffusion on the stability of active suspensions’. arXiv:1607.00316 .C. R. Holloway, et al. (2015). ‘Linear Taylor–Couette stability of a transversely isotropicfluid’. Proc. R. Soc. Lond. A (2178):20150141.E. L. Koschmieder (1993).
B´enard cells and Taylor vortices . Cambridge University Press.K. Kruse, et al. (2005). ‘Generic theory of active polar gels: A paradigm for cytoskeletaldynamics’.
Euro. Phys. J. E (1):5–16.J. C. Lagarias, et al. (1998). ‘Convergence properties of the Nelder–Mead simplex methodin low dimensions’. SIAM J. Control (1):112–147.M. E. M. Lee & H. Ockendon (2005). ‘A continuum model for entangled fibres’. Eur. J.Appl. Math. (2):145–160.R. Marrington, et al. (2005). ‘Validation of new microvolume Couette flow linear dichro-ism cells’. Analyst (12):1608–1616.J. R. A. McLachlan, et al. (2013). ‘Calculations of flow-induced orientation distributionsfor analysis of linear dichroism spectroscopy’.