Thin film extensional flow of a transversely isotropic viscous fluid

M. J. Hopwood School of Mathematical Sciences, University of Adelaide, Adelaide, SA, 5005, Australia School of Mathematics, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK B. Harding School of Mathematical Sciences, University of Adelaide, Adelaide, SA, 5005, Australia School of Mathematics and Statistics, Victoria University of Wellington, PO Box 600, Wellington, 6140, New Zealand J.E.F. Green School of Mathematical Sciences, University of Adelaide, Adelaide, SA, 5005, Australia R.J. Dyson School of Mathematics, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK

(July 29, 2025)

Abstract

Many biological materials such as cervical mucus and collagen gel possess a fibrous micro-structure. This micro-structure affects the emergent mechanical properties of the material, and hence the functional behaviour of the system. We consider the canonical problem of stretching a thin sheet of transversely-isotropic viscous fluid as a simplified version of the spinnbarkeit test for cervical mucus.

We propose a novel solution to the model constructed by Green & Friedman by manipulating the model to a form amenable to arbitrary Lagrangian–Eulerian (ALE) techniques. The system of equations, reduced by exploiting the slender nature of the sheet, are solved numerically and we discover that the bulk properties of the sheet are controlled by an effective viscosity dependent on the evolving angle of the fibres. In addition, we confirm a previous conjecture by demonstrating that the centre-line of the sheet need not be flat, and perform a short timescale analysis to capture the full behaviour of the centre-line.

1 Introduction

Fluids can be classified as isotropic, or anisotropic, depending upon whether the mechanical properties of the fluid are uniform in all directions, or not. Most common fluids such as water are isotropic. However, there are many examples of fluids that arise in biology and industry which contain fibres or elongated particles, for example, collagen gels [12], cervical mucus [6], and nematic liquid crystals [5]. The presence of fibres or particles within the fluid creates an underlying structure that causes the fluid to exhibit anisotropy. The anisotropy of biological fluids mean that they can display interesting behaviours and possess unusual and evolving mechanical properties, which influence how they perform their particular functions. For example, the anisotropy of cervical mucus is suspected to play a role in fertility, by controlling how easily sperm can migrate through to the egg [4].

One of the first models of an anisotropic viscous fluid was formulated by Ericksen [10]. He considered a type of anisotropy known as ‘transverse isotropy’ where the material possesses a single preferred direction which may vary both spatially and temporally; its physical properties are symmetric in all directions normal to this preferred direction. Whilst in general, for a fibrous material, one would need to consider the evolution (in time and space) of a probability distribution describing the alignment of a fibre along a particular direction, the Ericksen model simplified the problem by assuming that there is a unique fibre alignment at each point in space. Examples of materials that have been modelled using this approach include fibre–reinforced composites [28], entanglements of textile fibres used in the carding process [19] and a number of biological materials – collagen gels [12], the extra–cellular matrix [8], and primary plant cell walls [9]. A significant number of studies motivated by problems in composites manufacturing include the additional assumption that the fluid is inextensible in the fibre direction [18, 26, 28]; these are termed ‘ideal’ incompressible fibre–reinforced fluids.

In this paper, we consider the extensional flow of a thin sheet of transversely isotropic viscous fluid. One motivation for studying this problem is that it provides a simplified representation of the ‘spinnbarkeit’ or ‘spinnability’ test, which is applied to cervical mucus as a means of assessing fertility [11]. The test entails taking a sample of mucus and stretching it. Around ovulation, the mucus has a lower pH, a higher concentration of water (which has the effect of lowering the viscosity of the mucus), and the fibrous reinforcement takes a more parallel alignment that allows sperm to migrate. At this point the mucus can be stretched the furthest, i.e. the fluid exhibits maximum spinnbarkeit. Conversely, during the most infertile parts of the menstrual cycle, the mucus does not stretch and simply breaks [20, 29]. Understanding the dynamics of stretching a transversely isotropic sheet can provide some basic insights into how factors such as fibre concentration and alignment influence spinnbarkeit. A second motivation for studying this problem is that is provides a simple prototype situation for investigating how the feedback between the macroscopic fluid flow and the microstructure influences the mechanical behaviour of anisotropic fluids. In particular, as shown in [12], the model can be reduced to one–dimension, with the changing orientation of the fibres specified by a single angle.

Extensional thin film flows of incompressible Newtonian fluids, which arise in a number of industrial and biological applications, have been extensively studied. The problem of a two dimensional thin film was considered by Howell [17], who employed an asymptotic expansion of the Navier–Stokes and mass conservation equations in powers of an inverse aspect ratio, to obtain a reduced model (termed the Trouton model) which involves only the leading order longitudinal fluid velocity and sheet thickness. Additionally, the roles of inertia and surface tension were considered, as well as the complementary problem of a fluid thread (i.e. slender cylinder). In order to capture the full behaviours of the centre–line of the fluid of the sheet, a short timescale analysis is carried out [2, 16, 17]. These models are not valid for sheets undergoing deformation by bending, and led to the development of a general theory for thin viscous sheets that undergo deformation by stretching, bending, or an arbitrary combination of both, by Ribe [24, 25]. This work was able to explicitly quantify relations between bending and stretching of the sheet. A complimentary study considered sheets with an inhomogenous viscosity, [21]. One of their main results was that ‘necking’ of the sheet, where regions of the sheet thinned faster than others, could be induced by in–plane variations of viscosity throughout the fluid. However, for a transversely isotropic fluid, we are only aware of three studies of extensional flow: those of [3, 9], wherein both works modelled the primary plant cell wall as a thin axisymmetric fibre–reinforced viscous sheet supported between rigid end plates, and that of [12], who used a similar approach to that of Howell [16] to derive a transversely isotropic version of the Trouton model. The model presented in [12] is referred to throughout as the Green & Friedman model and is of particular relevance to this paper. They presented some analytical results for cases in which the equations simplify (e.g. when some of the anisotropic stress terms are negligible) but did not tackle the general case – a problem we pursue in this paper.

A number of studies have investigated transversely isotropic fluid flows in other geometries to gain insight into how their behaviour is affected by the microstructure. These include the studies of the stability of a transversely isotropic fluid in a Taylor–Couette device with the aim of understanding the behaviour of suspensions of bio-molecules, as well as treating the problem of Rayleigh–Bernard convection [14, 15]. It was found in both works that transversely isotropic effects delay the onset of instabilities, primarily through the incorporation of an anisotropic shear viscosity. In the context of modifying transversely isotropic fluid models to incorporate active swimming suspensions, as in [6, 13], transversely isotropic effects are also capable of increasing the size of a developing instability, in particular where translation diffusion is neglected. Other studies have focused on prototypical flows, to give more generic insights into fluid-fibre interactions. For example, Phan–Thien and Graham studied both the flow of a transversely isotropic fluid around a sphere [23], and the squeezing flow of a layer of fluid compressed between two fixed plates [22] (this latter problem was studied independently for the case of an ideal fibre reinforced fluid [26]).

This paper uses a combination of analysis and numerical simulations to significantly extend the previous work of [12] to include cases where all of the anisotropic terms in the fluid stress are non-negligible, and the fibre alignment within the sheet may vary with depth. One issue of particular interest was to verify their conjecture that, in contrast to the Newtonian case, the centre–line of a transversely isotropic fluid sheet need not always be straight. As preliminary simulation of the model indeed produced results with non-flat centre–lines, we sought to confirm this prediction by additionally investigating the short timescale behaviour of the sheet, similar to the work by [2, 17] on the Newtonian problem.
The paper is organised as follows. In Section 2 we briefly recap the thin film approach and governing equations as given in [12], as well as introducing a short timescale. In order to present the first solutions to the full Green & Friedman model, the equations are manipulated in order become amenable to numerical strategies in section 3. We then present the arbitrary Lagrangian–Eulerian techniques we use to solve the model in Section 4, and validate the numerical techniques by comparison with analytical results for short time, and further present new results primarily for a passive transversely isotropic fluid. In Section 5 we present results for the behaviour of the fluid on a short timescale. We conclude with a discussion and suggestions for future work in Section 6.

2 Mathematical model

We consider the extensional flow of an inertialess thin sheet of incompressible, transversely isotropic, viscous fluid. As shown in Figure 1, we use the 2D Cartesian coordinates $\left(x^{*},y^{*}\right)$ to describe the horizontal and vertical axis respectively, with $t^{*}$ denoting time (throughout this paper asterisks denote dimensional quantities). The upper and lower boundaries of the fluid sheet are denoted by $y^{*}=H^{\pm^{*}}=H^{*}\pm\frac{h^{*}}{2}$ , where $H^{*}(x^{*},t^{*})$ is the position of the centre–line and $h^{*}(x^{*},t^{*})$ is the thickness of the fluid sheet. The left– and right– hand side boundaries are given by $x^{*}=0,L^{*}(t^{*})$ . The right-hand of the sheet, at $x^{*}=L^{*}$ , is pulled in the $x^{*}$ direction; we prescribe either the speed of pulling, $\dot{L}^{*}$ , or the tension $T^{*}$ applied to the sheet.

We let $\boldsymbol{u}^{*}=\left(u^{*},v^{*}\right)$ be the fluid velocities in the $x^{*},y^{*}$ directions respectively, and denoted the stress tensor by $\boldsymbol{\sigma}^{*}$ . The equations of conservation of fluid mass and momentum are thus

	$\displaystyle\nabla^{}\cdot\boldsymbol{u}^{}=0,$		(2.1)
	$\displaystyle\nabla^{}\cdot\boldsymbol{\sigma}^{}=0.$		(2.2)

The constitutive law for $\boldsymbol{\sigma}^{*}$ is

\displaystyle\sigma^{*}_{ij}=-p^{*}\delta_{ij}+2\mu^{*}e^{*}_{ij}+\mu^{*}_{1}a_{i}a_{j}+\mu^{*}_{2}a_{i}a_{j}a_{k}a_{l}e^{*}_{kl}+2\mu_{3}^{*}(a_{i}a_{l}e^{*}_{jl}+a_{j}a_{m}e^{*}_{mi}),

(2.3)

where $p^{*}$ is the pressure, $\boldsymbol{a}$ is a unit vector describing the orientation of the fibres within the fluid and $\boldsymbol{e}^{*}$ is the rate of strain tensor. This relationship was derived by Ericksen [10], as the most general stress tensor that is linear in the rate of strain tensor $\boldsymbol{e}^{*}$ , invariant under the transformation $-\boldsymbol{a}\rightarrow\boldsymbol{a}$ , and satisfies $\boldsymbol{\sigma}^{*}=\boldsymbol{\sigma}^{*T}$ . The constants $\mu^{*},\mu^{*}_{2},\mu^{*}_{3}$ , are all viscosity-like parameters and $\mu^{*}_{1}$ is the active tension in the fibre direction [3, 12, 27].

We note first that by setting $\mu^{*}_{1}=\mu^{*}_{2}=\mu^{*}_{3}=0$ , one immediately recovers the stress tensor for an incompressible, isotropic Newtonian fluid, with $\mu^{*}$ as the familiar dynamic (shear) viscosity. The physical interpretations of $\mu_{2}^{*}$ and $\mu_{3}^{*}$ can be identified by considering three deformations of a $2D$ sheet of fibres in a Cartesian plane, as illustrated in [9, 12, 14]. In the same way as previous work, we interpret $\mu_{2}^{*}$ as the anisotropic extensional viscosity, and $\mu_{3}^{*}$ the anisotropic shear viscosity. Additionally, we observe that there is no velocity component to the $\mu_{1}^{*}$ term, indicating that there exists stress in the fluid even when at rest and that this stress must be a tensile stress as no stress is induced in the fibres when the fluid is compressed.

Figure 1: Schematic of extensional flow of a sheet. The fluid is fixed to two plates at

x^{*}=0,L

, by a no-slip boundary condition. The plate at

x^{*}=L^{*}

is moved at a prescribed speed

\dot{L}^{*}

. The centre–line of the fluid is given by

H^{*}(x^{*},t^{*})

and the thickness of the film by

h^{*}(x^{*},t^{*}),

so that the free boundaries of the film are located at

y^{*}=H^{*}\pm\frac{h^{*}}{2}

. Note that this figure is not to scale, as we model the thin film behaviour of the sheet.

In addition to the constitutive law for $\boldsymbol{\sigma}^{*}$ , we require an equation governing the evolution of the fibre direction. We use the form given by Green & Friedman (for derivation, see [12]), which results from allowing the fibres to advect with the fluid flow

\displaystyle\frac{\partial\boldsymbol{a}}{\partial t^{*}}+(\boldsymbol{u}^{*}\cdot\nabla^{*})\boldsymbol{a}+\zeta^{*}\boldsymbol{a}=(\boldsymbol{a}\cdot\nabla^{*})\boldsymbol{u}^{*},

(2.4)

where

\displaystyle\zeta^{*}(x^{*},y^{*},t^{*})=\boldsymbol{a}\cdot(\boldsymbol{a}\cdot\nabla^{*}\boldsymbol{u}^{*}),

(2.5)

is the fractional rate of extension of the fibres. The first two terms are the convective derivative of the fibre orientation, the third is the fractional rate of extension of the fibre in the direction of the fibres, all of which balances with effect of the fibres upon the flow velocity field on the right hand side. Equation (2.4) corresponds to the specific case of Ericksen’s equation in the long fibre limit [6, 10, 12]. Since $\boldsymbol{a}$ is a unit vector and our model is two dimensional we write $\boldsymbol{a}=\left(\cos\theta,\sin\theta\right)$ , where $\theta\left(x^{*},y^{*},t^{*}\right)$ is the angle the fibre direction makes with the $x^{*}$ -axis.

In order to close the model, we must impose suitable boundary and initial conditions. At the ends of the sheet, we set

	$\displaystyle u^{}\left(0,y^{},t^{*}\right)=0,$	$\displaystyle u^{}\left(L^{},y^{},t^{}\right)=\dot{L}^{*},$
	$\displaystyle H^{}\left(0,t^{}\right)=0,$	$\displaystyle H^{}\left(L^{},t^{*}\right)=0.$

On the upper and lower free surfaces, we apply a no-stress boundary condition

\displaystyle\boldsymbol{\sigma^{*}}\cdot\hat{\boldsymbol{n}}=0;\text{ on }y^{*}=H^{*}\pm\frac{1}{2}h^{*},

together with the usual kinematic condition

\displaystyle v^{*}=\frac{\partial H^{*}}{\partial t^{*}}\pm\frac{1}{2}\frac{\partial h^{*}}{\partial t^{*}}+u^{*}\left(\frac{\partial H^{*}}{\partial x^{*}}\pm\frac{1}{2}\frac{\partial h^{*}}{\partial x^{*}}\right);\text{ on }y=H^{*}\pm\frac{1}{2}h^{*}.

(2.6)

Initial conditions must also be prescribed, however, since the number of initial conditions required depends upon the timescale considered, this will be discussed later.

2.1 Thin film approximation

We now introduce the assumption that the sheet is thin, which allows considerable simplification of the governing equations. Full details of the derivation can be found in [12], but for the sake of completeness, we recapitulate the main points here. We let $L_{0}$ and $h_{0}$ be the initial length and typical initial thickness of the fluid sheet, respectively, and let $U$ be a typical value for the velocity of the fluid at the pulled boundary. We then introduce the parameter $\varepsilon=h_{0}/L_{0}\ll 1$ , which is the initial inverse aspect ratio of the sheet. We are interested in the behaviour of the sheet as it undergoes significant changes in length, and so consider the timescale $t^{*}\sim L_{0}/U$ . Following [12, 16] we nondimensionalise as follows:

	$\displaystyle\left(x^{},y^{}\right)=\left(xL_{0},\varepsilon yL_{0}\right),\quad\left(u^{},v^{}\right)$	$\displaystyle=\left(uU,\varepsilon vU\right),\quad p^{}=\dfrac{\mu^{}U}{L_{0}}p,\quad t^{*}=\frac{L_{0}}{U}t,$
	$\displaystyle\left(H,L,h\right)=$	$\displaystyle\left(\varepsilon L_{0}H^{},L_{0}L^{},\varepsilon L_{0}h^{*}\right).$

These scalings introduce the dimensionless material parameters

\displaystyle\mu_{1}=\frac{\mu_{1}^{*}L}{\mu^{*}U},\quad\mu_{2}=\frac{\mu_{2}^{*}}{\mu^{*}},\quad\mu_{3}=\frac{\mu_{3}^{*}}{\mu^{*}}.

At this point, we exploit the thin geometry of the sheet by expanding all of the dependent variables as power series in terms of the inverse aspect ratio $\varepsilon$ , that is

\displaystyle u=u^{(0)}+\varepsilon u^{(1)}+...

with similar expressions for the other dependent variables (here, unlike in Howell [16], we encounter terms involving odd powers of $\varepsilon$ ). After some lengthy algebra, full details of which are given in [12], we obtain a system of one dimensional equations for the quantities $h^{(0)},H^{(0)},u^{(0)},u^{(1)},v^{(0)},\theta^{(0)}$ . As a consequence of the analysis, it is found that the leading order longitudinal velocity satisfies $u^{(0)}=u^{(0)}(x,t)$ only (i.e. the flow is extensional). Conservation of mass yields:

\displaystyle\frac{\partial h^{(0)}}{\partial t}+

\displaystyle\frac{\partial}{\partial x}\left(h^{(0)}u^{(0)}\right)=0.

(2.7)

Taking a depth-averaged force balance over the sheet leads to the following equation for $u^{(0)}$

\displaystyle\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}4\left(1+\mu_{3}\right)\frac{\partial u^{(0)}}{\partial x}+

\displaystyle\mu_{1}\cos 2\theta^{(0)}+\mu_{2}\left(\cos^{2}2\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\frac{1}{4}\sin 4\theta^{(0)}\frac{\partial u^{(1)}}{\partial y}\right)\mathrm{d}y=0.

(2.8)

Consideration of the momentum equations and associated no-stress boundary conditions on the upper and lower boundaries of the fluid at higher order gives an equation governing the centre–line of the fluid, $H^{(0)}$

\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{y}4\left(1+\mu_{3}\right)\frac{\partial u^{(0)}}{\partial x}+\mu_{1}\cos 2\theta^{(0)}\\ +\mu_{2}\left(\cos^{2}2\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\frac{1}{4}\sin 4\theta^{(0)}\frac{\partial u^{(1)}}{\partial y}\right)\mathrm{d}s\mathrm{d}y=0.

(2.9)

We note first that $s$ is a dummy variable, we are integrating over the second argument of the functions in the integrand twice, so that the quantities in the integrand of equation (2.9) are $\theta^{(0)}\left(x,s,t\right),u^{(1)}\left(x,s,t\right)$ (to emphasise this, the dummy variable is not used in [12]). Additionally, we note that equation (2.9) is a particularly unusual form of an integro-differential equation, wherein one of the variables of integration appears in the limit of one of the integrals. The equation (2.4) yields

\displaystyle\frac{\partial\theta^{(0)}}{\partial t}+u^{(0)}\frac{\partial\theta^{(0)}}{\partial x}+v^{(0)}\frac{\partial\theta^{(0)}}{\partial y}=

\displaystyle-2\sin\theta^{(0)}\cos\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}-\sin^{2}\theta^{(0)}\frac{\partial u^{(1)}}{\partial y},

(2.10)

as an evolution equation governing the behaviour of the fibre director field. The leading order transverse velocity is found by integrating the incompressiblity condition and applying the kinematic condition (2.6)

\displaystyle v^{(0)}=

\displaystyle\frac{\partial H^{(0)}}{\partial t}+\frac{\partial}{\partial x}\left(H^{(0)}u^{(0)}\right)-y\frac{\partial u^{(0)}}{\partial x},

(2.11)

consideration of the $x$ -momentum equation at $\mathcal{O}\left(\varepsilon\right)$ with its associated boundary condition supplies

\frac{\partial u^{(1)}}{\partial y}=-\mu_{1}\frac{2\sin 2\theta^{(0)}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta^{(0)}}-\mu_{2}\frac{\sin 4\theta^{(0)}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta^{(0)}}\frac{\partial u^{(0)}}{\partial x},

(2.12)

as the next-order correction term for $u$ . In this form, one can observe that in the case of $\mu_{1}=\mu_{2}=\mu_{3}=0$ , we return to the Trouton model for a Newtonian fluid. In the case of $\mu_{1}=\mu_{2}=0$ , the fluid behaves very similarly to a Newtonian fluid, with the only modification being a modified Trouton ratio, which effectively only changes the tension required to be applied to the sheet to achieve the same behaviour as a Newtonian fluid. In these cases, the model can be solved by means of a Lagrangian transformation, as detailed in [12, 16].

The boundary conditions for the leading order longitudinal velocity are then

\displaystyle u^{(0)}(0,t)=0,\quad u^{(0)}(L,t)=\dot{L}.

(2.13)

We assume that the sheet is being extended in the $x$ direction and that the end points of the centre–line remain fixed to $y=0$ , thus

\displaystyle H^{(0)}(0,t)=H^{(0)}(L,t)=0.

(2.14)

The kinematic boundary condition yields

\displaystyle v^{(0)}=\frac{\partial H^{(0)}}{\partial t}\pm\frac{1}{2}\frac{\partial h^{(0)}}{\partial t}+u^{(0)}\left(\frac{\partial H^{(0)}}{\partial x}\pm\frac{1}{2}\frac{\partial h^{(0)}}{\partial x}\right);\text{ on }y=H^{(0)}\pm\frac{1}{2}h^{(0)}.

(2.15)

We will need to prescribe initial conditions for the thickness, $h^{(0)},$ and the fibre direction, $\theta^{(0)}$ in the sheet. As in the Newtonian problem (see [16]), we do not need to prescribe an initial condition for $H^{(0)}$ as we are unable to satisfy an arbitrary initial condition for $H^{(0)}$ , as noted in [12]. In order to study the behaviour of sheets that do not initially obey equation (2.9) we must consider a shorter timescale than $\dfrac{L_{0}}{U}$ . We turn to such a timescale in the next subsection.

2.2 Short timescale model

Similarly to the Newtonian problem, we have a singular perturbation problem for $H^{(0)}$ in $t$ , the outer solution of which is determined by the solution of (2.9). As suspected by Green & Friedman [12], we discover that the centre–line is not necessarily straight (unlike the Newtonian case). In order to study the evolution of a sheet with an arbitrary initial condition for $H^{(0)}$ , we follow a similar approach to Howell and Buckmaster et al. [2, 16] by considering a timescale of $\varepsilon^{2}\dfrac{L_{0}}{U}$

We note from (2.15), that in addition to rescaling time, we must also rescale the velocity $v$ ; hence, we introduce

\tau=\frac{t}{\varepsilon^{2}},\quad v=\frac{V}{\varepsilon^{2}}.

(2.16)

As in the previous section, we exploit the slender geometry of the sheet by expanding variables as a power series of $\varepsilon$ . At leading order the continuity equation gives

\displaystyle\frac{\partial V^{(0)}}{\partial y}=0,

(2.17)

which, combined with the kinematic condition at the same order

\displaystyle V^{(0)}=\frac{\partial H^{(0)}}{\partial\tau}\pm\frac{1}{2}\frac{\partial h^{(0)}}{\partial\tau};\text{ on }y=H^{\pm},

(2.18)

tells us

\displaystyle V^{(0)}=\frac{\partial H^{(0)}}{\partial\tau},\qquad\frac{\partial h^{(0)}}{\partial\tau}=0,

(2.19)

since if $V^{(0)}$ is independent of $y$ , it must take the same value on the upper and lower free surfaces. Hence, we have an expression for $V^{(0)}$ in terms of $H^{(0)}$ , and note that there is no thinning of the sheet on this timescale. Additionally we obtain the equation for $\theta^{(0)}$

\frac{\partial\theta^{(0)}}{\partial\tau}+\frac{\partial H^{(0)}}{\partial\tau}\frac{\partial\theta^{(0)}}{\partial y}=0.

(2.20)

Consideration of the $x$ -momentum equation and the associated no-stress boundary condition at $\mathcal{O}\left(\varepsilon\right)$ yields a result for $u^{(0)}$ :

\displaystyle u^{(0)}=\bar{u}\left(x,\tau\right)+\left(H^{(0)}-y\right)\frac{\partial^{2}H^{(0)}}{\partial\tau\partial x},

(2.21)

where $\bar{u}(x,t)$ is an as yet unknown function arising from integration. This is precisely the same result as for a Newtonian fluid [16]. Consideration of the $y$ -momentum equation at the same order yields an identity which is already satisfied. Next, the $y$ -momentum equation at $\mathcal{O}\left(\varepsilon^{2}\right)$ yields an expression for pressure

-\frac{\partial p^{(0)}}{\partial y}+\frac{\partial^{2}V^{(0)}}{\partial x^{2}}+\frac{\partial^{2}V^{(2)}}{\partial y^{2}}+\mu_{1}\frac{\partial}{\partial y}\left(\sin^{2}\theta^{(0)}\right)+\mu_{2}\frac{\partial}{\partial y}\left(\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}\right.\\ +\left.\cos\theta^{(0)}\sin^{3}\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)+\sin^{4}\theta\frac{\partial V^{(2)}}{\partial y}\right)\\ +2\mu_{3}\frac{\partial}{\partial y}\left(2\sin^{2}\theta^{(0)}\frac{\partial V^{(2)}}{\partial y}+\cos\theta^{(0)}\sin\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right)=0,

(2.22)

with the associated boundary condition

-p^{(0)}+2\frac{\partial V^{(2)}}{\partial y}+\mu_{1}\sin^{2}\theta^{(0)}+\mu_{2}\left(\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}\right.\\ \left.+\cos\theta^{(0)}\sin^{3}\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)+\sin^{4}\theta\frac{\partial V^{(2)}}{\partial y}\right)\\ +2\mu_{3}\left(2\sin^{2}\theta^{(0)}\frac{\partial V^{(2)}}{\partial y}+\sin\theta^{(0)}\cos\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right)=0;\text{ on }y=H^{(0)^{\pm}}.

(2.23)

A number of terms involving $H^{(1)}$ arise in the calculation of (2.23), these terms are multiplied by a term that is identically zero and are thus omitted. We can directly integrate (2.22) and apply (2.23) to obtain an equation for pressure

p^{(0)}=-2\frac{\partial u^{(0)}}{\partial x}+\mu_{1}\sin^{2}\theta^{(0)}+\mu_{2}\left(\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}\right.\\ \left.+\cos\theta^{(0)}\sin^{3}\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right.\\ \left.-\sin^{4}\theta\frac{\partial u^{(0)}}{\partial x}\right)+2\mu_{3}\left(-2\sin^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\sin\theta^{(0)}\cos\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right),

(2.24)

where $\mathcal{O}\left(\varepsilon^{2}\right)$ continuity has been used to eliminate the $V^{(2)}$ terms. This result for pressure is essentially a Newtonian pressure enhanced by the presence of the fibres; should we choose $\mu_{1}=\mu_{2}=\mu_{3}=0$ , we recover the Newtonian pressure. We have also introduced $u^{(1)}$ and $V^{(1)}$ terms, which must now be eliminated. The $\mathcal{O}\left(\varepsilon^{2}\right)$ $x$ -momentum equation is

\frac{\partial^{2}u^{(1)}}{\partial y^{2}}+\mu_{1}\frac{\partial}{\partial y}\left(\cos\theta^{(0)}\sin\theta^{(0)}\right)+\mu_{2}\frac{\partial}{\partial y}\left(\cos^{3}\theta^{(0)}\sin^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}\right.\\ +\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)+\left.\cos\theta^{(0)}\sin^{3}\theta^{(0)}\frac{\partial V^{(2)}}{\partial y}\right)\\ +\mu_{3}\frac{\partial}{\partial y}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)=0,

(2.25)

with the associated boundary condition

\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}+\mu_{1}\cos\theta^{(0)}\sin\theta^{(0)}+\mu_{2}\left(\cos^{3}\theta^{(0)}\sin^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}\right.\\ +\left.\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right.+\left.\cos\theta^{(0)}\sin^{3}\theta^{(0)}\frac{\partial V^{(2)}}{\partial y}\right)\\ +\mu_{3}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)=0,\text{ on }y=H^{(0)^{\pm}}.

(2.26)

Combining these yields the following compatibility condition

	$\displaystyle\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\left(1+\mu_{2}\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}+\mu_{3}\right)+\mu_{1}\cos\theta^{(0)}\sin\theta^{(0)}$
	$\displaystyle+\mu_{2}\frac{\partial u^{(0)}}{\partial x}\left(\cos^{3}\theta^{(0)}\sin\theta^{(0)}-\cos\theta^{(0)}\sin^{3}\theta^{(0)}\right)=0.$		(2.27)

We note that in the analysis, the terms $V^{(1)}$ , $u^{(1)}$ only appear together. We can thus view (2.27) as an expression which allows us to eliminate both $V^{(1)}$ and $u^{(1)}$ . Hence, we can now write leading order pressure in terms of $\theta^{(0)}$ and $u^{(0)}$ only.

In order to close the model, we must go to yet higher orders in order to obtain equations for $\bar{u}$ and $H^{(0)}$ . Our approach is similar to Green & Friedman, [12]: we integrate the relevant equations over the depth of the sheet and apply the no-stress boundary conditions at $y=H^{(0)^{\pm}}$ . The expressions involved are rather cumbersome, but substituting for previously-determined quantities produces considerable simplification; full details can be found in Appendix D. We finally obtain the following equation for $\bar{u}$

\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\dfrac{\mu_{1}\cos 2\theta^{(0)}+\left(4+4\mu_{3}+\mu_{2}\right)\left(\frac{\partial\bar{u}}{\partial x}+\left(H^{(0)}-y\right)\frac{\partial H^{{(0)}^{3}}}{\partial x^{2}\partial\tau}+\frac{\partial H^{(0)}}{\partial x}\frac{\partial^{2}H^{(0)}}{\partial x\partial\tau}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta^{(0)}}\mathrm{d}y=0.

(2.28)

The equation for $H^{(0)}$ is

\frac{\partial^{2}}{\partial x^{2}}\left(~\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\int\limits_{H^{(0)^{-}}}^{y}\dfrac{\mu_{1}\cos 2\theta^{(0)}+\left(4+4\mu_{3}+\mu_{2}\right)\left(\frac{\partial\bar{u}}{\partial x}+\left(H^{(0)}-y\right)\frac{\partial H^{{(0)}^{3}}}{\partial x^{2}\partial\tau}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta^{(0)}}\mathrm{d}s\mathrm{d}y\right)\\ =\frac{\partial}{\partial x^{2}}\left(H^{(0)^{+}}\right)\left(\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\dfrac{\mu_{1}\cos 2\theta^{(0)}+\left(4+4\mu_{3}+\mu_{2}\right)\left(\frac{\partial\bar{u}}{\partial x}+\left(H^{(0)}-y\right)\frac{\partial H^{{(0)}^{3}}}{\partial x^{2}\partial\tau}+\frac{\partial H^{(0)}}{\partial x}\frac{\partial^{2}H^{(0)}}{\partial x\partial\tau}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta^{(0)}}\mathrm{d}y\right).

(2.29)

We have now derived a system of equations for $\theta^{(0)},\bar{u}$ , and $H^{(0)}$ , namely (2.20),(2.28),(2.29) respectively. We note that (2.29) now includes a $\dfrac{\partial^{5}H^{(0)}}{\partial x^{4}\partial\tau}$ term. We must prescribe two more boundary conditions for $H^{(0)}$ , in addition to the those discussed in the previous section. We set

\displaystyle\frac{\partial H^{(0)}}{\partial x}(0,\tau)=\frac{\partial H^{(0)}}{\partial x}(L,\tau)=0.

(2.30)

Much like the Green & Friedman model, the equations (2.20),(2.28),(2.29) must be solved numerically. The approach to the short time-scale problem is more straightforward than we use for the Green & Friedman model (the details of which we will present in Section 3).

3 Arbitrary Lagrangian–Eulerian Methods

Although the Green & Friedman system, (2.7)–(2.12) and the short timescale system, (2.20), (2.28)–(2.29) are both significant simplifications compared to the full two-dimensional problem, they are too complex to allow significant analytical progress and must be solved numerically. In this section, we reformulate the Green & Friedman and short-timescale equations into an arbitrary Lagrangian–Eulerian (ALE) formulation. ALE methods involve the construction of a reference domain together with mappings from this domain to both the Lagrangian and Eulerian descriptions of the flow. Unlike numerical techniques based on either a purely Lagrangian description, where the nodes of the computational mesh follow an associated material particle throughout the motion, or upon a purely Eulerian description, where the computational mesh is fixed and the motion of the continuum is with respect to the grid, ALE methods allow freedom in moving the mesh in a way that is not necessarily fixed to a material particle. This can provide accurate solutions when modelling greater distortions of a flow problem than can ordinarily be handled by numerical techniques upon a Lagrangian description, with more resolution than is often attainable by a purely Eulerian description [7].

Our approach largely follows [7], but we outline the details here for completeness. We introduce Lagrangian, Eulerian, and reference domain variables which we denote by $\boldsymbol{X}=\left(X,Y\right)$ , $\boldsymbol{x}=\left(x,y\right)$ and $\boldsymbol{x^{\prime}}=\left(x^{\prime},y^{\prime}\right),$ respectively. Converting between Eulerian and Lagrangian descriptions of flow fields is well established, and has been employed in Newtonian extensional flow problems [12, 16, 30]. We define the map $\boldsymbol{\varphi}$ from the Lagrangian to the Eulerian descriptions such that

\left(\boldsymbol{x},t\right)=\left(\boldsymbol{\varphi}\left(\boldsymbol{X},t\right),t\right),

where the material velocity ${\boldsymbol{\upsilon}}$ is given by

\frac{\partial\boldsymbol{\varphi}}{\partial t}=\boldsymbol{\upsilon}\left(\boldsymbol{X},t\right),

(3.1)

so that the familiar material time derivative of an arbitrary scalar field (e.g pressure) is

\displaystyle\frac{Df\left(\boldsymbol{X},t\right)}{Dt}=\frac{\partial f\left(\boldsymbol{x},t\right)}{\partial t}+\frac{\partial f\left(\boldsymbol{x},t\right)}{\partial\boldsymbol{x}}\frac{\partial\boldsymbol{\varphi}}{\partial t}=\frac{\partial f}{\partial t}+\left(\boldsymbol{\upsilon}\cdot\nabla_{\boldsymbol{x}}\right)f.

(3.2)

Now, it remains only to define the map from the reference domain to the Eulerian domain, which we denote by ${\boldsymbol{\Phi}}$ . This satisfies

\left(\boldsymbol{x},t\right)=\left({\boldsymbol{\Phi}}\left(\boldsymbol{x^{\prime}},t\right),t\right).

The mesh velocity, $\boldsymbol{u}_{\text{mesh}}$ , is given by

\boldsymbol{u}_{\text{mesh}}\left(\boldsymbol{x^{\prime}},t\right)=\frac{\partial{\boldsymbol{\Phi}}}{\partial t}.

We obtain the relations between physical quantities, such as pressure, in the reference and Eulerian descriptions in the same way we do between Lagrangian and Eulerian descriptions. In particular, the time derivative of an arbitrary scalar field $f$ is

\displaystyle\frac{\partial f\left(\boldsymbol{x^{\prime}},t\right)}{\partial t}

\displaystyle=\frac{\partial f\left(\boldsymbol{x},t\right)}{\partial t}+\left(\boldsymbol{u}_{\text{mesh}}\cdot\nabla_{\boldsymbol{x}}\right)f,

(3.3)

so that the transformation from the reference to the Eulerian description behaves very much like a familiar transformation between Lagrangian and Eulerian frames. For convenience, we rewrite (3.3) as

\displaystyle\frac{\partial f\left(\boldsymbol{x^{\prime}},t\right)}{\partial t}-\frac{\partial\boldsymbol{x^{\prime}}}{\partial\boldsymbol{x}}\left(\boldsymbol{u}_{\text{mesh}}\cdot\nabla_{\boldsymbol{x^{\prime}}}\right)f=\frac{\partial f\left(\boldsymbol{x},t\right)}{\partial t}.

(3.4)

As a final note, if $\boldsymbol{\Phi}=\boldsymbol{\varphi}$ , then the reference description is the same as the Lagrangian description and thus $\boldsymbol{u}_{\text{mesh}}=\boldsymbol{\upsilon}$ . If $\boldsymbol{\Phi}=\boldsymbol{x^{\prime}}$ , then the reference description is the same as the Eulerian description, and $\boldsymbol{u}_{\text{mesh}}=\boldsymbol{0}$ . For further explanation of ALE techniques, we direct the reader to [7].

3.1 Employing ALE upon the Green & Friedman model

Henceforth we drop the superscript notation on leading order terms. We first employ ALE techniques upon the equations governing the $t\sim\dfrac{L_{0}}{U}$ timescale, equations (2.7) – (2.15). We define the reference domain, $D_{\text{ref}}$ , to be

\left(x^{\prime},y^{\prime}\right)\in D_{\text{ref}}=[0,1]\times\left[-\frac{1}{2},\frac{1}{2}\right].

(3.5)

We define the mapping from the reference variables to the Eulerian variables $\boldsymbol{\Phi}:D_{\text{ref}}\rightarrow\left[0,L\right]\times\left[H^{-},H^{+}\right]$ such that

\left(x,y\right)=\boldsymbol{\Phi}\left(x^{\prime},y^{\prime}\right)=\left(Lx^{\prime},H\left(\frac{Lx^{\prime}}{L(0)},t\right)+h\left(\frac{Lx^{\prime}}{L(0)},t\right)y^{\prime}\right).

(3.6)

This choice of map maintains equidistant spacing in the mesh in both horizontal and vertical directions. Therefore, our discretisation of $D_{\text{ref}}$ can be a simple equidistant grid. As already discussed, the mesh velocity is readily obtained by differentiating the mapping $\boldsymbol{\Phi}$ with respect to time. Written in terms of Eulerian variables, we have

\frac{\partial\boldsymbol{\Phi}}{\partial t}=\mathbf{u}_{\text{mesh}}=\left(\dot{L}\frac{x}{L},\dot{L}\frac{x}{L}\left(\frac{\partial H}{\partial x}(x,t)+\frac{y-H(x,t)}{h(x,t)}\frac{\partial h}{\partial x}(x,t)\right)+\frac{\partial H}{\partial t}(x,t)\right.\\ \left.+\frac{y-H(x,t)}{h(x,t)}\frac{\partial h}{\partial t}(x,t)\right).

(3.7)

We now account for the introduction of the moving mesh by using (3.4) and (3.6), we first modify the equation for $\theta$ , (2.10). Using (2.12) to eliminate $u^{(1)}$ , this equation becomes

\frac{\partial\theta}{\partial t}+\left(\frac{u-u_{\text{mesh}}}{L}\right)\frac{\partial\theta}{\partial x^{\prime}}+\left(\frac{v-v_{\text{mesh}}}{h}\right)\frac{\partial\theta}{\partial y^{\prime}}=\frac{\sin 2\theta\left(2\mu_{1}\sin^{2}\theta-\frac{1}{L}\left(1+2\mu_{2}\sin^{2}\theta\right)\frac{\partial u}{\partial x^{\prime}}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}.

(3.8)

This form of (3.8) allows us yet further simplification. As explicitly demonstrated in Appendix B, this equation corresponds to advection purely in a horizontal direction on the reference domain, which significantly eases implementation. We note that this also decouples $\theta$ from $H$ , and since $v$ does not arise in any other equation, we need not prescribe an initial $H$ and can simply compute the transverse velocity on demand. To summarise the model in the ALE form, we have

$\displaystyle\frac{\partial h}{\partial t}$	$\displaystyle+\left(\frac{u-u_{\text{mesh}}}{L}\right)\frac{\partial h}{\partial x^{\prime}}+\frac{1}{L}\frac{\partial u}{\partial x^{\prime}}=0,$	(3.9)
$\displaystyle\frac{\partial\theta}{\partial t}$	$\displaystyle+\left(\frac{u-u_{\text{mesh}}}{L}\right)\frac{\partial\theta}{\partial x^{\prime}}=\frac{\sin 2\theta\left(2\mu_{1}\sin^{2}\theta-\frac{1}{L}\left(4+4\mu_{3}+2\mu_{2}\sin^{2}\theta\right)\frac{\partial u}{\partial x^{\prime}}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta},$	(3.10)
$\displaystyle\frac{\partial}{\partial x^{\prime}}$	$\displaystyle\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\frac{\mu_{1}\cos 2\theta+\frac{1}{L}\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial u}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}y^{\prime}=0,$	(3.11)
$\displaystyle\frac{\partial^{2}}{\partial x^{\prime 2}}$	$\displaystyle\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{y^{\prime}}\frac{\mu_{1}\cos 2\theta+\frac{1}{L}\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial u}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}s^{\prime}\mathrm{d}y^{\prime}$
	$\displaystyle\qquad=\left(\frac{\partial^{2}H}{\partial x^{\prime 2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{\prime 2}}\right)\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\frac{\mu_{1}\cos 2\theta+\frac{1}{L}\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial u}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}y^{\prime},$	(3.12)

for $\left(x^{\prime},y^{\prime}\right)\in D_{\text{ref}}$ . The associated boundary and initial conditions now become

$\displaystyle u(0,t)=0,$	$\displaystyle u(1,t)=1,$	(3.13)
$\displaystyle H(0,t)=0,$	$\displaystyle H(1,t)=0,$	(3.14)
$\displaystyle h(x^{\prime},0)=h_{i}\left(x^{\prime}\right),$	$\displaystyle\theta\left(x^{\prime},y^{\prime},0\right)=\theta_{i}\left(x^{\prime},y^{\prime}\right).$	(3.15)

We exclude the equation for the transverse velocity $v$ , since this quantity is readily calculated from (2.11) at any $(x^{\prime},y^{\prime},t)$ once the model (3.9)–(3.12) is solved.

In order to solve this model we use the following algorithm. Given initial conditions, equation (3.11) is solved to obtain $u\left(x^{\prime},0\right)$ , using the trapezoidal rule and a finite difference approximation. We note that this particular equation can be approached in either ALE or Eulerian variables, since applying the trapezoidal rule to solve the integral (3.11) effectively decouples $H$ from the system (see Appendix C). Once solved, equations (3.9) and (3.10) can be used to update $h,\theta$ to the next time–step, at which point the process is repeated until the required end time is reached. The quantities $v,H$ can then be computed at any time using (2.11) and (3.12), respectively, as well as any other quantities of interest. In our implementation, we solved equations (3.9)–(3.12) in the ALE form given above. Since the ALE framework acts like a simple substitution on the integral equations, it is equally possible to discretise and solve equations (3.11) and (3.12) in their Eulerian forms, and indeed the discretisation is similar. We include the discretisations of equations (3.11),(3.12) in Appendix C.

3.2 Validation of numerical method

In order to validate the approach, we consider the early–time behaviour of the sheet for which an analytical expression is available [12]. Introducing the short timescale $t^{\prime}=\delta^{-1}t$ , where $\varepsilon\ll\delta\ll 1$ , assuming constant initial film thickness, $h_{i}$ , and alignment angle, $\theta_{i}$ , and performing a Taylor expansion on the variables so that

	$\displaystyle h=h_{i}+\hat{h}t^{\prime}\delta+\mathcal{O}\left(\delta^{2}\right),$
	$\displaystyle\theta=\theta_{i}+\hat{\theta}t^{\prime}\delta+\mathcal{O}\left(\delta^{2}\right),$

where $\hat{h},\hat{\theta}$ are the changes in thickness and fibre direction respectively, it is possible to derive an analytical expression for the evolution of the fibre angle in the sheet for constant initial conditions for thickness and fibre direction, see [12]. The result is:

\hat{\theta}=-2\sin 2\theta_{i}+\frac{2\sin^{2}\theta_{i}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta_{i}}\left[\mu_{1}\sin 2\theta_{i}+\frac{\mu_{2}}{2}\sin 4\theta_{i}\right],

(3.16)

where $t^{\prime}\hat{\theta}$ is the total change in the fibre angle over the short time $\tau$ , and $\theta_{i}$ is the constant initial condition. In [12], the authors note that a consequence of this analysis is that so long as $\mu_{1},\mu_{2},T$ are all sufficiently small, the fibres will align with the direction of pulling as long as $\theta_{i}\neq\frac{\pi}{2},\frac{3\pi}{2}$ . In Figure 2 we plot the numerical solution of (3.9) – (3.12) against equation (3.16) up to $t=0.05$ . The maximum absolute error in Figure 2b is 0.0033, and occurs at $\theta_{i}=1.4556,1.6860$ , for $\mu_{2}=10$ .

\begin{overpic}[width=212.47617pt]{Figures/Section3/Mu1AnalyticalVsNumerical.eps} \put(3.0,73.0){a)} \put(40.0,15.0){\vector(0,1){50.0}} \put(28.0,66.0){$\mu_{1}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section3/Mu2AnalyticalVsNumerical.eps} \put(3.0,73.0){b)} \put(62.0,42.0){\vector(-1,1){15.0}} \put(28.0,58.0){$\mu_{2}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section3/AbsoluteErrorMu1Thesis300.eps} \put(3.0,73.0){c)} \put(40.0,20.0){\vector(1,2){10.0}} \put(40.0,42.0){$\mu_{1}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section3/AbsoluteErrorMu2Thesis300.eps} \put(3.0,73.0){d)} \put(30.0,30.0){\vector(0,1){38.0}} \put(15.0,70.0){$\mu_{2}$ increasing} \end{overpic}

Figure 2: Comparison of the analytical result for

\hat{\theta}

t^{\prime}=0.05

, with

\theta_{i}\in[0,\pi]

for a)

\mu_{1}=0,1,10

(blue, red, yellow respectively) with

\mu_{2}=\mu_{3}=1

fixed and b)

\mu_{2}=0,1,10

(blue, red, yellow respectively) with

\mu_{1}=\mu_{3}=1

fixed. The analytical results from (3.16) are solid lines, the dotted lines are numerical results obtained by solving (3.9)–(3.12). The absolute errors between the numerical and analytical results are given for c) varied

\mu_{1}

and d) varied

\mu_{2}

3.3 Employing ALE to the short timescale

We again define the reference domain and mapping from the reference to the Eulerian domain by (3.5) and (3.6). In Section 2.2 we determined that on the short timescale there is no extension or thinning of the sheet, which gives the simpler mesh velocity

\displaystyle\mathbf{u}_{\text{mesh}}=\left(0,\frac{\partial H}{\partial t}(x,t)\right).

(3.17)

The fibre director equation on the short timescale, (2.20), becomes

\displaystyle\frac{\partial\theta}{\partial\tau}=0,

(3.18)

hence there is no evolution in $\theta$ in the ALE framework on this timescale. Next, equation (2.28) gives,

\frac{\partial}{\partial x^{\prime}}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\dfrac{\mu_{1}\cos 2\theta+\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial\bar{u}}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}y^{\prime}+\frac{1}{L}\frac{\partial}{\partial x^{\prime}}\left(\frac{\partial H}{\partial x^{\prime}}\dfrac{\partial^{2}H}{\partial x^{\prime}\partial\tau}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\dfrac{\left(4+4\mu_{3}+\mu_{2}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}y^{\prime}\right.\\ \left.-\dfrac{\partial^{3}H}{\partial x^{\prime 2}\partial\tau}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\dfrac{\left(4+4\mu_{3}+\mu_{2}\right)y^{\prime}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}y^{\prime}\right)=0,

(3.19)

with the momentum equation (2.29) becoming

\frac{\partial^{2}}{\partial x^{\prime 2}}\left(\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{\tilde{y}}\dfrac{\mu_{1}\cos 2\theta+\left(4+4\mu_{3}+\mu_{2}\right)\frac{1}{L}\frac{\partial\bar{u}}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}\tilde{y}^{\prime}\mathrm{d}\tilde{y}\right.\\ \left.+\frac{\partial H}{\partial x^{\prime}}\frac{\partial^{2}H}{\partial x^{\prime}\partial\tau}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{\tilde{y}}\dfrac{4+4\mu_{3}+\mu_{2}}{4+\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}\tilde{y}^{\prime}\mathrm{d}\tilde{y}\right.\left.-\frac{\partial^{3}H}{\partial x^{\prime 2}\partial\tau}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{\tilde{y}}\dfrac{(4+4\mu_{3}+\mu_{2})\tilde{y}^{\prime}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{3}\mathrm{d}\tilde{y}^{\prime}\mathrm{d}\tilde{y}\right)\\ =\left(\frac{\partial^{2}H}{\partial x^{\prime 2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{\prime 2}}\right)\left(\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\dfrac{\mu_{1}\cos 2\theta+\left(4+4\mu_{3}+\mu_{2}\right)\frac{1}{L}\frac{\partial\bar{u}}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}\tilde{y}\right.\\ \left.+\frac{\partial H}{\partial x^{\prime}}\frac{\partial^{2}H}{\partial x^{\prime}\partial\tau}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\dfrac{4+4\mu_{3}+\mu_{2}}{4+4\mu_{2}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}\tilde{y}\right.+\left.-\frac{\partial^{3}H}{\partial x^{\prime 2}\partial\tau}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\dfrac{\left(4+4\mu_{3}+\mu_{2}\right)\tilde{y}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}\tilde{y}\right).

(3.20)

As the equation for $\theta$ , (3.18), tells us that there is no evolution of $\theta$ in the ALE framework, we now need only to solve (3.19), (3.20) for $\bar{u}$ and $H$ respectively. As $\dfrac{\partial\bar{u}}{\partial x}$ depends on only $x^{\prime},\tau$ , it can be removed from the integrals, the integral coefficients are pre–computable at each time step, which significantly eases the numerical implementation.

The strategy is as follows: given an initial condition for $\theta,H$ , we simultaneously solve for $\bar{u}$ at the current time step, and $H$ at the next time step using a forward time centred space finite difference discretisation. We repeat this until we reach the desired time. We give the discretisations of equations (2.20),(3.19),(3.20) and details of the construction of the resulting linear system in Appendix E.

4 Results on the extensional flow timescale

We first examine the behaviours of the sheet on the $t\sim\frac{L_{0}}{U}$ timescale. Green & Friedman considered special cases where further analytical insight into the behaviour of the model is possible . They considered two cases: where the sheet is not extending $(\dot{L}=0)$ , in which the fibres tend to align parallel to the $y$ -axis, and the case of an extensional flow with $\mu_{1}=\mu_{2}=0$ . In this case, the equations imply that $H\equiv 0$ , and are amenable to further progress via Lagrangian transformation. It is found that the fibres align with the direction of extension (i.e parallel to the $x$ -axis) [12]. Throughout the next section we assume that $\mu_{1}=0$ in order to examine the contributions of the anisotropic viscosities $\mu_{2},\mu_{3}$ to the behaviour of the sheet.

4.1 Solutions for initially-uniform transversely-isotropic sheets with $\mu_{1}=0$

We now aim to study the effect of the anisotropic viscosities $\mu_{2},\mu_{3}$ upon the fluid as the sheet is stretched. We first turn our attention to transversely isotropic sheets that have an initially constant thickness and $\mu_{1}=0$ . For clarity of exposition, we introduce the quantity

\displaystyle G\left(x^{\prime},y^{\prime},t\right)=hG_{1}+\frac{h}{L}\frac{\partial u}{\partial x^{\prime}}G_{2}=\int\limits_{-\frac{1}{2}}^{y^{\prime}}\frac{\mu_{1}\cos 2\theta+\left(4+4\mu_{3}+\mu_{2}\right)\frac{1}{L}\frac{\partial u}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}y^{\prime},

(4.1)

so that

	$\displaystyle G_{1}(x^{\prime},y^{\prime},t)=\int\limits_{-\frac{1}{2}}^{y^{\prime}}\frac{\mu_{1}\cos 2\theta}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}\mathrm{d}y^{\prime},$		(4.2)
	$\displaystyle G_{2}(x^{\prime},y^{\prime},t)=\int\limits_{-\frac{1}{2}}^{y^{\prime}}\frac{\left(4+4\mu_{3}+\mu_{2}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}\mathrm{d}y^{\prime}.$		(4.3)

The definitions of $G_{1}$ and $G_{2}$ are readily obtained in the physical domain by performing the inverse of the transformation $\boldsymbol{\Phi}$ and reintroducing dimensional parameters . We start by considering sheets with no tension in the fibre direction when the fluid is at rest, that is, $\mu_{1}=0$ , in which case the equation for $u$ in ALE variables, (3.11), becomes

\displaystyle\frac{\partial}{\partial x^{\prime}}\left(h\frac{\partial u}{\partial x^{\prime}}G_{2}\left(x^{\prime},\frac{1}{2},t\right)\right)=0.

(4.4)

This is of the same form as the longitudinal momentum equation in the Trouton model, with $G_{2}$ playing the role of a spatially-varying viscosity (setting $\mu_{2}=\mu_{3}=0$ gives $G_{2}=4$ , the Trouton ratio for a Newtonian thin sheet [16]). We interpret $G_{2}$ as a heterogeneous, time-dependent, ‘effective viscosity’. As we shall show, we see the effect of $G_{2}$ is to induce ‘necking’ in the sheet (the sheet undergoes thinning at a greater rate in some areas of the sheet than others, generating a ‘neck’). This behaviour has been observed in Newtonian fluids which possess an inhomogenous viscosity.

We also note that if the fibre angle is independent of $x^{\prime}$ , or if $\mu_{1}=\mu_{2}=0$ [12], then $G_{2}$ as it appears in equation (4.4) does not possess $x^{\prime}$ dependence. Hence, a fluid that does not possess tension in the fibre direction when at rest, and has a uniform fibre direction will behave like a Newtonian fluid, with a modified viscosity. In this case, the centre–line of the fluid is always a straight line, and the model can be solved by transforming to Lagrangian variables, as detailed in [12, 16].

A first integral of (4.4) yields,

\displaystyle h\frac{\partial u}{\partial x^{\prime}}G_{2}\left(x^{\prime},\frac{1}{2},t\right)=T\left(t\right),

(4.5)

where $T$ is the tension applied to the sheet. Directly from (4.5) we see that $G_{2}$ may induce non-linear behaviours in $u$ . For a Newtonian fluid, or a transversely isotropic fluid where only $\mu_{3}\neq 0$ , choosing an initial condition of constant thickness across the sheet yields that $u$ must be linear in $x^{\prime}$ . In that case, as shown by Howell, $u=\dot{L}x$ for all time, and $h$ is a function of $t$ only [16] (using a similar approach, the same result was shown for a fluid with $\mu_{1}=\mu_{2}=0,$ [12]). However, for a transversely isotropic sheet with $\mu_{2}\neq 0$ , the existence of the trigonometric terms inside $G_{2}$ prohibits this. When $\mu_{2}>0$ , with $h\left(x^{\prime},0\right)$ constant, if $\theta(x^{\prime},y^{\prime},0)$ depends upon $x^{\prime}$ , this will result in $u$ becoming nonlinear in $x^{\prime}$ , and by (3.9), this will cause $h$ to become spatially non-uniform. In Figure 3, we illustrate the early evolution of the thickness of the sheet and the behaviour of ${1}/{G_{2}}$ , for $h\left(x^{\prime},0\right)=1,\theta\left(x^{\prime},y^{\prime},0\right)=\cos\left(4\pi x^{\prime}y^{\prime}\right)-0.1,\mu_{1}=\mu_{3}=0,\mu_{2}=5$ . Notice that the sheet thickness immediately becomes non-uniform, and the location of the peaks and troughs in $1/G_{2}$ correlates with the locations of local minima and maxima of the thickness of the sheet when $t\neq 0$ .

\begin{overpic}[width=212.47617pt]{Figures/Section4/ShortEvolution/shorttimeshapeofsheet2NL.eps} \put(3.0,73.0){a)} \put(50.0,9.0){\vector(0,1){15.0}} \put(35.0,25.0){$t$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/ShortEvolution/shorttimeinverseG22NL.eps} \put(3.0,73.0){b)} \put(52.0,20.0){\vector(0,1){22.0}} \put(38.0,44.0){$t$ increasing} \end{overpic}

Figure 3: The evolution of a) the thickness of the sheet and b) the function

\dfrac{1}{G_{2}}

at the times

t=0,0.125,0.25,0.375,0.5

, pulling at

L(t)=1+t

, with the initial conditions

h\left(x^{\prime},0\right)=1,\theta\left(x^{\prime},y^{\prime},0\right)=\cos\left(4\pi x^{\prime}y^{\prime}\right)-0.1,\mu_{1}=\mu_{3}=0,\mu_{2}=5

4.1.1 Effect of varying the extensional and shear viscosities $\mu_{2},\mu_{3}$ with $\mu_{1}=0.$

In this section, we continue to use the initial conditions $h\left(x^{\prime},0\right)=1,\>\theta\left(x^{\prime},y^{\prime},0\right)=\cos\left(4\pi x^{\prime}y^{\prime}\right)-0.1,\>\mu_{1}=0$ , but now vary $\mu_{2}$ and $\mu_{3}$ and compare the state of the sheet at $t=5$ . First, we note that varying $\mu_{3}$ with $\mu_{2}=0$ simply changes the value of the constant obtained from $G_{2}$ . We plot in Figure 4 the thickness and velocity $u$ in the sheet for varied $\mu_{2}$ and notice that for these choices of $\theta(x^{\prime},y^{\prime},0)$ and $h(x^{\prime},0)$ that we see a global increase in the longitudinal velocity for increasing $\mu_{2}$ . Additionally, we see that there are regions of the sheet than thin more quickly for increasing $\mu_{2}$ , and other regions that thin more slowly. Intuition based upon the behaviour of a pipe flow would lead one to expect that in regions where the sheet is thicker, the velocity would be lower. This is not true here, and as in the previous subsection, this behaviour is linked to the behaviour of $G_{2}$ as we shall now demonstrate. Integrating (4.5) and using that $u(L)=1$ , we may write

1=T\int_{0}^{L}\dfrac{1}{hG_{2}(x^{\prime},\frac{1}{2},t)}\mathrm{d}x^{\prime}.

(4.6)

We note that by (4.3) and (4.6), as the fibres within the fluid sheet flatten out and $G_{2}$ increases everywhere and hence so does the tension, $T$ , and that the tension will increase with increasing $\mu_{2}$ . Eliminating tension from equation (4.6) gives a result for $u$ ,

\displaystyle u=\dfrac{\int\limits_{0}^{x^{\prime}}\dfrac{1}{hG_{2}(s,\frac{1}{2},t)}\mathrm{d}s}{\int\limits_{0}^{L}\dfrac{1}{hG_{2}(x^{\prime},\frac{1}{2},t)}\mathrm{d}x^{\prime}},

(4.7)

and so $u$ behaves as a cumulative integral. If $G_{2}$ is fixed, the longitudinal velocity does behave as a pipe flow intuition would expect (in regions where $h$ is small, $1/h$ is large, and $u$ will increase). However, when $G_{2}$ is not fixed (i.e $\mu_{2}\neq 0$ ), the behaviour of the longitudinal velocity is more complex. The derivative $u_{x}$ is dependent on the behaviour of $hG_{2}$ . We see this in Figure 4, where we give results for $\mu_{3}=0$ and see that for increasing $\mu_{2}$ , $hG_{2}$ decreases (despite the increase in $G_{2}$ ), which leads to an greater $u_{x}$ on the left hand side of the domain. Where $hG_{2}$ is larger, around $x=3$ , we see the gradient of $u$ decrease below that of $\mu_{2}=0$ .

If we allow $\mu_{2},\mu_{3}\neq 0$ , we find that increasing $\mu_{3}$ has the effect of globally increasing the value of $G_{2}$ and hence the tension applied to the ends of the sheet, and inhibits the uniformity-breaking behaviour of $\mu_{2}$ . That is, increasing the $\mu_{3}$ term drives the fluid to behave in a ‘more Newtonian’ manner, albeit with a higher tension, whilst increasing $\mu_{2}$ drives the non-Newtonian behaviour of breaking the uniformity of the sheet previously discussed. As an illustrative example, in Figure 5 we plot the behaviour of the sheet at $t=5$ for varying $\mu_{3}$ with $\mu_{2}=5$ and initial conditions $h(x^{\prime},0)=1,\>\theta(x^{\prime},y^{\prime},0)=\cos\left(4\pi x^{\prime}y^{\prime}\right)-0.1$ . Notice that increasing $\mu_{3}$ causes the thickness of the sheet to exhibit less deviation from uniformity, and the longitudinal velocity to tend towards $u=x^{\prime}$ , the solution for a Newtonian fluid.

\begin{overpic}[width=212.47617pt]{Figures/Section4/VariedMu2/VariedThicknessM2ChapterNL.eps} \put(3.0,73.0){a)} \put(50.0,35.0){\vector(0,1){30.0}} \put(48.0,66.0){$\mu_{2}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/VariedMu2/VariedVelocityM2ChapterNL.eps} \put(3.0,73.0){b)} \put(52.0,30.0){\vector(0,1){35.0}} \put(50.0,67.0){$\mu_{2}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/VariedMu2/VariedMu2G2ChapterNL.eps} \put(3.0,73.0){c)} \put(70.0,10.0){\vector(0,1){55.0}} \put(68.0,66.0){$\mu_{2}$ increasing} \end{overpic}

Figure 4: Comparison of a) thickness, b) longitudinal velocity and c)

G_{2}

t=5

, pulling with

L(t)=1+t

, for

\mu_{2}=0,1,2,3,4,5,10,15

, with the conditions

h\left(x^{\prime},0\right)=1,\>\theta\left(x^{\prime},y^{\prime},0\right)=\cos\left(4\pi x^{\prime}y^{\prime}\right)-0.1,\>\mu_{1}=\mu_{3}=0

. Notice that more extreme behaviour in

G_{2}

correlates with greater change in the thickness and velocity profile across the sheet and increases with

\mu_{2}

\begin{overpic}[width=212.47617pt]{Figures/Section4/VariedMu3/VariedMu3ThicknessNL.eps} \put(3.0,73.0){a)} \put(50.0,60.0){\vector(0,-1){25.0}} \put(47.5,32.0){$\mu_{3}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/VariedMu3/VariedMu3VelocityNL.eps} \put(3.0,73.0){b)} \put(50.0,50.0){\vector(0,-1){15.0}} \put(47.5,32.0){$\mu_{3}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/VariedMu3/VariedMu3G22NL.eps} \put(3.0,73.0){c)} \put(50.0,65.0){\vector(0,-1){44.0}} \put(47.5,18.0){$\mu_{3}$ increasing} \end{overpic}

Figure 5: Comparison of a) thickness of the sheet, and b) longitudinal velocity of the sheet at

t=5

and c)

G_{2}

t=5

\theta(x^{\prime},y^{\prime},0)=\cos\left(4\pi x^{\prime}y^{\prime}\right)-0.1

, for

\mu_{1}=0,\>\mu_{2}=5

, and

\mu_{3}=0,1,2,3,4,5

4.2 Behaviour of the fibres

4.2.1 Uniform initial fibre director field

We now turn our attention to the fibre behaviour within the sheet. In [12], Green & Friedman considered two special cases that validate our results. It is found that for the case of an extensional flow with $\mu_{1},\mu_{2}=0$ the fibres align with the direction of extension (i.e parallel to the $x$ -axis) and that for the special case of not extending the sheet ( $\dot{L}=0$ ) with $\mu_{1},\mu_{2},\mu_{3}\neq 0$ , the fibres tend to align parallel to the $y$ -axis.

For $\mu_{1}=0$ , we start by extending the results obtained for early-time with a constant initial fibre direction described in Section 6 of Green & Friedman [12]. If $h_{i},\theta_{i}$ are uniform, then we can see from equations (3.10) and (3.11) that $\dfrac{\partial u}{\partial x^{\prime}}=1$ , and $\theta$ and $h$ must be functions of time only. Moreover, (3.10) yields

\frac{\mathrm{d}\theta}{\mathrm{d}t}=\frac{-\sin 2\theta\left(4+4\mu_{3}+2\mu_{2}\sin^{2}\theta\right)\frac{1}{L}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta},

(4.8)

which implies that the behaviour of the fibres is determined by the sign of $\sin 2\theta$ . For $0<\theta<\dfrac{\pi}{2}$ , we have $\dfrac{\partial\theta}{\partial t}<0$ , whilst for $\dfrac{\pi}{2}<\theta<\pi$ , $\dfrac{\partial\theta}{\partial t}>0$ . Therefore, much like the case $\mu_{1}=\mu_{2}=0$ studied in [12], the fibres tend to orient themselves with the direction of extension, regardless of the value of $\mu_{2}$ , given uniform initial conditions for $h$ and $\theta$ .
Additionally, we note that $\theta=\pm\frac{\pi}{2}$ are unstable fixed points of equation (4.8), whilst $\theta=0,\pi$ are stable fixed points.

We now consider the case $\mu_{1}\neq 0$ , for an initially constant $h_{i}$ and $\theta_{i}$ . We again have that $\dfrac{\partial u}{\partial x^{\prime}}=1$ and $h,\theta$ remain uniform for all time. However, equation (3.10) now yields

\displaystyle\frac{\mathrm{d}\theta}{\mathrm{d}t}=\frac{\sin 2\theta\left(2\mu_{1}\sin^{2}\theta-\left(4+4\mu_{3}+2\mu_{2}\sin^{2}\theta\right)\frac{1}{L}\right)}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta},

(4.9)

and so the evolution of the fibre direction is less clear. By again considering the cases $0<\theta<\dfrac{\pi}{2}$ , $\dfrac{\pi}{2}<\theta<\pi$ , $-\dfrac{\pi}{2}<\theta<0$ , $-\pi<\theta<-\dfrac{\pi}{2}$ we find that in order for the fibres to align along the $x$ axis with time, we require

\displaystyle\sin^{2}\theta\left(L\mu_{1}-\mu_{2}\right)<2+2\mu_{3}.

(4.10)

If this condition is satisfied, fibres with angles between $-\dfrac{\pi}{2}<\theta<\dfrac{\pi}{2}$ will rotate towards the positive $x$ -axis, with fibres outside of this range rotating towards the negative $x$ -axis, similarly to the above. However, since this expression includes $L(t)$ , it is possible for fibres that initially rotate towards the longitudinal orientation can reverse their evolution as $L\mu_{1}$ grows with time to violate (4.10). We illustrate this behaviour in Figure 6. For the choices of $\theta_{i}=\frac{\pi}{4},L_{i}=h_{i}=1,\mu_{1}=5,\mu_{2}=0,\mu_{3}=1$ , we give the evolution of the fibre angle and evolution of the left hand side (4.10).

\begin{overpic}[width=212.47617pt]{Figures/Section4/RHSUniformThetaTheta.eps} \put(3.0,73.0){a)} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/RHSUniformTheta.eps} \put(3.0,73.0){b)} \end{overpic}

Figure 6: Evolution of a) the fibre angle

\theta

, and b) the left hand side of (4.10) (blue), with the threshold

2+2\mu_{3},

(red), for the choices of

\theta_{i}=\frac{\pi}{4},L_{i}=h_{i}=1,\mu_{1}=5,\mu_{2}=0,\mu_{3}=1

. Note the reversal in the direction of rotation when

L\mu_{1}\sin^{2}\theta\geq 2+2\mu_{3}

4.2.2 Non-uniform initial fibre angles

Consider now the extension of the sheet with $\theta(x^{\prime},y^{\prime},0)=\cos(4\pi x^{\prime}y^{\prime})-0.1$ , $h(x^{\prime},0)=1$ and $L(t)=1+t$ . We see in Figure 7 that the fibres have a tendency to align in the direction of the sheet when $\mu_{1}=0$ . We also note that the rate at which the fibres align in this direction is enhanced in the regions of the sheet that undergoes fastest thinning of the fluid.

\begin{overpic}[width=212.47617pt]{Figures/Section4/Quiver/ThesisInitialQuiver.eps} \put(0.0,75.0){a)} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/Quiver/ThesisFinalQuiverBLU2.eps} \put(3.0,75.0){b)} \end{overpic}

Figure 7: a) The initial thickness of the sheet and orientation of the fibres, and b) the thickness of the sheet and orientation of the fibres at

t=5

, for the initial conditions

h(x^{\prime},0)=1,\theta(x^{\prime},y^{\prime},0)=\cos(4\pi x^{\prime}y^{\prime})-0.1

and choices of

\mu_{1}=\mu_{3}=0,\mu_{2}=5.

4.3 Special cases for sheets possessing active behaviour

Here, we briefly examine two special cases for the sheet for which $\mu_{1}>0$ . We start by prescribing $L(t)=1+t$ , $h_{i}=1$ , $\mu_{1}=5$ , $\mu_{2}=\mu_{3}=0$ , and we test two choices for the initial fibre angle, $\theta_{i}=0$ , corresponding to the fibres being aligned in the direction of extension, and $\theta_{i}=\dfrac{\pi}{2}$ , corresponding to the fibres being aligned in the transverse direction in the sheet. In Figure 8 we plot the tension required to be applied to the sheet to achieve the prescribed rate of pulling for these two choices. First, we note that there is no evolution in $\theta$ with time. Next, we notice that for fibres arranged in the transverse direction of the sheet, that the tension applied to the sheet is negative. To explain this we begin with the equation for $u$ in spatial variables,

\displaystyle T=\int_{H^{-}}^{H^{+}}\frac{\mu_{1}\cos 2\theta+\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial u}{\partial x}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}\mathrm{d}y,

(4.11)

now since $\theta$ has no evolution for these particular choices of $\theta_{i}$ (as they are fixed points of equation (3.8)), then for $\theta_{i}=\dfrac{\pi}{2}$ ,

\displaystyle T=\frac{1}{4\left(1+\mu_{3}\right)}\left(\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial u}{\partial x}-\mu_{1}\right)h.

(4.12)

It is possible that $\mu_{1}>\left(4+4\mu_{3}+\mu_{2}\right)u_{x}$ and hence $T<0$ . We interpret this as being due to the fibres in the transverse direction attempting to contract the sheet, which due to mass conservation, would generate a compression in the longitudinal direction, as the sheet attempts to extend longitudinally. If we define the total tension required to move the sheet at the prescribed speed as $T_{L}$ , the tension caused by the fibres as $T_{f}$ and the tension applied to the sheet as $T$ , then we expect that

T_{L}=T+T_{f}.

(4.13)

If the rate of extension is too slow to compensate for the compression generated by the fibres pulling in the transverse direction, then $T<0$ . As discussed by Howell for the Newtonian case, [16], when the sheet is in compression we expect buckling to occur, and that the curvature of the centre–line will in time become significant. Thus, in this case, the nearly-straight centre–line scaling assumed in the Green & Friedman model may be violated, and the model will no longer be valid.

\begin{overpic}[width=212.47617pt]{Figures/Section4/Tension/TensionTheta0Chapter.eps} \put(3.0,73.0){a)} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section4/Tension/TensionPiBy2.eps} \put(3.0,73.0){b)} \end{overpic}

Figure 8: Evolution of the tension applied to the sheet for two initial choices of fibre direction, a)

\theta_{i}=0

, and b)

\theta_{i}=\dfrac{\pi}{2}

up to

t=5

with

h(x^{\prime},0)=1,\>\mu_{1}=5,\>\mu_{2}=\mu_{3}=0

, and prescribed pulling

L=1+t

4.4 Behaviour of the centre–line

4.4.1 Conditions for a flat centre–line

In the Newtonian problem, it was shown by Howell [16] that the centre–line of the sheet straightens on a timescale shorter than $\dfrac{L_{0}}{U}$ [2, 16] (and is therefore generally taken to simply be $H=0$ ). This is not necessarily true for a transversely-isotropic fluid. It is noted in [12] that should $\theta=\theta(x^{\prime},t)$ and $\mu_{1}=\mu_{2}=0$ , then equations (3.11) and (3.12) together with the requirement that $H(0,t)=H(L(t),t)=0$ , imply that the centre–line must be flat. We will now demonstrate that this is also true when $\mu_{1},\mu_{2}\neq 0$ . Supposing $\theta=\theta(x^{\prime},t)$ , then the integrand of equations (3.11) and (3.12) can be evaluated explicitly, yielding

	$\displaystyle\frac{\partial}{\partial x^{\prime}}\left(hf\right)$	$\displaystyle=0,$		(4.14)
	$\displaystyle\frac{\partial^{2}}{\partial x^{\prime 2}}\left(\frac{1}{2}h^{2}f\right)$	$\displaystyle=\left(\frac{\partial^{2}H}{\partial x^{\prime 2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{\prime 2}}\right)hf,$		(4.15)

where

f(x^{\prime},t)=\frac{\mu_{1}\cos 2\theta+\frac{1}{L}\left(4+4\mu_{3}+\mu_{2}\right)\frac{\partial u}{\partial x^{\prime}}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}.

(4.16)

Expanding the second derivative on the LHS of equation (4.15) and using (4.14) yields that

\displaystyle\frac{\partial^{2}H}{\partial x^{\prime 2}}=0,

(4.17)

and therefore $H$ must be a straight line. It is then possible to choose our co-ordinate system such that $H(x^{\prime},t)\equiv 0$ .

Now, let us consider what happens when $\theta(x^{\prime},y^{\prime},0)=\theta_{i}(y^{\prime})$ , and $\mu_{1}=0$ . In this case the equation of the centre–line can be written as

\displaystyle\left(\frac{\partial^{2}H}{\partial x^{\prime 2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{\prime 2}}\right)hG_{2}\left(\frac{1}{2},t\right)=\frac{\partial^{2}}{\partial x^{\prime 2}}\left(h^{2}\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}G_{2}(y^{\prime},t)\mathrm{d}y^{\prime}\right).

(4.18)

We recall that, from Section 4.1, if we choose $\theta_{i}=\theta_{i}(y^{\prime})$ only and the initial thickness $h_{i}$ to be uniform, then $G_{2}$ will not possess $x^{\prime}$ -dependence and $h$ will also remain uniform for all time. As a result we once again obtain (4.17), so that asymmetry in the fibre angles over the centre–line is not sufficient to cause an initially uniform sheet with $\mu_{1}=0$ to deflect. We will now consider the conditions under which the centre–line of the sheet will not be straight.

4.4.2 Conditions for a nonzero centre–line

For sheets possessing $\mu_{1}=0$ and a uniform condition for the thickness of the sheet, $h_{i}$ , the choice of $\theta(x^{\prime},y^{\prime},0)=\theta_{i}(x^{\prime},y^{\prime})$ being a prescribed function that is non-symmetric (i.e. $\frac{\partial\theta}{\partial y^{\prime}}|_{y^{\prime}=0}\neq 0$ ) is required to obtain centre–line deflection. In Figure 9 we give a plot of the centre–line evolution for the initial conditions $\theta(x^{\prime},y^{\prime},0)=\sin\left(4\pi x^{\prime}y^{\prime}\right)-0.1$ with $\mu_{2}=5,\>\mu_{1}=\mu_{3}=0$ . We see that the centre–line initially possesses deflection and tends to $0$ . This behaviour is driven by the right hand side of (4.18) becoming small. Once this occurs, the behaviour of the centre–line is then dominated by the $\frac{\partial^{2}h}{\partial x^{\prime 2}}$ term, which itself possesses an implicit dependence on $G_{2}$ through $u$ . Additionally, for a transversely isotropic fluid with $\mu_{1}=0$ , the deflection is small and decays quickly. Since the centre–line does not instantaneously collapse to $0$ , this suggests that considering the behaviour of the fluid on a short timescale may yield some interesting behaviour that is markedly different from a Newtonian fluid.

If we consider a condition for $h(x^{\prime},0)$ that is not uniform, we find that there will exist a small deflection when $\theta(x^{\prime},y^{\prime},0)=\theta_{i}(y^{\prime})$ only, since endowing $h$ with $x^{\prime}$ -dependence will impart $x^{\prime}$ -dependence on $\frac{\partial u}{\partial x}$ through equation (3.11), and hence $\theta$ will gain $x^{\prime}$ -dependence after the first time step through equation (3.10).

\begin{overpic}[width=212.47617pt]{Figures/Section4/Centreline/Centreline800ChapterNL.eps} \put(3.0,73.0){a)} \put(19.0,17.0){\vector(0,1){50.0}} \put(18.0,68.0){$t$ increasing} \end{overpic}

Figure 9: Evolution of the centre–line of the sheet under the initial conditions

h(x^{\prime},0)=1,\>\theta(x^{\prime},y^{\prime},0)=\sin\left(4\pi x^{\prime}y^{\prime}\right)-0.1

with

\mu_{2}=5,\>\mu_{1}=\mu_{3}=0

. Plotted at

t=0,0.25,0.5,0.75,1

5 Results on the short timescale

We now examine the behaviours of the sheet over a short timescale of $t\sim\varepsilon^{2}\dfrac{L_{0}}{U}$ . For a Newtonian fluid, it was demonstrated by Howell [16] that the leading order equations for the extensional flow of a slender, viscous, Newtonian sheet predict that the centre–line of the sheet is straight. Hence, the model cannot satisfy an initial condition in which the centre–line is not straight. In order to study the behaviour of initially curved sheets, a short timescale analysis is performed. The required timescale must be $\varepsilon^{2}\frac{L_{0}}{U}$ [2]. In the transversely isotropic problem it is similarly impossible to satisfy an arbitrary initial condition for the centre–line and, as we have demonstrated, there exist choices of the key parameters and initial distribution of fibre angles that give rise to a centre–line that is non-zero on the flow timescale. This indicates that there may exist interesting behaviours, different from the Newtonian case, over a short timescale.

5.1 Short-time evolution of the centre–line

First, we check that the long-time behaviour of the short time model is consistent with the Green & Friedman model. In Figure 10, we give the evolution of the centre–line over the short timescale for the initial conditions of $\mu_{1}=\mu_{3}=0,\>\mu_{2}=5$ , $L(0)=1$ , $\theta(x,y,0)=\sin(4\pi xy)-0.1$ , ${H^{(0)}(0,\tau)=H^{(0)}(L,\tau)=\dfrac{\partial H^{(0)}}{\partial x}(L,\tau)=\dfrac{\partial H^{(0)}}{\partial x}H(L,\tau)=0}$ , up to $\tau=200000$ , and compare the values of $H(x,\tau=200000)$ with $H(x,t=0)$ . We see that the short time-scale result closely matches the result produced by the solver for the Green & Friedman model. In Figure 10b we give the evolution of the centre–line over $\tau$ . Notice in Figure 10b the centre–line converges to the Green & Friedman model fairly quickly, and there is a small absolute difference between the centre–line at $\tau=10000$ , and $\tau=200000$ .

For a Newtonian fluid, Howell was able to obtain a analytical expression for the decay of the centre–line of an initially curved sheet undergoing stretching by use of eigenfunction expansions [16], assuming the same boundary conditions we use in this section, $H(0,\tau)=H(1,\tau)=\dfrac{\partial H}{\partial x}(0,\tau)=\dfrac{\partial H}{\partial x}(1,\tau)=0$ . It is found that the centre–line decays exponentially to $H=0$ . The behaviour for a transversely isotropic fluid is more complex. The results plotted in Figure 10 began with the initial condition $H(x,\tau=0)=0$ , and we immediately see from Figure 10c that the convergence to the Green & Friedman centre–line is not exponential for all time. Indeed, there is an initial lag, as the centre–line evolves away from the flat initial condition to adopt a similar shape of the centre–line predicted by the Green & Friedman model, before decaying to the expected result. As an illustrative example, we include Figure 10d. This figure shows the formation of the peaks and troughs of the general shape of the centre–line given by the Green & Friedman model.

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/Comparison.eps} \put(3.0,73.0){a)} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/EvolFig.eps} \put(3.0,73.0){b)} \put(52.5,55.0){\vector(0,-1){25.0}} \put(51.5,56.0){$\tau$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/RelDiff.eps} \put(3.0,73.0){c)} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/EarlyEvolFig.eps} \put(3.0,73.0){d)} \put(40.0,30.0){\vector(0,1){40.0}} \put(39.0,27.5){$\tau$ increasing} \end{overpic}

Figure 10: Position of the centre–line a) at

t=0

(red) and

\tau=200000

(blue) and b) at

\tau=100,500,1000,2500,10000,200000

, where

\tau=0

is omitted as

H(x,0)=0

, In c) we plot the maximum of the absolute difference between the evolving centre–line on the short timescale and the Green & Friedman result and in d) we give the centre–line at

\tau=2.5,5,10,15,25

. The conditions for all of the above are

H(x,0)=0,\>h(x,0)=1,\>\theta(x,y,0)=\sin(4\pi xy)-0.1,\>\mu_{1}=\mu_{3}=0,\>\mu_{2}=5

L(0)=1

5.2 Effect of key parameters upon convergence

In this subsection we discuss the effect of increasing $\mu_{2},\mu_{3}$ upon the convergence of the short timescale centre–line to the result from the Green & Friedman model. First we note that if $\mu_{2}=0$ , $\mu_{3}$ has no effect upon convergence. In Figure 11 we give a plot of the decay of the maximum value of $H(x,\tau)$ for the initial conditions of $H(x,0)=\frac{1}{20}e^{-\frac{\left(x-0.5\right)^{2}}{0.01}},\theta(x,y,0)=\sin(4\pi xy),h(x,0)=1,L=1$ . There is no difference in the decay of the centre–line to flat between the Newtonian case and $\mu_{3}=10$ .

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/Convergence/Mu3GaussianDecayNL.eps} \end{overpic}

Figure 11: The decay of the maximum value of

H(x,\tau)

for

\mu_{3}=0

(blue) and

\mu_{3}=10

(red) with the conditions

H(x,0)=\frac{1}{20}e^{-\frac{\left(x-0.5\right)^{2}}{0.01}},\>\mu_{1}=\mu_{2}=0,\>h(x,0)=1,\>L=1

Now choosing $\mu_{2}>0$ , in Figure 12 we plot the maximum absolute difference across the sheet between the centre–line on the short timescale and the result obtained by solving equation (3.12) for varied values of $\mu_{2},\mu_{3}$ with the conditions $h(x,0)=1,\theta(x,y,0)=\sin(4\pi xy)-0.1,\mu_{1}=0,\mu_{2}=5,L(0)=1$ . First, in Figure 12a, we fix $\mu_{3}=0$ and vary $\mu_{2}$ (note the case of $\mu_{2}=5$ corresponds to the example given above). We see that as $\mu_{2}$ increases, the initial difference between the flat initial condition of the centre–line and the result of the Green & Friedman model also increases. This is due to the deepening of the trough around $x=0.5$ seen in Figure 10a. We note that changing $\mu_{2}$ does not appear to affect the length of the lag as the centre–line is adopting the correct general shape before decaying, but as $\mu_{2}$ increases, the decay is faster. In Figure 12b, we fix $\mu_{2}=5$ and vary $\mu_{3}$ (as before, the case of $\mu_{3}=0$ corresponds to the example given above). Once again, we see that the effect of increasing $\mu_{3}$ has little affect upon the convergence of the centre–line, and that the behaviour we see here is a consequence of the role of $\mu_{3}$ in moderating the effects of $\mu_{2}$ in the Green & Friedman model as discussed in Section 4.

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/Convergence/ConvergenceFigureMu2.eps} \put(3.0,73.0){a)} \put(20.0,20.0){\vector(0,1){50.0}} \put(18.0,72.0){$\mu_{2}$ increasing} \end{overpic}

\begin{overpic}[width=212.47617pt]{Figures/Section5/BNT/Convergence/ConvergenceComparisonMu3.eps} \put(3.0,73.0){b)} \put(20.0,63.0){\vector(0,-1){46.0}} \put(15.5,14.5){$\mu_{3}$ increasing} \end{overpic}

Figure 12: Maximum absolute difference across the sheet between the centre–line on the short timescale and the result obtained from the Green & Friedman model (3.12) for a)

\mu_{3}=0

\mu_{2}=1,2,3,4,5

and b)

\mu_{2}=5

\mu_{3}=0,1,2,3,4,5

, with the conditions

H(x,0)=0,\>h(x,0)=1,\>\theta(x,y,0)=\sin(4\pi xy)-0.1,\>\mu_{1}=0,\>L(0)=1

6 Discussion

In this paper we have constructed and employed a numerical strategy to solve the model proposed by Green & Friedman for the extensional flow of a thin two-dimensional sheet of a fibre-reinforced fluid, first reducing the model by eliminating $u^{(1)}$ and then employing a arbitrary Lagrangian–Eulerian method. We have shown how the distribution of fibres within the fluid can cause interesting non–Newtonian behaviours such as driving non-uniformity in the development of the thickness of an initially uniform sheet and deflection of the centre–line even with the implicit assumption that the centre–line is nearly straight. Our results also show that the bulk properties of a passive transversely isotropic fluid sheet are controlled largely by the behaviour of a derived ‘effective viscosity’.

As far as the behaviour of an active transversely isotropic fluid is concerned, preliminarily we have seen that allowing $\mu_{1}\neq 0$ allows the fibres to develop towards alignments that are not in the direction of extension of the fluid. However, if the fibres are aligned in the longitudinal direction of the sheet and possesses active behaviour, the tension within the sheet is increased. Active behaviour giving rise to greater tension has been observed in the seeding of hydrogels with a suspension of self-aligning cells [1]. Future work in this area could include constructing more biologically realistic multiphase model that incorporates the work in this paper as a fibrous extracellular matrix or hydrogel, with the cells exhibiting active behaviour instead of the fibres. This could result in a model of how different experimental setups lead to different alignment patterns of cells and could determine the best conditions to grow neural tissue.

There are a number of other avenues for further work related to this paper. We model sheets that are nearly straight, with the employment of a Cartesian co-ordinate system restricting the model to examining sheets which are initially slightly curved, i.e $\dfrac{H}{L_{0}}$ is small. Where this is not the case, future work could entail the use of a curvilinear co-ordinate system to approach sheets with curvature in the centre–line, in works similar to Ribe [25]. As the sheet starts to become very thin, there may be a new regime where the $\mu_{1}$ term in (3.11) dominates the $\dfrac{\partial u}{\partial x^{\prime}}$ term. Perhaps the behaviour of the sheet in this regime may shed light upon how the active behaviour of the fluid may drive breakup of the sheet. Simpler modifications could include prescribing the tension applied to the ends of the sheet, rather than prescribing the length. Furthermore we could also modify the model to include the effects of surface tension, inertia, and body forces.

Declaration of competing interests
None.

Appendix A Model equations in full

To briefly summarise, the dimensionless equations are

\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,

(A.1)

for conservation of mass, with the momentum equation (2.2) yielding in the $x$ -direction:

$\displaystyle-\varepsilon^{2}\frac{\partial p}{\partial x}$	$\displaystyle+\frac{\partial^{2}u}{\partial y^{2}}+\varepsilon^{2}\frac{\partial^{2}u}{\partial x^{2}}+\varepsilon^{2}\mu_{1}\frac{\partial}{\partial x}(\cos^{2}{\theta})+\varepsilon\mu_{1}\frac{\partial}{\partial y}(\cos{\theta}\sin{\theta})$
	$\displaystyle+\mu_{2}\frac{\partial}{\partial x}\left[\varepsilon^{2}\cos^{4}{\theta}\frac{\partial u}{\partial x}+\cos^{3}{\theta}\sin{\theta}\left(\varepsilon\frac{\partial u}{\partial y}+\varepsilon^{3}\frac{\partial v}{\partial x}\right)+\varepsilon^{2}\cos^{2}{\theta}\sin^{2}{\theta}\frac{\partial v}{\partial y}\right]$
	$\displaystyle+2\mu_{3}\frac{\partial}{\partial x}\left[2\varepsilon^{2}\cos^{2}{\theta}\frac{\partial u}{\partial x}+\cos{\theta}\sin{\theta}\left(\varepsilon\frac{\partial u}{\partial y}+\varepsilon^{3}\frac{\partial v}{\partial x}\right)\right]$
	$\displaystyle+\mu_{2}\frac{\partial}{\partial y}\left[\cos{\theta}\sin{\theta}\left(\varepsilon\cos^{2}{\theta}\frac{\partial u}{\partial x}+\cos{\theta}\sin{\theta}\left(\frac{\partial u}{\partial y}+\varepsilon^{2}\frac{\partial v}{\partial x}\right)+\varepsilon\sin^{2}{\theta}\frac{\partial v}{\partial y}\right)\right]$
	$\displaystyle+\mu_{3}\frac{\partial}{\partial y}\left[\frac{\partial u}{\partial y}+\varepsilon^{2}\frac{\partial v}{\partial x}\right]=0,$	(A.2)

whilst in the $y$ direction we have:

$\displaystyle-\varepsilon\frac{\partial p}{\partial y}$	$\displaystyle+\varepsilon^{3}\frac{\partial^{2}v}{\partial x^{2}}+\varepsilon\frac{\partial^{2}v}{\partial y^{2}}+\varepsilon\mu_{1}\frac{\partial}{\partial y}(\sin^{2}{\theta})+\varepsilon^{2}\mu_{1}\frac{\partial}{\partial x}(\cos{\theta}\sin{\theta})$
	$\displaystyle+\mu_{2}\frac{\partial}{\partial y}\left[\varepsilon\sin^{2}{\theta}\cos^{2}{\theta}\frac{\partial u}{\partial x}+\cos{\theta}\sin^{3}{\theta}\left(\frac{\partial u}{\partial y}+\varepsilon^{2}\frac{\partial v}{\partial x}\right)+\varepsilon\sin^{4}{\theta}\frac{\partial v}{\partial y}\right]$
	$\displaystyle+2\mu_{3}\frac{\partial}{\partial y}\left[2\varepsilon\sin^{2}{\theta}\frac{\partial v}{\partial y}+\cos{\theta}\sin{\theta}\left(\frac{\partial u}{\partial y}+\varepsilon^{2}\frac{\partial v}{\partial x}\right)\right]$
	$\displaystyle+\mu_{2}\frac{\partial}{\partial x}\left[\varepsilon^{2}\sin{\theta}\cos^{3}{\theta}\frac{\partial u}{\partial x}+\cos^{2}{\theta}\sin^{2}{\theta}\left(\varepsilon\frac{\partial u}{\partial y}+\varepsilon^{3}\frac{\partial v}{\partial x}\right)+\varepsilon^{2}\cos{\theta}\sin^{3}{\theta}\frac{\partial v}{\partial y}\right]$
	$\displaystyle+\mu_{3}\frac{\partial}{\partial x}\left[\varepsilon\frac{\partial u}{\partial y}+\varepsilon^{3}\frac{\partial v}{\partial x}\right]=0,$	(A.3)

with the fibre director field being given by

\displaystyle\varepsilon\frac{\partial\theta}{\partial t}+\varepsilon u\frac{\partial\theta}{\partial x}+\varepsilon v\frac{\partial\theta}{\partial y}=-\varepsilon\sin{\theta}\cos{\theta}\frac{\partial u}{\partial x}-\sin^{2}{\theta}\frac{\partial u}{\partial y}+\varepsilon^{2}\cos^{2}{\theta}\frac{\partial v}{\partial x}+\varepsilon\sin{\theta}\cos{\theta}\frac{\partial v}{\partial y}.

(A.4)

Appendix B Simplification of the equation for $\theta$

In the main text, we claimed that equation (3.8) permitted great simplification by noting that the equation corresponded only to advection in a purely horizontal direction the reference domain. To demonstrate this simplification, suppose $\tilde{\theta}\left(x^{\prime},y^{\prime},t\right)$ is a function defined over $D_{\text{ref}}$ that satisfies the advection equation

\frac{\partial\tilde{\theta}}{\partial t}+\tilde{u}\frac{\partial\tilde{\theta}}{\partial x^{\prime}}=\tilde{f},

(B.1)

where $\tilde{u}\left(x^{\prime},t\right)$ and $\tilde{f}\left(x^{\prime},t\right)$ are a horizontal advection velocity and forcing term respectively. We can relate $\theta(x,y,t)=\tilde{\theta}(x(x^{\prime},t),y(x^{\prime},y^{\prime},t),t)$ , and using the mapping $\boldsymbol{\Phi}$ , which gives

	$\displaystyle\frac{\partial\tilde{\theta}}{\partial t}=\frac{\partial\theta}{\partial t}+\frac{\dot{L}x}{L}\frac{\partial\theta}{\partial x}+\left(\frac{\partial H}{\partial t}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial t}+\frac{\dot{L}x}{L}\left(\frac{\partial H}{\partial x}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial x}\right)\right)\frac{\partial\theta}{\partial y},$		(B.2)
	$\displaystyle\frac{\partial\tilde{\theta}}{\partial x^{\prime}}=L\frac{\partial\theta}{\partial x}+L\left(\frac{\partial H}{\partial x}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial x}\right)\frac{\partial\theta}{\partial y}.$		(B.3)

Substituting (B.2)-(B.3) into (B.1), then yields

\frac{\partial\theta}{\partial t}+\frac{\partial\theta}{\partial x}\left(L\tilde{u}+\frac{\dot{L}x}{L}\right)+\frac{\partial\theta}{\partial y}\left(\frac{\partial H}{\partial t}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial t}+\left(\frac{\partial H}{\partial x}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial x}\right)\left(L\tilde{u}+\frac{\dot{L}x}{L}\right)\right)=\tilde{f},

(B.4)

we now choose $u\left(x,t\right)=L\tilde{u}+\dfrac{\dot{L}x}{L}$ , in order to recover the correct coefficient of $\dfrac{\partial\theta}{\partial x}$ . Examining the coefficient of the $\dfrac{\partial\theta}{\partial y}$ term we note that

	$\displaystyle\frac{\partial H}{\partial t}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial t}+u\left(\frac{\partial H}{\partial x}+\left(\frac{y-H}{h}\right)\frac{\partial h}{\partial x}\right)$		(B.5)
	$\displaystyle=\frac{\partial H}{\partial t}+u\frac{\partial H}{\partial x}+\left(\frac{y-H}{h}\right)\left(-h\frac{\partial u}{\partial x}\right)$		(B.6)
	$\displaystyle=\frac{\partial H}{\partial t}+\frac{\partial}{\partial x}\left(uH\right)-y\frac{\partial u}{\partial x}=v$		(B.7)

where we have used the equation for conservation of mass, (2.7), to obtain (B.6). Here, we have demonstrated that the coefficient of $\theta_{y}$ is precisely $v$ when mapping back from $D_{ref}$ to the original domain. Therefore, we have shown that the advection of $\theta$ is purely horizontal upon the reference domain, with velocity $\tilde{u}\left(x^{\prime},t\right)=\dfrac{u-\dot{L}x^{\prime}}{L}=\dfrac{u-u_{\text{mesh}}}{L}$ .

Appendix C Discretisation of the Green and Friedman integral equations

In what follows, the treatment of the integral equations is in the Eulerian framework. The integral equations (3.11),(3.12) require further treatment before being discretised and solved. Introduce

\displaystyle F\left(x,y,t\right)=\frac{\mu_{1}\cos 2\theta+\left(4+4\mu_{3}+\mu_{2}\right)u_{x}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta},

(C.1)

it will be convenient to write $F=F_{1}+u_{x}F_{2}$ , where

\displaystyle F_{1}=\frac{\mu_{1}\cos 2\theta}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta},

\displaystyle F_{2}=\frac{4+4\mu_{3}+\mu_{2}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta},

(C.2)

in much the same way, we also introduce notation for the integrals of $F$ , by defining $G=G_{1}+u_{x}G_{2},$ where

G_{m}\left(x,y,t\right)=\int\limits_{H^{-}}^{y}F_{m}\left(x,s,t\right)\mathrm{d}s,

(C.3)

for $m=1,2$ . Equation (3.11) can now be written as

\displaystyle 0=\frac{\partial}{\partial x}\left(G_{1}\left(x,H^{+},t\right)+u_{x}\left(x,t\right)G_{2}\left(x,H^{+},t\right)\right),

(C.4)

using the trapezoidal rule,

G_{m}\left(x,H^{+},t\right)=\frac{h\left(x,t\right)}{N-1}\left(\frac{F_{m}\left(x,y_{0},t\right)+F_{m}\left(x,y_{N-1},t\right)}{2}+\sum_{i=1}^{N-2}F_{m}\left(x,y_{i},t\right)\right),

(C.5)

where $N$ is the number of nodes in the $y$ -direction. We note that upon substituting the trapezoidal rule into (C.4), the result does not depend on $H$ . That is, its presence in the integration limits is redundant and effectively just describes a vertical translation. Therefore, $H$ is decoupled from the rest of the system and we can solve the equations for $u$ and $\theta$ and then compute $H$ as required at a desired time. For completeness, we include the discretisation for equation (C.4). First, introduce the notation

\displaystyle\left[G_{m}\right]_{i,N-1}^{k}=G_{m}\left(x_{i},H_{i}^{k}+h_{i}^{k}/2,t_{k}\right),

(C.6)

now (C.4) gives us, through centred finite differences,

\frac{\left[G_{1}\right]^{k}_{i-1,N-1}-\left[G_{1}\right]^{k}_{i+1,N-1}}{2L/\left(M-1\right)}=\frac{}{}\frac{\left[G_{2}\right]^{k}_{i+1,N-1}-\left[G_{2}\right]^{k}_{i-1,N-1}}{2L/\left(M-1\right)}\frac{u_{i+1}^{k}-u_{i}^{k}}{2L/\left(M-1\right)}\\ +\frac{u_{i}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{\left(L/\left(M-1\right)\right)^{2}}\left[G_{2}\right]_{i,N-1}^{k},

(C.7)

where $N,M$ are the number of nodes in the vertical and horizontal directions respectively so that $i=0:M-1,j=0:N-1$ . Noting that $u_{0}^{k}=0$ , $u_{M-1}^{k}=\dot{L}\left(t_{k}\right)$ , and that $G_{m}$ is readily precomputed at each time-step $k$ , yields a tri-diagonal system for $u^{k}$ . If we now consider the equation for $H$ , (3.9), this is the only equation in the model that is indeed easier to treat in the Eulerian framework than the ALE. Using equation (2.12) and the Leibniz rule we may write equation (2.9) as

0=-\left(H_{xx}+\frac{h_{xx}}{2}\right)G\left(x,H^{+},t\right)+\frac{\partial^{2}}{\partial x^{2}}\int\limits_{H^{-}}^{H^{+}}G(x,y,t)\mathrm{d}y,

(C.8)

in order to proceed, one must apply the trapezoidal rule twice to each $G_{m}$ . Applying it once yields

\int\limits_{H^{-}}^{H^{+}}G_{m}\left(x,y,t\right)\mathrm{d}y=\frac{h\left(x,t\right)}{N-1}\left(\frac{G_{m}\left(x,y_{0},t\right)+G_{m}\left(x,y_{N-1},t\right)}{2}+\sum_{j=1}^{N-2}G_{m}\left(x,y_{j},t\right)\right),

(C.9)

then, for each $j>0$ ,

G_{m}\left(x,y_{j},t\right)=\frac{h\left(x,t\right)}{N-1}\left(\frac{F_{m}\left(x,y_{0},t\right)+F_{m}\left(x,y_{j},t\right)}{2}+\sum_{i=1}^{j-1}F_{m}\left(x,y_{i},t\right)\right),

(C.10)

for the case $j=0$ , $G_{m}\left(x,y_{0},t\right)=0$ . Substitution of (C.10) into (C.9) yields

\int\limits_{H^{-}}^{H^{+}}G_{m}\left(x,y,t\right)\mathrm{d}y=\left(\frac{h}{N-1}\right)^{2}\left(\frac{2N-3}{4}F_{m}\left(x,y_{0},t\right)+\frac{1}{4}F_{m}\left(x,y_{N-1},t\right)\right.\\ \left.+\sum_{j=1}^{N-2}\left(N-1-j\right)F_{m}\left(x,y_{j},t\right)\right).

(C.11)

Finally, we require the introduction of the notation

\displaystyle\left[IG_{m}\right]^{k}_{i}=\int\limits_{H^{-}}^{H^{+}}G_{m}\left(x_{i},y,t_{k}\right)\mathrm{d}y,

(C.12)

for $m=1,2$ . Clearly, we use (C.11) to precompute $[IG_{m}]$ at the required nodes as necessary. The discretisation of equation (C.8) is then

\frac{H_{i-1}^{k}-2H_{i}^{k}+H_{i+1}^{k}}{\left(L(t_{k})/\left(M-1\right)\right)^{2}}\left(\left[G_{1}\right]_{i}^{k}+\frac{u_{i+1}^{k}-u_{i-1}^{k}}{2L\left(t_{k}\right)/\left(M-1\right)}\left[G_{2}\right]_{i}^{k}\right)\\ =-\frac{h_{i-1}^{k}-2h_{i}^{k}+h_{i+1}^{k}}{2\left(L(t_{k})/\left(M-1\right)\right)^{2}}\left(\left[G_{1}\right]_{i}^{k}+\frac{u_{i+1}^{k}-u_{i-1}^{k}}{2L\left(t_{k}\right)/\left(M-1\right)}\left[G_{2}\right]_{i}^{k}\right)+\frac{\left[IG_{1}\right]_{i-1}^{k}-2\left[IG_{1}\right]_{i}^{k}+\left[IG_{1}\right]_{i+1}^{k}}{\left(L\left(t_{k}\right)/\left(M-1\right)\right)^{2}}\\ +\frac{u_{i+2}^{k}-2u_{i+1}^{k}+2u_{i-1}^{k}-u_{i-2}^{k}}{2\left(L\left(t_{k}\right)/\left(M-1\right)\right)^{3}}\left[IG_{2}\right]_{i}^{k+1}+\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{\left(L\left(t_{k}\right)/\left(M-1\right)\right)^{2}}\frac{\left[IG_{2}\right]_{i+1}^{k}-\left[IG_{2}\right]_{i-1}^{k}}{2L\left(t_{k}\right)/\left(M-1\right)}\\ +\frac{\left[IG_{2}\right]_{i+1}^{k}-2\left[IG_{2}\right]_{i}^{k}+\left[IG_{2}\right]_{i-1}^{k}}{\left(L\left(t_{k}\right)/\left(M-1\right)\right)^{2}}\frac{u_{i+1}^{k}-u_{i-1}^{k}}{2L\left(t_{k}\right)/\left(M-1\right)}.

(C.13)

Equation (C.13) contains wholly precomputable quantities on the RHS. Therefore, similar to the equation for $u$ , this creates a tri-diagonal system to be solved for $H$ . For the specific cases of $i=0,M-1$ , we have the boundary condition that $H_{0}^{k}=H_{M-1}^{k}=0$ . For the cases of $i=1,M-2$ , the discretisation must be modified slightly as the stencil for $u_{xxx}$ is too wide. This can be done in a number of ways and is omitted.

Appendix D Derivation of the short timescale integral equations

In order to close the model, we must go to yet higher orders in order to obtain equations for $u^{(0)},H^{(0)}$ . Our approach is to integrate the relevant equations over the depth of the sheet and apply the no-stress boundary conditions at $y=H^{(0)^{\pm}}$ . We find that a number of higher order terms in the boundary conditions will be eliminated by substitution of previously obtained quantities. The $\mathcal{O}\left(\varepsilon^{3}\right)$ $x$ -momentum equation is

\frac{\partial}{\partial x}\left[-p^{(0)}+2\frac{\partial u^{(0)}}{\partial x}+\mu_{1}\cos^{2}\theta^{(0)}\right.+\mu_{2}\left(\cos^{4}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\cos^{3}\theta^{(0)}\sin\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right.\\ \left.+\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}\frac{\partial V^{(2)}}{\partial y}\right)+\left.2\mu_{3}\left(2\cos^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\cos\theta^{(0)}\sin\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right)\right]-\frac{\partial^{2}u^{(0)}}{\partial x^{2}}\\ =-\frac{\partial}{\partial y}\left[\frac{\partial u^{(2)}}{\partial y}+\frac{\partial V^{(2)}}{\partial x}+\mu_{1}\left(\theta^{(1)}\cos 2\theta^{(0)}\right)+\mu_{2}\left(\cos^{3}\theta^{(0)}\sin\theta^{(0)}\frac{\partial u^{(1)}}{\partial x}\right.\right.\\ +\theta^{(1)}\left(\cos^{4}\theta^{(0)}-3\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}\right)\frac{\partial u^{(0)}}{\partial x}+\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}\left(\frac{\partial u^{(2)}}{\partial y}+\frac{\partial V^{(2)}}{\partial x}\right)\\ +\frac{1}{2}\theta^{(1)}\sin 4\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)+\cos\theta^{(0)}\sin^{3}\theta^{(0)}\frac{\partial V^{(3)}}{\partial y}\\ \left.\left.+\theta^{(1)}\left(3\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}-\sin^{4}\theta^{(0)}\right)\frac{\partial V^{(2)}}{\partial y}\right)\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\mu_{3}\left(\frac{\partial u^{(2)}}{\partial y}+\frac{\partial V^{(2)}}{\partial x}\right)\right]\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\frac{\partial^{2}V^{(2)}}{\partial x\partial y},

(D.1)

and its associated no-stress boundary condition is

\frac{\partial u^{(2)}}{\partial y}+\frac{\partial V^{(2)}}{\partial x}+\mu_{1}\theta^{(1)}\cos 2\theta^{(0)}+\mu_{2}\left(\cos^{3}\theta^{(0)}\sin\theta^{(0)}\frac{\partial u^{(1)}}{\partial x}+\theta^{(1)}\left(\cos^{4}\theta^{(0)}-3\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}\right)\frac{\partial u^{(0)}}{\partial x}\right.\\ \left.+\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}\left(\frac{\partial u^{(2)}}{\partial y}+\frac{\partial V^{(2)}}{\partial x}\right)+\frac{1}{2}\theta^{(1)}\sin 4\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)+\cos\theta^{(0)}\sin^{3}\theta^{(0)}\frac{\partial V^{(3)}}{\partial y}\right.\\ \left.+\theta^{(1)}\left(3\sin^{2}\theta^{(0)}\cos^{2}\theta^{(0)}-\sin^{4}\theta^{(0)}\right)\frac{\partial V^{(2)}}{\partial y}\right)+\mu_{3}\left(\frac{\partial V^{(2)}}{\partial x}+\frac{\partial u^{(2)}}{\partial y}\right)\\ =\left(\frac{\partial H^{(0)}}{\partial x}\pm\frac{1}{2}\frac{\partial h^{(0)}}{\partial x}\right)\left[-p^{(0)}+2\frac{\partial u^{(0)}}{\partial x}+\mu_{1}\cos^{2}\theta^{(0)}\right.\\ +\mu_{2}\left(\cos^{4}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\cos^{3}\theta^{(0)}\sin\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)+\cos^{2}\theta^{(0)}\sin^{2}\theta^{(0)}\frac{\partial V^{(2)}}{\partial y}\right)\\ \left.2\mu_{3}\left(2\cos^{2}\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\cos\theta^{(0)}\sin\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right)\right];\text{ on }y=H^{(0)^{\pm}}.

(D.2)

We note that the remaining terms outside of the derivatives in (D.1) cancel due to the continuity equation at $\mathcal{O}(\varepsilon^{2})$ . Integrating equation (D.1) over the depth of the sheet yields

\displaystyle\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\frac{\partial}{\partial x}\left(-p^{(0)}+g_{1}(x,\tau)\right)\mathrm{d}y=-\left[\frac{\partial u^{(2)}}{\partial y}+\frac{\partial V^{(2)}}{\partial x}+g_{2}(x,\tau)\right]_{y=H^{(0)^{-}}}^{y=H^{(0)^{+}}},

(D.3)

where $g_{1},g_{2}$ are functions containing the collected $\mu_{1},\mu_{2},\mu_{3}$ terms from (D.1) and are readily obtained by inspection. Application of (D.2) now yields

\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\frac{\partial}{\partial x}\left(-p^{(0)}+g_{1}\right)\mathrm{d}y=\left(\frac{\partial H^{(0)}}{\partial x}-\frac{1}{2}\frac{\partial h^{(0)}}{\partial x}\right)\left(-p^{(0)}+g_{1}\right)_{y=H^{(0)^{-}}}\\ -\left(\frac{\partial H^{(0)}}{\partial x}+\frac{1}{2}\frac{\partial h^{(0)}}{\partial x}\right)\left(-p^{(0)}+g_{1}\right)_{y=H^{(0)^{+}}},

(D.4)

and hence by use of the Leibniz rule, we now obtain an equation for $\bar{u}$

\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\left(4\frac{\partial u^{(0)}}{\partial x}+\mu_{1}\cos 2\theta^{(0)}+\mu_{2}\left(\cos^{2}2\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\frac{1}{4}\sin 4\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right)\right.\\ \left.+4\mu_{3}\frac{\partial u^{(0)}}{\partial x}\right)\mathrm{d}y=0.

(D.5)

A similar process at $\mathcal{O}\left(\varepsilon^{4}\right)$ yields an equation for $H$ ,

\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{H^{(0)^{+}}}\frac{\partial}{\partial x}\int\limits_{H^{(0)^{-}}}^{y}\left(4\frac{\partial u^{(0)}}{\partial x}+\mu_{1}\cos 2\theta^{(0)}+\mu_{2}\cos^{2}2\theta^{(0)}\frac{\partial u^{(0)}}{\partial x}+\frac{\mu_{2}}{4}\sin 4\theta^{(0)}\left(\frac{\partial u^{(1)}}{\partial y}+\frac{\partial V^{(1)}}{\partial x}\right)\right.\\ \left.+4\mu_{3}\frac{\partial u^{(0)}}{\partial x}\right)\mathrm{d}y^{\prime}\mathrm{d}y=0.

(D.7)

We note that these equations are of the same form as equations (2.8) and (2.9), with the difference being that $u^{(0)}$ now possesses $y$ -dependence. Substitution of (2.12) and (2.27) into both (D.5) and (D.7) and use of the Liebniz rule upon (D.7) yield equations (2.28) and (2.29).

Appendix E Discretisation of the short timescale model

In this section, we give the discretisation of the integral equations in the system of the short timescale equations, namely (3.19), (3.20). As before, we drop the superscript notation for leading-order quantities. We discretise with

\theta(x_{i},y_{j},\tau_{k})=\theta_{i,j}^{k},\text{ etc}

(E.1)

and for functions that contain subscripts:

Z_{1}(x_{i},y_{j},\tau_{k})=[Z_{1}]_{i,j}^{k},\text{ etc}

(E.2)

where $i=1:M-1,j=1:N-1,$ where $M,N$ are the number of nodes in the $x,y$ -directions.Similarly to the discretisation of the Green and Friedman model equations, we define

	$\displaystyle Z_{1}(x,y,\tau)=\int\limits_{-\frac{1}{2}}^{y}\frac{\mu_{1}\cos 2\theta}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}\tilde{y},$		(E.3)
	$\displaystyle Z_{2}(x,y,\tau)=\int\limits_{-\frac{1}{2}}^{y}\frac{4+4\mu_{3}+\mu_{2}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h\mathrm{d}\tilde{y},$		(E.4)
	$\displaystyle J(x,t)=\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\frac{\left(4+4\mu_{3}+\mu_{2}\right)h^{2}y^{\prime}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}\mathrm{d}\tilde{y},$		(E.5)

we note that the definitions of $Z_{1},Z_{2}$ here are similar to $G_{1},G_{2}$ in appendix C, but are not precisely the same. It is possible to obtain $G_{1},G_{2}$ from $Z_{1},Z_{2}$ by undoing both the transformation $y=H+h\tilde{y}$ and the short-timescale. Since the ALE transformation in Section 3 acts as a linear transformation on the integral equations, one may treat $Z_{1},Z_{2}$ as the integrals in ALE form on the short-timescale. We may now rewrite equation (3.19) as

Z_{2}(x,\frac{1}{2},\tau)\frac{\partial^{2}\bar{u}}{\partial x^{2}}+\frac{\partial Z_{2}(x,\frac{1}{2},\tau)}{\partial x}\frac{\partial\bar{u}}{\partial x}+\frac{\partial^{2}H}{\partial x\partial\tau}\left(\frac{\partial Z_{2}(x,\frac{1}{2},\tau)}{\partial x}\frac{\partial H}{\partial x}+Z_{2}(x,\frac{1}{2},\tau)\frac{\partial^{2}H}{\partial x^{2}}\right)\\ +\frac{\partial^{3}H}{\partial x^{2}\partial\tau}\left(Z_{2}(x,\frac{1}{2},\tau)\frac{\partial H}{\partial x}-\frac{\partial J}{\partial x}\right)-J\frac{\partial^{4}H}{\partial x^{3}\partial\tau}=-\frac{\partial Z_{1}(x,\frac{1}{2},\tau)}{\partial x}.

(E.6)

We choose to use a FTCS finite difference method, hence the discretisation of (E.6) is

[Z_{2}]_{i,N-1}^{k}\frac{\bar{u}_{i+1}^{k}-2\bar{u}_{i}^{k}+\bar{u}_{i-1}^{k}}{\Delta x^{2}}+\frac{[Z_{2}]_{i+1,N-1}^{k}-[Z_{2}]_{i-1,N-1}^{k}}{2\Delta x}\frac{\bar{u}_{i+1}^{k}-\bar{u}_{i-1}^{k}}{2\Delta x}+\\ \frac{\left(H_{i+1}^{k+1}-H_{i-1}^{k+1}\right)-\left(H_{i+1}^{k}-H_{i-1}^{k}\right)}{2\Delta\tau\Delta x}\left(\frac{[Z_{2}]_{i+1,N-1}^{k}-[Z_{2}]_{i-1,N-1}^{k}}{2\Delta x}\frac{H_{i+1}^{k}-H_{i-1}^{k}}{2\Delta x}\right.\\ \left.+[Z_{2}]_{i,N-1}^{k}\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}\right)\\ +\frac{\left(H_{i+1}^{k+1}-2H_{i}^{k+1}+H_{i-1}^{k+1}\right)-\left(H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}\right)}{\Delta\tau\Delta x^{2}}\left([Z_{2}]_{i,N-1}^{k}\frac{H_{i+1}^{k}-H_{i-1}^{k}}{2\Delta x}\right.\\ \left.-\frac{J_{i+1}^{k}-J_{i-1}^{k}}{2\Delta x}\right)-J_{i}^{k}\frac{\left(H_{i+2}^{k+1}-2H_{i+1}^{k+1}+2H_{i-1}^{k+1}-H_{i-2}^{k+1}\right)-\left(H_{i+2}^{k}-2H_{i+1}^{k}+2H_{i-1}^{k}-H_{i-2}^{k}\right)}{2\Delta\tau\Delta x^{3}}\\ =-\frac{[Z_{1}]_{i+1,N-1}^{k}-[Z_{1}]_{i-1,N-1}^{k}}{2\Delta x},

(E.7)

This discretisation can be written in the matrix form

\Delta\tau\Delta x^{2}\mathbf{b}+\mathbf{M_{H}}\mathbf{H}^{k}=\mathbf{M_{H}}\mathbf{H}^{k+1}+\Delta\tau\Delta x\mathbf{M_{\bar{U}}}\mathbf{\bar{u}}^{k},

(E.8)

where $\boldsymbol{H}^{k}=\left(H_{1}^{k},H_{2}^{k},\dots,H_{N+1}^{k}\right)^{T}$ , $\boldsymbol{\bar{u}}^{k}=\left(\bar{u}_{1}^{k},\bar{u}_{2}^{k},\dots,\bar{u}_{N+1}^{k}\right)^{T}$ , and $\boldsymbol{M_{H}},\boldsymbol{M_{\bar{U}}}$ are matrices whose entries are the coefficients of the $H^{k+1}$ and $\bar{u}^{k}$ terms respectively, and are dependent upon the choice of discretisation of the $x$ -derivatives of $H,\bar{u}$ . We note that the left hand side of (E.8) is known and precomputable at each time step. Due to the $\dfrac{\partial H^{4}}{\partial x^{3}\partial\tau}$ term, the stencil must be adjusted at the ends of the domain, by using biased finite differences.
Next, we define more functions for notational convenience:

	$\displaystyle IZ_{1}(x,\tau)=\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{\tilde{y}}\frac{\mu_{1}\cos 2\theta}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}\tilde{y}^{\prime}\mathrm{d}\tilde{y},$		(E.9)
	$\displaystyle IZ_{2}(x,\tau)=\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{\tilde{y}}\frac{4+4\mu_{3}+\mu_{2}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}h^{2}\mathrm{d}\tilde{y}^{\prime}\mathrm{d}\tilde{y},$		(E.10)
	$\displaystyle K(x,\tau)=\int\limits_{-\frac{1}{2}}^{\frac{1}{2}}\int\limits_{-\frac{1}{2}}^{\tilde{y}}\frac{\left(4+4\mu_{3}+\mu_{2}\right)h^{3}\tilde{y}^{\prime}}{4+4\mu_{3}+\mu_{2}\sin^{2}2\theta}\mathrm{d}\tilde{y}^{\prime}\mathrm{d}\tilde{y}.$		(E.11)

With the introduced functions, we may express (3.20) as

IZ_{2}\frac{\partial^{3}\bar{u}}{\partial x^{3}}+2\frac{\partial IZ_{2}}{\partial x}\frac{\partial^{2}\bar{u}}{\partial x^{2}}+\frac{\partial\bar{u}}{\partial x}\left(\frac{\partial^{2}IZ_{2}}{\partial x^{2}}-Z_{2}\left(\frac{\partial^{2}H}{\partial x^{2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{2}}\right)\right)\\ -K\frac{\partial^{5}H}{\partial x^{4}\partial\tau}+\frac{\partial^{4}H}{\partial x^{3}\partial\tau}\left(IZ_{2}\frac{\partial H}{\partial x}-2\frac{\partial K}{\partial x}\right)\\ +\frac{\partial^{3}H}{\partial x^{2}\partial\tau}\left(2IZ_{2}\frac{\partial^{2}H}{\partial x^{2}}+2\frac{\partial IZ_{2}}{\partial x}\frac{\partial H}{\partial x}+J\left(\frac{\partial^{2}H}{\partial x^{2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{2}}\right)-\frac{\partial^{2}K}{\partial x^{2}}\right)\\ +\frac{\partial^{2}H}{\partial x\partial\tau}\left(IZ_{2}\frac{\partial^{3}H}{\partial x^{3}}+2\frac{\partial IZ_{2}}{\partial x}\frac{\partial^{2}H}{\partial x^{2}}+\frac{\partial^{2}IZ_{2}}{\partial x^{2}}\frac{\partial H}{\partial x}-Z_{2}\frac{\partial H}{\partial x}\left(\frac{\partial^{2}H}{\partial x^{2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{2}}\right)\right)\\ =Z_{1}\left(\frac{\partial^{2}H}{\partial x^{2}}+\frac{1}{2}\frac{\partial^{2}h}{\partial x^{2}}\right)-\frac{\partial^{2}IZ_{1}}{\partial x^{2}}.

(E.12)

The discretisation of (E.12) is:

[IZ_{2}]_{i}^{k}\frac{\bar{u}_{i+2}^{k}-2\bar{u}_{i+1}^{k}+2\bar{u}_{i-1}^{k}-\bar{u}_{i-2}^{k}}{2\Delta x^{3}}+\frac{[IZ_{2}]_{i+1}^{k}-[IZ_{2}]_{i-1}^{k}}{\Delta x}\frac{\bar{u}_{i+1}^{k}-2\bar{u}_{i}^{k}+\bar{u}_{i-1}^{k}}{\Delta x^{2}}\\ +\frac{\bar{u}_{i+1}^{k}-\bar{u}_{i-1}^{k}}{2\Delta x}\left(\frac{[IZ_{2}]_{i+1}^{k}-2[IZ_{2}]_{i}^{k}+[IZ_{2}]_{i-1}^{k}}{\Delta x^{2}}-Z_{2_{i}}\left(\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}+\frac{1}{2}\frac{h_{i+1}^{k}-2h_{i}^{k}+h_{i-1}^{k}}{\Delta x^{2}}\right)\right)\\ +\frac{\left(H_{i+1}^{k+1}-H_{i-1}^{k+1}\right)-\left(H_{i+1}^{k}-H_{i-1}^{k}\right)}{2\Delta\tau\Delta x}\left(\frac{[IZ_{2}]_{i+1}^{k}-2[IZ_{2}]_{i}^{k}+[IZ_{2}]_{i-1}^{k}}{\Delta x^{2}}\frac{H_{i+1}^{k}-H_{i-1}^{k}}{2\Delta x}\right.\\ \left.+\frac{[IZ_{2}]_{i+1}^{k}-[IZ_{2}]_{i-1}^{k}}{\Delta x}\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}\right.\\ \left.+[IZ_{2}]_{i}^{k}\frac{H_{i+2}^{k}-2H_{i+1}^{k}+2H_{i-1}^{k}-H_{i-2}^{k}}{2\Delta x^{3}}-Z_{2_{i}}^{k}\frac{H_{i+1}^{k}-H_{i-1}^{k}}{2\Delta x}\left(\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}+\frac{1}{2}\frac{h_{i+1}^{k}-2h_{i}^{k}+h_{i-1}^{k}}{\Delta x^{2}}\right)\right)\\ +\frac{\left(H_{i+1}^{k+1}-2H_{i}^{k+1}H_{i-1}^{k+1}\right)-\left(H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}\right)}{\Delta\tau\Delta x^{2}}\left(2[IZ_{2}]_{i}^{k}\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}+\frac{IZ_{2_{i+1}}-IZ_{2_{i-1}}}{\Delta x}\frac{H_{i+1}^{k}-H_{i-1}^{k}}{2\Delta x}\right.\\ \left.-\frac{K_{i+1}^{k}-2K_{i}^{k}+K_{i-1}^{k}}{\Delta x^{2}}+J_{i}^{k}\left(\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}+\frac{1}{2}\frac{h_{i+1}^{k}-2h_{i}^{k}+h_{i-1}^{k}}{\Delta x^{2}}\right)\right)\\ +\frac{\left(H_{i+2}^{k+1}-2H_{i+1}^{k+1}+2H_{i-1}^{k+1}-H_{i-2}^{k+1}\right)-\left(H_{i+2}^{k}-2H_{i+1}^{k}+2H_{i-1}^{k}-H_{i-2}^{k}\right)}{2\Delta\tau\Delta x^{3}}\left(IZ_{2_{i}}^{k}\frac{H_{i+1}^{k}-H_{i-1}^{k}}{2\Delta x}-\frac{K_{i+1}^{k}-K_{i-1}^{k}}{2\Delta x}\right)\\ -K_{i}^{k}\frac{\left(H_{i+2}^{k+1}-4H_{i+1}^{k+1}+6H_{i}^{k+1}-4H_{i-1}^{k+1}+H_{i-2}^{k+1}\right)-\left(H_{i+2}^{k}-4H_{i+1}^{k}+6H_{i}^{k}-4H_{i-1}^{k}+H_{i-2}^{k}\right)}{\Delta\tau\Delta x^{4}}\\ =Z_{1_{i}}^{k}\left(\frac{H_{i+1}^{k}-2H_{i}^{k}+H_{i-1}^{k}}{\Delta x^{2}}+\frac{1}{2}\frac{h_{i+1}^{k}-2h_{i}^{k}+h_{i-1}^{k}}{\Delta x^{2}}\right)-\frac{IZ_{2_{i+1}}^{k}-2IZ_{2_{i}}^{k}+IZ_{2_{i-1}}^{k}}{\Delta x^{2}}.

(E.13)

As before, we may write (E.13) in the form

\Delta\tau\Delta x^{2}\mathbf{c}+\mathbf{M_{H}^{\prime}}\mathbf{H}^{k}=\mathbf{M_{H}^{\prime}}\mathbf{H}^{k+1}+\Delta\tau\Delta x\mathbf{M_{\bar{U}}^{\prime}}\mathbf{\bar{u}}^{k},

(E.14)

where the coefficients of all $H^{k+1},\bar{u}^{k}$ terms are entries within $\boldsymbol{M^{\prime}_{H}},\boldsymbol{M^{\prime}_{\bar{U}}}$ respectively.
In our implementation, the matrices $\boldsymbol{M^{\prime}_{\bar{U}}},\boldsymbol{M^{\prime}_{H}}$ are a quintuple banded matrices. The additional derivative in $H$ does not change the size of the stencil. We include the boundary conditions for $H$ in the first two and final two lines of both matrices, so that it is not necessary to adjust the stencil near the endpoints of the domain. Hence, we can construct the linear system

\displaystyle\begin{pmatrix}\mathbf{M_{H}}&\Delta\tau\Delta x\mathbf{M_{\bar{U}}}\\ \mathbf{M_{H}^{\prime}}&\Delta\tau\Delta x\mathbf{M_{\bar{U}}^{\prime}}\end{pmatrix}\begin{pmatrix}\mathbf{H}^{k+1}\\ \mathbf{\bar{u}}^{k}\end{pmatrix}=\begin{pmatrix}\Delta\tau\Delta x^{2}\mathbf{b}+\mathbf{M_{H}}\mathbf{H}^{k}\\ \Delta\tau\Delta x^{2}\mathbf{c}+\mathbf{M_{H}^{\prime}}\mathbf{H}^{k}\end{pmatrix}.

(E.15)

This linear system is solved at each $k$ , with equation $\eqref{FirstTheta}$ requring $H^{k+1}$ in order to update $\theta^{k}\rightarrow\theta^{k+1}$ .

References

[1] Merav Antman-Passig, Shahar Levy, Chaim Gartenberg, Hadas Schori, and Orit Shefi. Mechanically oriented 3d collagen hydrogel for directing neurite growth. Tissue Engineering Part A, 23(9-10):403–414, 2017.
[2] JD Buckmaster, A Nachman, and L Ting. The buckling and stretching of a viscida. J Fluid Mech, 69(1):1–20, 1975.
[3] Jeevanjyoti Chakraborty, Jingxi Luo, and Rosemary J Dyson. Lockhart with a twist: Modelling cellulose microfibril deposition and reorientation reveals twisting plant cell growth mechanisms. Journal of Theoretical Biology, 525:110736, 2021.
[4] FC Chretien and G David. Temporary obstructive effect of human cervical mucus on spermatozoa throughout reproductive life: a scanning electron microscopic study. Eur J Obstet Gynecol Reprod Biol, 8(6):307–321, 1978.
[5] LJ Cummings. Evolution of a thin film of nematic liquid crystal with anisotropic surface energy. Eur J Appl Math, 15(6):651–677, 2004.
[6] G Cupples, RJ Dyson, and DJ Smith. Viscous propulsion in active transversely isotropic media. J Fluid Mech, 812:501–524, 2017.
[7] J Donea, A Huerta, JP Ponthot, and A Rodríguez-Ferran. Arbitrary Lagrangian–Eulerian methods. Encyclopedia of Computational Mechanics Second Edition, pages 1–23, 2017.
[8] RJ Dyson, JEF Green, JP Whiteley, and HM Byrne. An investigation of the influence of extracellular matrix anisotropy and cell–matrix interactions on tissue architecture. J. Math. Biol, 72(7):1775–1809, 2016.
[9] RJ Dyson and OE Jensen. A fibre-reinforced fluid model of anisotropic plant cell growth. J Fluid Mech, 655:472–503, 2010.
[10] JL Ericksen. Transversely isotropic fluids. Kolloid Z, 173(2):117, 1960.
[11] E Evans-Hoeker, DA Pritchard, DL Long, AH Herring, JB Stanford, and AZ Steiner. Cervical mucus monitoring prevalence and associated fecundability in women trying to conceive. Fertil Steril, 100(4):1033–1038, 2013.
[12] JEF Green and A Friedman. The extensional flow of a thin sheet of incompressible, transversely isotropic fluid. Eur J Appl Math, 19(3):225–257, 2008.
[13] CR Holloway, G Cupples, DJ Smith, JEF Green, RJ Clarke, and RJ Dyson. Influences of transversely isotropic rheology and translational diffusion on the stability of active suspensions. Roy Soc Open Sci, 5(8):180456, 2018.
[14] CR Holloway, RJ Dyson, and DJ Smith. Linear Taylor–Couette stability of a transversely isotropic fluid. Proc. R. Soc. A, 471(2178):20150141, 2015.
[15] CR Holloway, DJ Smith, and RJ Dyson. Linear Rayleigh–Bénard stability of a transversely isotropic fluid. Eur J Appl Math, pages 1–23, 2018.
[16] PD Howell. Extensional thin layer flows. PhD thesis, University of Oxford, 1994.
[17] PD Howell. Models for thin viscous sheets. Eur J Appl Math, 7(4):321–343, 1996.
[18] BD Hull, TG Rogers, and AJM Spencer. Theory of fibre buckling and wrinkling in shear flows of fibre-reinforced composites. Compos Manuf, 2(3-4):185–191, 1991.
[19] MEM Lee and H Ockendon. A continuum model for entangled fibres. Eur J Appl Math, 16(2):145–160, 2005.
[20] KS Moghissi. The function of the cervix in fertility. Fertil Steril, 23(4):295, 1972.
[21] G Pfingstag, B Audoly, and A Boudaoud. Thin viscous sheets with inhomogeneous viscosity. Phys Fluids, 23(6):063103, 2011.
[22] N Phan-Thien and AL Graham. The squeezing flow of a model suspension fluid. Rheologica acta, 29(5):433–441, 1990.
[23] N Phan-Thien and AL Graham. A new constitutive model for fibre suspensions: flow past a sphere. Rheologica acta, 30(1):44–57, 1991.
[24] NM Ribe. Bending and stretching of thin viscous sheets. J Fluid Mech, 433:135–160, 2001.
[25] NM Ribe. A general theory for the dynamics of thin viscous sheets. J Fluid Mech, 457:255–283, 2002.
[26] TG Rogers. Squeezing flow of fibre-reinforced viscous fluids. J Eng Math, 23(1):81–89, 1989.
[27] AJM Spencer. Deformations of fibre-reinforced materials. Oxford Univ Press New York, 1972.
[28] AJM Spencer. Fibre-streamline flows of fibre-reinforced viscous fluids. Eur J Appl Math, 8(2):209–215, 1997.
[29] DP Wolf, L Blasco, MA Khan, and M Litt. Human cervical mucus. IV. Viscoelasticity and sperm penetrability during the ovulatory menstrual cycle. Fertil Steril, 30(2):163–169, 1978.
[30] JJ Wylie, BH Bradshaw-Hajek, and YM Stokes. The evolution of a viscous thread pulled with a prescribed speed. J Fluid Mech, 795:380–408, 2016.

	$\displaystyle u^{}\left(0,y^{},t^{*}\right)=0,$	$\displaystyle u^{}\left(L^{},y^{},t^{}\right)=\dot{L}^{*},$
	$\displaystyle H^{}\left(0,t^{}\right)=0,$	$\displaystyle H^{}\left(L^{},t^{*}\right)=0.$

	$\displaystyle\left(x^{},y^{}\right)=\left(xL_{0},\varepsilon yL_{0}\right),\quad\left(u^{},v^{}\right)$	$\displaystyle=\left(uU,\varepsilon vU\right),\quad p^{}=\dfrac{\mu^{}U}{L_{0}}p,\quad t^{*}=\frac{L_{0}}{U}t,$
	$\displaystyle\left(H,L,h\right)=$	$\displaystyle\left(\varepsilon L_{0}H^{},L_{0}L^{},\varepsilon L_{0}h^{*}\right).$