# 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

###### 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.

$\dot{L}^{*}$$x^{*}$$y^{*}$$y^{*}=H^{*}+\frac{h^{*}}{2}$$y^{*}=H^{*}-\frac{h^{*}}{2}$$y^{*}=H^{*}(x^{*},t^{*})$$L^{*}(t^{*})$$\theta$$u^{*}$$v^{*}$

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}$ at
$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}$
at $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}$ at $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.