# Peeling of linearly elastic sheets using complex fluids at low Reynolds
numbers

Anirudh Venkatesh  Vishal Anand  Vivek Narsimhan<EMAIL_ADDRESS>[
Davidson School of Chemical Engineering, Purdue University, West Lafayette,
Indiana 47907, USA

###### Abstract

We investigate the transient, fluid structure interaction (FSI) of a non-
Newtonian fluid peeling two Hookean sheets at low Reynolds numbers ($Re\ll
1$). Two different non-Newtonian fluids are considered – a simplified Phan-
Thien-Tanner (sPTT) model, and an inelastic fluid with shear thinning
viscosity (generalized Newtonian fluid). In the limit of small gap between the
sheets, we invoke a lubrication approximation and numerically solve for the
gap height between the two sheets during the start-up of a pressure-controlled
flow. What we observe is that for an impulse pressure applied to the sheet
inlet, the peeling front moves diffusively ($x_{f}\sim t^{1/2}$) toward the
end of the sheet when the fluid is Newtonian. However, when one examines a
complex fluid with shear thinning, the propagation front moves sub-diffusively
in time ($x_{f}\sim t^{\alpha},\alpha<0.5$), but ultimately reaches the end
faster due to an order of magnitude larger prefactor for the propagation
speed. We provide scaling analyses and similarity solutions to delineate
several regimes of peeling based on the sheet elasticity, flow Weissenberg
number (for sPTT fluid), and shear thinning exponent (for generalized
Newtonian fluid). To conclude, this study aims to afford to the
experimentalist a system of knowledge to apriori delineate the peeling
characteristics of a certain class of complex fluids.

###### keywords:

fluid–structure interactions , viscoelastic fluid , power-law fluid , plate
theory , transient deformations , microfluidics

††journal: Journal of Non-Newtonian Fluid Mechanics

url]https://viveknarsimhan.wixsite.com/website

##  1 Introduction

Microfluidics- the scientific study of flows at the microscale and their
concomitant phenomena [1] \- has managed to hold and sustain the interests of
researchers [2], engineers [3, 4] , and entrepreneurs [5] for some decades.
Even though the governing principles of fluid mechanics at the macroscale are
applicable unchanged at the microscale, the high surface area to volume ratio
of microscale devices [6] perpetuates a cascade of interesting, even
unexpected, physical phenomena including, slip at the interface [7], flows
actuated by the ‘weak force’ of surface tension [8], electrohydrodynamics [9],
inertialess compressibility [10, 11] among others. For instance, and relevant
to the context of this paper, the softness of the material (e.g., PDMS) used
in microfluidic conduits affects the fluid flow considerably [12] and has been
leveraged in applications such as actuators[13], valves[14], bio-mimetic micro
devices, [15] and soft robotics[16]. Inside these conduits, a slowly-moving
viscous fluid exerts pressure on the soft walls constituted of PDMS. This
deformation alters the domain of fluid flow thereby influencing the flow field
non-trivially. This two-way coupled interaction between the fluid flow and the
structure is analysed under the umbrella term of fluid structure interactions
(FSI). [12, 17].

One of the earliest efforts that brought the idea of a deformable surface in
microfluidics forward was by Gervais et al.[18] where they conducted
experiments to delineate the relationship between flow rate and pressure drop
for a channel with a deformable top wall. They coupled the local hydrodynamic
pressure and the local deformation using an empirically derived fitting
parameter. Being empirically derived, their model could not a priori predict
the deformation and the corresponding relationship between flow rate and
pressure drop for the deformable channel. Later, Kiran Raj et al.[19]
presented a better estimate of the fitting parameter mentioned by Gervais et
al.[18] by using a thick plate theory. However, the empiricism and lack of
completeness in the theory posed a drawback to predict these relations a
priori. Christov et al.[20] presented the first detailed theory that used a
two-way coupled FSI model– they did not assume any fitting relationship
between the pressure and the deformation in the channel wall. This study
helped in explaining all the experimental observations from the past
deformable channel studies [18, 21, 22, 23, 24, 19]. The theory given by
Christov et al. was limited to steady state problems, and it did not deal with
the transient problem. The transient counterpart to the steady state problem
of Christov et al.[20] was solved by Martinez-Calvo et al.[25] when they
studied the start-up flow of a Newtonian fluid in a deformable micro-channel.
They provided a characteristic time scale for the fluid-structure interactions
that depended on the geometry and the material of the channel walls. In all of
these studies, the discussion was limited to the inflation of micro channels
with clamped walls. An in-depth understanding of the dynamics of sheets with
free side walls getting peeled off a surface is much less understood at a
micro scale. In relation to free sheet peeling, Hosoi and Mahadevan [26]
studied the time taken for a sheet to peel off from a surface under the
influence of the flow of a thin layer of glycerine competing against gravity
and the attractive van Der Waals forces, albeit at a macro scale.

Beyond the realm of microfluidics, transient FSI has revealed itself in
several other domains in engineering as well as in nature. In geophysics, the
process of hydraulic fracturing, which is employed to enhance the recovery of
oil and gas in reservoirs, is contingent upon the physics of transient FSI
[27, 28]. Closely related to hydraulic fracturing by way of the underlying
physics, is the process of peeling, wherein the low Re transient FSI is
leveraged to separate layers of soft materials, which may or may not be
previously adhered to each other. The emerging area of soft robotics [10, 29]
provides another fascinating platform for researchers and engineers to employ
transient FSI, to push the boundaries of human endeavor. In biological flows,
on the other hand, transient FSI has been investigated in the realms of both
pulmonary systems (airways) [30] and cardiovascular system (hemodynamics)
[31]. In the former, the structure (tube) is frequently modeled as
collapsible, where the radius /area of the tube decreases as a result of FSI
due to negative transmural pressure [32]. In the latter, however, inflation is
the norm, and the increase in radius of the arteries and capillaries due to
pulsatile flow of blood conveyed by them is a phenomena widely studied and
extensively researched for several decades [33]. More recently, the focus in
biological flows FSI has shifted towards the investigation of brain aneurysms
where the blood vessel inside the brain softens itself and develops a bulge
locally to relieve the high blood pressure [34].

Newtonian rheology is an idealization, and deviations from this idealization
abound both in nature (biological fluids [35]) and in laboratory (high
molecular weight polymers [36]). Such deviations include, but are not limited
to, shear dependent viscosity, nonlinear relationship between stress and
strain rate, memory of the stress field, among others. One of the first forays
into the FSI of non-Newtonian fluids was made by the Anand et al. [37], who
investigated shear thinning fluids inside both thin and thick walled
structures. They deduced that compared to the case of the Newtonian fluids,
the shear thinning fluids show smaller deformation, whilst, counter-
intuitively, the thick-walled channel show higher deformation compared to the
thin-walled case. This study has been further extended to make use of even
more realistic models that incorporate viscoelasticity [38].

The fluid structure interaction in the creeping flow regime has also been
analysed from the point of view of fluid instabilities. Inertialess flow of a
Newtonian fluid in a rigid channel is stable because there are no time
dependent terms; the introduction of softness in the walls alters this status
quo. In general, wall viscosity stabilizes the flow, while wall elasticity
destabilizes it [39, 40]. On the other hand, the non Newtonian rheology of the
fluid further complicates the situation [41]. For the plane Couette flow of a
viscoelastic fluid past a linearly elastic solid, the system becomes
unconditionally stable beyond a critical relaxation time of the fluid: no
amount of solid elasticity can make the flow unstable if the fluid elasticity
is beyond a critical value [42]. For the plane Couette flow of a shear
thinning fluid, shear thinning has a stabilizing effect if the ratio of solid
to fluid thickness is small, while the converse is true for large thickness
ratios [41]. Recently, research into instability of creeping flow FSI has
forayed into the domain of finite (but still low) Reynolds numbers [43].

To summarize, a bird’s eye prospect of the extant literature pertaining to
microscale FSI reveals that we are at the foot of an interesting cusp. On one
hand, a fair amount of ground has been covered in Newtonian FSI research,
where both steady state and transient problems across a wide array of
geometries- thin/thick-walled channels, peeled sheets- have been investigated
with aplomb. On the other hand, research in non-Newtonian FSI has lagged
wherein only a minutiae of steady state problems fluid instabilities have been
analyzed in simple geometries. Therefore in this paper, we attempt to plug
this knowledge gap in non Newtonian FSI literature and answer the following
questions: What are the transient characteristics of peeling actuated by non-
Newtonian fluids ? How do non-Newtonian attributes such as shear thinning and
viscoelasticity affect the deformation, rate of deformation, rate of front
propagation and pressure drop during the peeling of an elastic sheet? What are
the quantitative scales of these dependent variables which may be deduced a a
priori and be made available in the service of the experimentalist and
designer?

The outline of this paper is as follows. Section 2 describes the problem setup
of a non-Newtonian fluid peeling two linearly elastic sheets. This consists of
the geometry, structural mechanics, and the rheological models. We discuss the
properties described by the simplified Phan-Thien Tanner (sPTT) constitutive
relationship, as this model captures many qualitative features observed in
polymeric solutions such as shear thinning, extensional thickening, normal
stress differences, and viscoelasticity. We also discuss the rheological
characteristics of a classic power law model in this section. We conclude
section 2 with a discussion on the physical parameters that could be realized
via experiments. Section 3 discusses the numerical and similarity solutions
for an sPTT fluid. For the case of startup, transient peeling under pressure
controlled flow, we provide scaling relationships and similarity solutions for
the peeling time and the propagation front under the limits of strong
viscoelasticity and weak viscoelasticity, as well as weak and strong peeling
deformations. In section 4, we repeat the same analysis for a generalized
power law fluid, which allows one to examine a much wider range of shear
thinning behavior than the sPTT model. Section 5 gives the final takeaway
messages from this study to facilitate future work in the domain.

(a)

(b)

Figure 1:  Problem schematic. (a) side view and (b) isometric view above the
line of symmetry.

##  2 Problem formulation

###  2.1 Problem geometry

The schematic of the system is similar to that in Elbaz et al.[16] and is
shown in 1. In our discussion, we will denote all quantities with * as
dimensional, while unstarred quantities are dimensionless. We consider an
incompressible, non-Newtonian fluid with density $\rho^{*}$ flowing through a
deformable conduit. The equations that describe the fluid flow are the Cauchy
momentum equations and fluid incompressibility the equation of continuity:

$\rho^{*}\left(\frac{\partial\bm{u^{*}}}{\partial
t^{*}}+\bm{u^{*}}\cdot\nabla^{*}\bm{u}^{*}\right)=-{\nabla^{*}p^{*}}+{\nabla^{*}}\cdot{\tau^{*}}$
(1) $\bm{\nabla^{*}}\cdot\bm{u^{*}}=0$ (2)

In the above equations, $\bm{\tau}^{*}$ is the stress tensor for the fluid.
The non-Newtonian fluids we examine are discussed in Sections 2.3.1 and 2.3.2.

The coordinate system chosen for this problem is such that the $x_{1}^{*}$ and
$x^{*}_{2}$ axes align with the flow and shear gradient directions,
respectively. The origin is at the center of the flow inlet, and thus the
structure is symmetric about the $x_{2}^{*}$ axis. The velocity of the fluid
in the $x_{1}^{*}$ and $x_{2}^{*}$ directions is denoted by
$\bm{u^{*}}=(v_{1}^{*},v_{2}^{*})$, and the characteristic flow speed in the
$x_{1}^{*}$ direction is $U^{*}$. The dimensional gauge pressure is given by
$p^{*}$. The initial undeformed fluid gap between the sheets is $2h_{0}^{*}$
and the length of the channel along the positive $x_{1}^{*}$-axis is $l^{*}$.
We will consider a long, slender conduit such that the inverse channel aspect
ratio $\delta=h_{0}^{*}/l^{*}\ll 1$.

###  2.2 Structural Mechanics

The sheet is treated as a linearly elastic solid with Young’s modulus $E^{*}$,
Poisson’s ratio $\nu$, thickness $t_{0}^{*}$, and width $w^{*}$. Here, we lump
all these parameters into a single elastic stiffness $\kappa^{*}$ that relates
the sheet’s deformation to the local pressure:

$h^{*}=h_{0}^{*}+p^{*}/\kappa^{*}$ (3)

This relationship indicates that the channel walls mimic the behavior of a
Hookean spring. This statement at first may seem crude, but there are several
reasons that justify this expression:

  1. 1.

Thin sheets undergoing peeling predominantly exhibit deformation in transverse
directions, with negligible displacement in flow wise direction [20, 44, 16].
The same behavior is accurately portrayed by the set of linear springs, which
deform only vertically and not horizontally.

  2. 2.

In the limit of the structure being slender, more complicated structural
mechanics models like the classical plate theory, Reissner Mindlin plate
theory, and Donnell shell theory yield a linear relationship between the local
deformation and the hydrodynamic pressure[38, 12]. It has also been shown in a
recent review paper by Christov[12], that several types of peeling geometries,
including, but not limited to, rectangular channels and two dimensional
channels can be modeled as Hookean springs.

Of course, the linear relationship between the local hydrodynamic pressure
load and local deformation has a long and rich history. Starting from the
pioneering work of Coyle [45], which was later “revived” by the twin soft
lubrication papers of Stockthein and Mahadevan [46, 47], the relationship was
also used in for modeling lubrication flows of Newtonian fluids in elastic
cavities [48, 49] In short, the simplicity of the Hookean model coupled with
its capacity to mimic a myriad of structural mechanical models makes it an
obvious choice for this problem. Table 1 summarizes expressions for the
stiffness parameter $\kappa^{*}$ for different types of plates.

An important dimensionless parameter is the FSI parameter $\beta$. It
characterizes the extent of the sheet’s deformation, which is the ratio
between the lubrication pressure and the elastic resistance of the sheet:

$\beta=\frac{P^{*}_{lub}}{h^{*}_{0}\kappa^{*}}$ (4)

The lubrication pressure for the different fluid models will be discussed in
Sections 3.2 and 3.3. The resultant timescale that arises from the sheet’s
fluid-structure interaction is an FSI timescale $t_{FSI}^{*}$:

$t_{FSI}^{*}=\beta\frac{l^{*}}{U^{*}}$ (5)

This quantity is the time it takes for the sheet to deform to its steady
profile under the limit of small deformations.

| Thin sheet $\left(\frac{t_{0}^{*}}{w^{*}}\ll 1\right)$ | Thick sheet $\left(\frac{t_{0}^{*}}{w^{*}}\sim O(1)\right)$
---|---|---
$\kappa^{*}$ | $\dfrac{480E^{*}t_{0}^{*3}}{7w^{*4}(1-\nu^{2})}$ | $1.842455\dfrac{E^{*}}{(1-\nu^{2})w^{*}}$
Table 1: Elastic stiffness $\kappa^{*}$ for various width-averaged structural
mechanics models where $t_{0}^{*}$ is thickness, $w^{*}$ is the width, $E^{*}$
is the Young’s modulus, and $\nu$ is the Poisson’s ratio [50].

###  2.3 Fluid mechanics – constitutive relationship

####  2.3.1 Simple Phan Thien Tanner (sPTT model)

The constitutive model chosen for Section 4 is the simple-Phan Thien Tanner
(sPTT) model. This model utilises the Lodge-Yamamoto type of network theory to
predict the effective slip of a polymeric chain network under the dynamic
action of junction breakage and creation [51]. At constant temperature, the
model captures many qualitative features for the rheology of polymeric
solutions. For instance, under steady shear flow, the model predicts a
positive first normal stress difference, as well shear thinning for both the
viscosity and normal stress differences. Under small amplitude oscillatory
shear, one obtains a response characteristic of a viscoelastic liquid with one
relaxation time.[52]

The following equation describes the constitutive relationship of the sPTT
fluid[52]:

$f\left(\text{tr}\bm{\left(\bm{\tau^{*}}\right)}\right)\bm{\tau^{*}}+\lambda^{*}\bm{\widehat{\tau^{*}}}=2\eta^{*}\bm{D^{*}}$
(6)

In the above equation, $f$ is a function of the trace of the stress tensor and
$\bm{D}^{*}=\frac{1}{2}\bm{(}\bm{\nabla^{*}u^{*}}+\bm{\nabla^{*}u^{*T}}\bm{)}$
is the rate of deformation tensor. The zero shear viscosity of the polymeric
fluid is $\eta^{*}$, the longest relaxation time is $\lambda^{*}$, and
$\bm{\widehat{\tau^{*}}}$ is the Gordon Schowalter derivative:

$\bm{\widehat{\tau^{*}}}=\frac{D\bm{\tau^{*}}}{Dt^{*}}-\left(\nabla\bm{u}^{*}\right)^{T}\cdot\bm{\tau^{*}}-\bm{\tau^{*}}\cdot\left(\nabla\bm{u}^{*}\right)$
(7)

In the above expression, $\frac{D}{Dt^{*}}$ is the material derivative – i.e.,
$\frac{D}{Dt^{*}}=\frac{\partial}{\partial t^{*}}+\bm{u}^{*}\cdot\nabla^{*}$.

The function $f$ can take many forms as long as it satisfies the properties
that $f(0)=1$, $f$ is increasing, and $f(x)\approx x$ for $x\ll 1$. Two of the
most popular forms for $f$ are given by[52]:

$\displaystyle f(\text{tr}(\bm{\tau^{*}}))$
$\displaystyle=\exp\left(\frac{\varepsilon\lambda^{*}}{\eta^{*}}\text{tr}(\bm{\tau^{*}})\right)$
(8a) $\displaystyle f(\text{tr}(\bm{\tau^{*}}))$
$\displaystyle=1+\frac{\varepsilon\lambda^{*}}{\eta^{*}}\text{tr}(\bm{\tau^{*}})$
(8b)

Here, $\varepsilon$ is the elongation parameter that acts as a representative
quantity to understand the elongational stretching and relaxation of the
polymer macro-molecules. The sPTT model given by Eq.(8a) is known as the
exponential sPTT model whereas the linear model in Eq.(8b) gives the
relationship for small (limiting) values of the trace of the stress tensor. We
will consider the linear sPTT model Eq.(8b) in this paper.

One important dimensionless quantity that arises from this model is the
Weissenberg number, which is the product of the polymer relaxation time and
the characteristic shear rate $\dot{\gamma}^{*}=U^{*}/h_{0}^{*}$:

$Wi\equiv\frac{\lambda^{*}U^{*}}{h_{0}^{*}}$ (9)

A value of $Wi\sim O(1)$ indicates that the viscoelastic contributions to the
stress tensor are significant. We will consider potentially large Weissenberg
numbers $Wi\sim O(1)-O(10)$ in this paper but cases where $\delta
Wi\equiv\lambda^{*}U^{*}/l^{*}\ll 1$ – in other words, the polymer relaxation
time is smaller than the average residence time $t^{*}_{conv}=l^{*}/U^{*}$ in
the conduit. In recent times, there have been studies where $\delta\sim
O(10^{-4})$ and lower[53] thereby allowing this simplification. We also note
that for larger values of $Wi$, the onset of viscoelastic turbulence could
affect the estimated predictions [54]. A more detailed stability analysis
(linear or nonlinear) is needed to explain these effects [55].

####  2.3.2 Generalized Newtonian fluid

The constitutive model chosen for Section 4 is a generalized Newtonian fluid,
an inelastic fluid that exhibits shear thinning. In its most general form, the
dimensional form of the stress tensor is:

$\bm{\tau}^{*}=2^{n}K^{*}\left(\sqrt{\bm{D}^{*}:\bm{D}^{*}}\right)^{(n-1)}\bm{D}^{*}$
(10)

where $\tau^{*}$ is the stress tensor, $\bm{D}^{*}$ is rate of strain tensor,
$n$ is a power law index, and $K^{*}$ is a consistency index. In a simple
shear flow $u_{x}^{*}=\dot{\gamma}^{*}y^{*}$, the shear stress behaves as a
Newtonian fluid with an effective viscosity that changes with the shear rate –
i.e., $\tau_{xy}^{*}=\eta_{eff}^{*}\dot{\gamma}^{*}$, where
$\eta_{eff}^{*}=K^{*}|\dot{\gamma}|^{n-1}$. Note – this model works well for
describing the steady shear behavior of many fluids where one observes a power
law scaling for viscosity versus shear rate. Some examples include Xanthan gum
and aqueous polyacrylamide solutions [56, 57]. We also note that this model
describes a wider range of shear thinning behavior than the sPTT model
discussed in section 2.3.1, at the expense of neglecting other effects such as
elasticity and normal stress differences.

###  2.4 Parametric analysis and summary of dimensionless numbers

Tables 2 and 3 list values of the dimensionless parameters we will examine in
this study, which correspond to typical values one can find in microfluidic
experiments with polymeric solutions and PET and PDMS sheets[53, 58]. The
values of geometrical parameters we consider are of the order of the
dimensions mentioned in the study performed by Mehboudi et al.[53] where the
channel height, width, and length are $h_{0}^{*}=2.5$ $\mu$m, $w^{*}=2000$
m$\mu m$ and $l^{*}=1$ cm, respectively. The typical flow rate through such
conduits is $q^{*}=0.1$ $\mu$L/min, which gives the average flow speed through
the channel being $U^{*}=q^{*}/(2h_{0}^{*}w^{*})=1.667\times 10^{-4}$ m/s. The
elastic sheet that we consider is PDMS with a $E^{*}=1.6$ MPa, Poisson’s ratio
$\nu=0.5$, and thickness $t_{0}^{*}\approx 0.6$ mm, giving an elastic
stiffness $\kappa^{*}=1.97\times 10^{9}$ Pa/m according to the equation for
$t_{0}^{*}/w^{*}\ll 1$ given in Table 1.

We consider three possible fluids whose shear rheology fit the linear sPTT
model. The first example is an aqueous solution of polyvinylpyrolidone (PVP).
At 8 $\%$ weight and molecular weight 360 kDa, the density of the solution is
$\rho^{*}\approx 1000$ kg/$\text{m}^{3}$, zero shear viscosity $\eta^{*}\sim
O(0.1)$ Pa.s and relaxation time $\lambda^{*}\sim O(1)$ ms. Similarly, a
solution of actatic polystyrene (a-PS, 0.008$\%$ wt., 6 MDa) in dioctyl
phthalate (DOP) exhibits a density $\rho^{*}\approx 1000$ kg/$\text{m}^{3}$,
zero shear viscosity $\eta^{*}\sim O(0.1)$ Pa.s and relaxation time
$\lambda^{*}\sim O(100)$ ms. Lastly, a solution of actatic polystyrene (a-PS,
0.14 $\%$ wt., 6.9 MDa) in tricresyl phosphate (TCP) has density
$\rho^{*}\approx 1000$ kg/$\text{m}^{3}$, zero shear viscosity $\eta^{*}\sim
O(0.1)$ Pa.s and relaxation time $\lambda^{*}\sim O(10)$ ms. All three
examples exhibit Newtonian behavior at low shear rates, but shear thinning
with a power law exponent $n\approx 1/3$ at high shear rates consistent with
the linear sPTT model[59, 60].

For power-law fluids, multiple experimental studies have been performed to
analyse flow dynamics in complex processes. [24, 61, 62]. These studies
generally deal with fluids like $0.1\%$ Xanthan gum (molecular weight -
$10^{6}$ g/mol) and $0.1\%$ polyethylene oxide(PEO) (molar weight- $4\times
10^{6}$ g/mol).

We see from the Tables 2 and 3 that the proposed geometry, elastic sheet, and
fluid properties satisfy the following conditions – (a) $\delta\ll 1$, $\delta
Re\ll 1$, and $\delta Wi\ll 1$ for the sPTT fluid, where
$Re=\rho^{*}U^{*}h_{0}^{*}/\eta^{*}$ is the channel Reynolds number that
dictates the inertial contribution; (b) $\delta\ll 1$ and $\delta Re_{eff}\ll
1$ for the generalized Newtonian fluid, where
$Re_{eff}=\rho^{*}U^{*}h_{0}^{*}/\eta^{*}_{eff}$ is the effective Reynolds
number. Under these limits, one can perform a lubrication analysis to
determine how non-Newtonian rheology affects the sheet’s peeling front and
time during startup, pressure controlled flow. Sections 3 and 4 will discuss
the lubrication model and results.

Variable |  Name |  Definition |  Order of Magnitude
---|---|---|---
$\delta$ |  Inverse channel aspect ratio |  $h_{0}^{*}/l^{*}$ |  $O(10^{-4}-10^{-3})$
$Re$ |  Reynolds number |  $\rho^{*}U^{*}h_{0}^{*}/\eta^{*}$ |  $O(10^{-2}-1)$
$Wi$ |  Weissenberg number |  $\lambda^{*}U^{*}/h_{0}^{*}$ |  $O(1-10)$
$\varepsilon$ |  Elongation parameter |  – |  $O(10^{-2})$
$\beta$ |  FSI parameter |  $P_{lub}^{*}/(h_{0}^{*}\kappa^{*})$ |  $O(10^{-2}-10)$
Table 2: Dimensionless numbers for channel geometry, channel elasticity, and sPTT fluid along with typical values for experiments. The lubrication pressure is $P_{lub}^{*}=\frac{\eta^{*}U^{*}}{\delta h_{0}^{*}}$. Variable |  Name |  Definition |  Order of Magnitude
---|---|---|---
$\delta$ |  Inverse channel aspect ratio |  $h_{0}^{*}/l^{*}$ |  $O(10^{-3}-10^{-2})$
$\eta_{eff}^{*}$ |  Effective viscosity |  $K^{*}|\dot{\gamma^{*}}|^{n-1}$ |  $O(10^{-3})-O(0.1)$ $Pa.s^{n}$
$n$ |  Power law index |  – |  $0.3-0.8$
$\beta$ |  FSI Parameter |  $P_{lub}^{*}/(h_{0}^{*}\kappa^{*})$ |  $O(10^{-2}-10)$
Table 3: Key parameters for a power-law fluid along with typical values for
experiments. The characteristic shear rate is
$\dot{\gamma^{*}}=U^{*}/h_{0}^{*}$ and the lubrication pressure is
$P_{lub}^{*}=\frac{\eta_{eff}^{*}U^{*}}{\delta h_{0}^{*}}$.

##  3  Lubrication analysis – sPTT fluid

###  3.1 Lubrication Equations

We invoke the lubrication approximation to solve for the fluid flow and the
height profile for the deformable channel when the fluid satisfies the sPTT
constitutive relationship. The lubrication approximation assumes that the
channel is long and slender such that $\delta\ll 1$ and $\delta Re\ll 1$. This
approximation allows us to neglect the acceleration terms (both Eulerian as
well as convective) in the momentum equation, so the flow is nearly
unidirectional:

$-\nabla^{*}p^{*}+\nabla^{*}\cdot\bm{\tau}^{*}=0$ (11)

Let us proceed with non-dimensionalization. We scale the lengths in in the
$x_{1}^{*}$ and $x_{2}^{*}$ direction by $l^{*}$ and $h_{0}^{*}$ respectively,
and the velocity in the $x_{1}^{*}$ direction by $U^{*}$. Using equation of
continuity, the velocity scale in the $x_{2}^{*}$ is revealed to be $\delta
U^{*}$. Time is scaled by the average residence time
$t^{*}_{conv}=l^{*}/U^{*}$ in the conduit. The pressure is scaled by a
lubrication pressure scale $P^{*}_{lub}$, and the stresses are scaled by
$\delta P^{*}_{lub}$. Depending on the problem at hand (pressure controlled or
flowrate controlled), one of the two scales $P^{*}_{lub}$ or $U^{*}$ will be
specified, with the other quantity obtained through
$P^{*}_{lub}=\eta^{*}U^{*}/(\delta h_{0}^{*})$. Below are the definitions of
the non-dimensional variables:

$\displaystyle x=\frac{x_{1}^{*}}{l^{*}};\quad
y=\frac{x_{2}^{*}}{h_{0}^{*}};\quad h=\frac{h^{*}}{h_{0}^{*}}$ (12a)
$\displaystyle u=\frac{v_{1}^{*}}{U^{*}};\quad v=\frac{v_{2}^{*}}{\delta
U^{*}};\quad t=\frac{U^{*}t^{*}}{l^{*}}$ (12b) $\displaystyle
p=\frac{p^{*}}{P^{*}_{lub}};\quad\bm{\tau}=\frac{\bm{\tau}^{*}}{\delta
P^{*}_{lub}};$ (12c) $\displaystyle P^{*}_{lub}=\frac{\eta^{*}U^{*}}{\delta
h_{0}^{*}}$ (12d)

Using the above definitions, we write the Cauchy momentum equations Eq.(11) in
non-dimensional form and discard terms of $O(\delta)$ or higher. The equations
become:

$\displaystyle\frac{\partial\tau_{xy}}{\partial y}=\frac{\partial p}{\partial
x}$ (13a) $\displaystyle\frac{\partial p}{\partial y}=0$ (13b)

To solve for the flow field in the channel, one needs to relate the shear
stress $\tau_{xy}$ to the velocity field. We proceed to inspect the linear
sPTT model in Eq.(6), Eq.(7), and Eq.(8b). In dimensionless form, the $xx$,
$yy$, and $xy$ components of the stress equations are written below,
neglecting terms $O(\delta)$ and higher as well as terms $O(\delta Wi)$ as
stated previously:

$\displaystyle f(\tau_{xx}+\tau_{yy})\tau_{xx}-2Wi\left[\frac{\partial
u}{\partial y}\tau_{xy}\right]=0$ (14a) $\displaystyle\tau_{yy}=0$ (14b)
$\displaystyle f(\tau_{xx}+\tau_{yy})\tau_{xy}-Wi\left[\frac{\partial
u}{\partial y}\tau_{yy}\right]=\frac{\partial u}{\partial y}$ (14c)

Because we neglected the $O(\delta Wi)$ terms, there is no transient term in
the constitutive viscoelastic relationship. This means that one can treat the
flow field as quasi-steady over timescales comparable to the residence time in
the channel as well as the deformation of the channel.

On inspecting the equation for function $f$ in Eq.(8b), one obtains:

$f(\tau_{xx}+\tau_{yy})=1+{\color[rgb]{0,0,1}\varepsilon}Wi(\tau_{xx}+\tau_{yy})$
(15)

Furthermore, dividing Eq.(14a) by Eq.(14c) yields a relationship between the
normal ($\tau_{xx}$) and shear stress ($\tau_{xy}$):

$\tau_{xx}=2Wi\left(\tau_{xy}\right)^{2}$ (16)

The above relationship illustrates the nonlinear characteristics of the
rheological model lucidly; there is nonzero normal stress present in the flow
$\tau_{xx}$, even though there is no normal strain rate $\frac{\partial
u}{\partial x}=0$. This normal stress does not directly contribute a force on
the channel wall – however, it indirectly alters the pressure distribution by
affecting the shear thinning of the fluid. We also remark that in the limit of
$Wi\to 0$, the Newtonian limit is retrieved and the normal stresses are zero.
Substituting expressions Eq.(16) and Eq.(15) into Eq.(14c) gives the
relationship between the shear stress and the velocity field:

$\left(1+2\varepsilon Wi^{2}\tau_{xy}^{2}\right)\tau_{xy}=\frac{\partial
u}{\partial y}$ (17)

The above captures some key features observed in the rheology of polymeric
solutions. At low shear rates ($Wi\ll 1$), the shear stress depends linearly
on the shear rate: $\tau_{xy}\sim\frac{\partial u}{\partial y}$. At large
shear rates ($Wi\gg 1$), the shear stress $\tau_{xy}\sim(\frac{\partial
u}{\partial y})^{1/3}$. This is equivalent to a shear thinning fluid with
power-law exponent $n=\frac{1}{3}$. We also note that the shear thinning is
direct consequence of nonlinearity of the rheological model. Had
$f(\tau_{xx}+\tau_{yy})=1$, the model would have been linear and a linear
relationship between shear stress and strain rate from Eq. (14c) been obtained
with no shear thinning.

Using Eq.(17), we can now solve for the flow field in the channel. We first
integrate the $x$-momentum equation Eq.(13a) over the $y$ direction and
enforce the symmetry condition at the channel center ($y=0$). This procedure
yields $\tau_{xy}=\frac{\partial p}{\partial x}y$. Substituting this
expression into the stress relationship Eq.(17) gives the shear rate in terms
of the pressure gradient.

$\frac{\partial u}{\partial y}=\left(\frac{\partial p}{\partial
x}\right)y+2\varepsilon Wi^{2}\left(\frac{\partial p}{\partial
x}\right)^{3}y^{3}$ (18)

We integrate the above expression and enforce the no-slip boundary condition
$u(y=h)=0$ at the wall. This procedure yields the velocity field in the
channel:

$u=\left(\frac{\partial p}{\partial
x}\right)\left(\frac{y^{2}-h^{2}}{2}\right)+\varepsilon
Wi^{2}\left(\frac{\partial p}{\partial
x}\right)^{3}\left(\frac{y^{4}-h^{4}}{2}\right)$ (19)

The above equation matches the expression given in [63]. On further inspecting
this velocity profile, we can make some comments. When the value of
$\varepsilon Wi^{2}=0$, the velocity profile collapses to a Newtonian fluid.
Now, this could mean that either $Wi=0$ or $\varepsilon=0$. The former case
deals with a fluid that does not have any elasticity thereby giving us a
Newtonian profile. The latter case, however, pertains to an Oldroyd-B (Boger)
fluid where there are no shear thinning effects. In this case, we would expect
the normal stresses to not contribute heavily to the velocity distribution up
to leading order. We point out an important argument about the scaling of the
stresses. In the scaling given in Eq.(12c), we take the scales for the
stresses to be the same for all the three components. This assumption holds
true for extremely shallow channels– $\delta\sim O(10^{-4})$. Under this
limit, the terms of order $\epsilon$ and higher can be neglected to give us an
accurate description of the viscous stresses. However, we do note that if the
slenderness becomes larger $\delta\sim O(0.01)-O(1)$, the normal stress terms
would have considerable effect. In such analyses, the scaling is different for
the three stress components and one can solve for the flow profile using a
regular perturbation expansion in $\delta Wi$ using techniques such as the
reciprocal theorem [64, 65].

We now determine how the channel height evolves over time. We start with the
continuity equation:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (20)

and integrate it over the height of the channel. Enforcing the kinematic
boundary condition at the sheet interface $v(y=h)=\frac{\partial h}{\partial
t}$ yields:

$\frac{\partial h}{\partial t}=\frac{\partial}{\partial
x}\left(\frac{h^{3}}{3}\left(\frac{\partial p}{\partial
x}\right)+\frac{2\varepsilon Wi^{2}h^{5}}{5}\left(\frac{\partial p}{\partial
x}\right)^{3}\right)$ (21)

This is the nonlinear partial differential equation that governs the height of
the conduit. As a consistency check, we remark that in the limit of
Weissenberg number going to zero $Wi\to 0$, which is the Newtonian limit, we
retrieve the same equation which was earlier reportedly used for modeling of
peeling of linearly elastic sheets conveying Newtonian fluids (see Eq.
($4.11$) in [66]) In order to close this equation, we note that the
relationship between the local pressure and channel height is given by:

$h=1+\beta p$ (22)

where $\beta$ is the FSI parameter discussed in Section 2.2. The following
sections will solve the PDE for different cases and describe how non-Newtonian
fluid rheology alters the start up effects for channel deformation.

###  3.2 Numerical and analytical solution – flow rate controlled case

For a fixed flowrate per unit width $Q^{*}$, we will set the characteristic
velocity and lubrication pressure scales for the non-dimensional variables in
eqns Eq.(12a)-Eq.(12c) as follows:

$U^{*}=\frac{Q^{*}}{2h_{0}^{*}};\qquad
P^{*}_{lub}=\frac{\eta^{*}Q^{*}}{2\delta h_{0}^{*2}}$ (23)

Using Table 2, the definitions for the Weissenberg number and the FSI
parameter for this problem are:

$Wi=\frac{\lambda^{*}Q^{*}}{2h_{0}^{*2}};\qquad\beta=\frac{\eta^{*}Q^{*}}{2\delta\kappa^{*}h_{0}^{*3}}$
(24)

Below, we will write the differential equation for the non-dimensional
pressure profile along the conduit. In non-dimensional form, the volumetric
flowrate per unit width is $q=\int_{-h}^{h}udy$, where $u$ is the velocity
profile. By the definition of our non-dimensionalization, the flowrate is also
$q=1$. Plugging the velocity profile in Eq.(19) and integrating over $y$, we
obtain.

${q}=1=-\frac{2h^{3}}{3}\left(\frac{\partial p}{\partial
x}\right)-\frac{4\varepsilon Wi^{2}h^{5}}{5}\left(\frac{\partial p}{\partial
x}\right)^{3}$ (25)

The above is the nonlinear ODE for pressure for a fixed flow rate at steady
state. If one specifies the pressure at either the inlet or outlet, one can
solve for the pressure profile in the conduit. Furthermore, on using the
relationship Eq.(22) between height and pressure, one can solve for the steady
deformation profile $h(x)$.

Fig.2 and 3 plot the conduit’s steady state deformation profile for different
values of Weissenberg number ($Wi$) and FSI parameter ($\beta$). The curves
are numerical solutions to the above ODE, while dots are the analytical
solution derived in the Supporting Info (Appendix A). Fig.2 plots steady state
deformation profiles for a fixed value of $\beta$ while varying the values of
$Wi$. This case represents a situation where one fixes the flowrate ($Q^{*}$),
zero shear viscosity ($\eta^{*}$), and sheet’s elastic properties
($\kappa^{*}$), but varies the fluid’s relaxation time ($\lambda^{*}$). As the
fluid’s viscoelasticity increases ($Wi$ increases) for a fixed flowrate, the
deformation decreases. Note – a similar analysis has been discussed before for
a deformable channel using an sPTT fluid and a Newtonian fluid [38, 20]. The
main reason cited for this reduction in deformation is the increased shear
thinning nature of the fluid.

(a)  Figure 2: Steady deformation profile of the elastic sheet for different
values of Weissenberg number ($Wi$) at fixed $\beta=0.4$ in an sPTT fluid. The
flow rate is $q=1$. The sPTT elongation parameter $\varepsilon=0.01$. The
pressure at the outlet is kept at a normalized value of 0. Definitions of the
Weissenberg number and FSI parameter are given in Eq.(24).

Fig.3 plots steady state deformation profiles when one fixes the ratio
$\alpha=Wi/\beta$, but varies $Wi$. This situation represents the case when
one fixes the sheet’s elastic properties ($\kappa^{*}$) and fluid’s
rheological properties (viscosity $\eta^{*}$ and relaxation time
$\lambda^{*}$), but one varies the flowrate ($Q^{*}$). For a fixed value of
$\alpha=1$, we notice that increasing the Weissenberg number (proportional to
flow rate) increases the maximum deformation of the channel, although the
increase is nonlinear.

(a)  Figure 3: Steady deformation profile of the elastic sheet in an sPTT
fluid for different values of Weissenberg number ($Wi$) and fixed
$\alpha=Wi/\beta=1$. The sPTT elongation parameter $\varepsilon=0.01$. The
pressure at the outlet is kept at a normalized value of 0. Definitions of the
Weissenberg number and FSI parameter are given in Eq.(24).

###  3.3 Numerical solution – pressure controlled case

(a)

(b)

(c)

Figure 4:  Numerical solution for transient pressure profiles at different
Weissenberg numbers $\left(Wi\right)$. For each Weissenberg number, the
pressure profiles are plotted at five different time instances
$t_{1}=0.001\beta;$ $t_{2}=0.01\beta;$ $t_{3}=0.05\beta;$ $t_{4}=0.1\beta;$
$t_{5}=2.5\beta$. The FSI parameter is $\beta=0.4$ for all plots. The
definitions of the FSI parameter $\beta$, Weissenberg number $Wi$, and non-
dimensional time $t$ are given by Eq. (26) and Eq.(27).

####  3.3.1 Definitions

We will examine an initially static conduit subject to a specified pressure
drop $\Delta P^{*}$ at time $t^{*}=0$. The inlet is clamped at pressure
$\Delta P^{*}$, while the outlet is clamped at zero pressure. The goal is to
determine how the pressure profile evolves over time. We will set the
lubrication pressure and velocity scales for the nondimensional variables in
eqns Eq.(12a)-Eq.(12c) to be:

$P^{*}_{lub}=\Delta P^{*};\qquad U^{*}=\frac{\delta\Delta
P^{*}h_{0}^{*}}{\eta^{*}}$ (26)

Using Table 2, the definitions for the Weissenberg number and the FSI
parameter are:

$Wi=\frac{\delta\lambda^{*}\Delta P^{*}}{\eta^{*}}\qquad\beta=\frac{\Delta
P^{*}}{\kappa^{*}h_{0}^{*}}$ (27)

Similarly, the definition of the non-dimensional time is:

$t=\frac{t^{*}U^{*}}{l^{*}}=\frac{t^{*}\delta^{2}\Delta P^{*}}{\eta^{*}}$ (28)

Given these definitions, the non-dimensional pressure distribution in the
conduit will satisfy PDE Eq.(21) with the condition $h=1+\beta p$. The initial
pressure is $p(x,t=0)=0$ with boundary conditions $p(x=0,t)=1$ and
$p(x=1,t)=0$. We propose two methods of solving this problem where the first
is the complete numerical solution of Eq.(21) and the second is a similarity
solution that aids in our understanding of the peeling rates.

####  3.3.2 Numerical solution

The PDE for pressure (Eq.(21)) is a nonlinear partial differential equation
with Dirichlet boundary conditions. These class of equations can be solved
using a nonlinear finite difference scheme [67, 68]. We divide the domain
$(x,t)$ into $N$ points and $M$ time steps. Each spatial step has a size
$\Delta x$ and each time step has a size $\Delta t$. We use a central
difference scheme in space and a backward implicit difference scheme in time.
On discretizing the equations, we get $N$ nonlinear algebraic equations for
the pressure values at the next time step ($p_{i}^{j+1}$) in terms of the
pressure values at the current time step ($p_{i}^{j}$), where $i=1,...,N$ are
the indices for the spatial points and $j$ is the index for time. These
nonlinear equations are solved iteratively at each time step using Newton’s
method. We typically use time steps of $\Delta t=10^{-8}-10^{-6}$ and grid
size $\Delta x$ between $0.01$ and $0.0223$. Our code has been tested for
convergence.

Figure 4 plots typical pressure profiles obtained from our numeric scheme at
different values of Weissenberg number ($Wi$). As mentioned previously, we set
the initial pressure to be the static pressure $p(x,t=0)=0$, and we specify
the inlet and outlet to be at $p(x=0,t)=1$ and $p(x=1,t)=0$. We observe two
major trends. First, we see that the pressure profile behaves as a front that
propagates towards a steady state. Second, we see that increasing the
Weissenberg number ($Wi$) hastens the time to reach steady state.

Figure 5: Front width $x_{f}$ plotted as a function of $t$ for values of $Wi$.
The value of $\beta$ is fixed at $0.4$. The definitions of the FSI parameter
$\beta$, Weissenberg number $Wi$, and non-dimensional time $t$ are given by
Eq. (26) and Eq.(27) .

To be more quantitative, we define the front width $x_{f}$ as coordinate at
which the pressure drops below $10^{-3}$. In Fig.5, we plot how $x_{f}$ varies
in time for different $Wi$. We see an interesting trend. For low $Wi$
($0<Wi<0.1$), the front appears to increase as $x_{f}\sim t^{1/2}$, while for
high $Wi$ ($Wi\sim O(10)$) the front appears to increase as $x_{f}\sim
t^{1/4}$. At first glance, this data seems at odds with the observation that
viscoelasticity hastens the time to reach equilibrium. However, we see from
the plot that the prefactor in front of the the power law is an order of
magnitude larger for strongly viscoelastic fluids. We will obtain analytical
estimates for such scaling relationships in the next section.

###  3.4  Scaling and similarity solution for pressure controlled flow

Let us look at PDE Eq.(21) with the substitution $h=1+\beta p$:

$\displaystyle\beta\frac{\partial p}{\partial t}=\frac{\partial}{\partial
x}\left(\frac{1}{3}h^{3}\frac{\partial p}{\partial x}+\frac{2\varepsilon
Wi^{2}}{5}h^{5}\left(\frac{\partial p}{\partial x}\right)^{3}\right);\qquad
h=1+\beta p$ (29a) $\displaystyle p(x,t=0)=0;\quad p(x=0,t)=1;\quad
p(x=1,t)=0$ (29b)

Under specific conditions, the PDE admits a similarity solution where one can
obtain scaling relationships for the front width ($x_{f}$) and the peeling
time ($t_{f}$) – i.e, the time it takes for the front to reach the conduit
exit. Tables 4 and 5 summarize the scaling laws for different parameter
regimes. Below describes how one obtains such expressions.

|  Newtonian regime ($\varepsilon Wi^{2}\ll 1$) |  non-Newtonian regime ($\varepsilon Wi^{2}\gg 1$)
---|---|---
Non-dimensional peeling time ($t_{f}$) |  $\beta$ |  $\beta/(\varepsilon Wi^{2})$
Non-dimensional front width ($x_{f}$) |  $\sqrt{t/t_{f}}$ |  $\left(t/t_{f}\right)^{1/4}$
Dimensional peeling time ($t_{f}^{*}$) |  $t_{FSI}^{*}$ |  $t_{FSI}^{*}/(\varepsilon Wi^{2})$
Dimensional front width ($x_{f}^{*}$) |  $l^{*}\sqrt{t^{*}/t_{f}^{*}}$ |  $l^{*}\left(t^{*}/t_{f}^{*}\right)^{1/4}$
Table 4: Scaling relationships for front width and peeling time for an sPTT fluid with moderate conduit deformations $\beta\sim O(1)$. The definitions of the FSI parameter $\beta$, Weissenberg number $Wi$, and non-dimensional time $t$ are given by Eq.(27) and Eq. (28). The Newtonian FSI timescale in the small deformation limit is is $t_{FSI}^{*}=\frac{\eta^{*}}{\delta^{2}h_{0}^{*}\kappa^{*}}$. |  Newtonian regime ($\varepsilon Wi^{2}\beta^{2}\ll 1$) |  non-Newtonian regime ($\varepsilon Wi^{2}\beta^{2}\gg 1$)
---|---|---
Non-dimensional peeling time ($t_{f}$) |  $\beta^{-2}$ |  $1/(\beta^{4}\varepsilon Wi^{2})$
Non-dimensional front width ($x_{f}$) |  $\sqrt{t/t_{f}}$ |  $\left(t/t_{f}\right)^{1/4}$
Dimensional peeling time ($t_{f}^{*}$) |  $t_{FSI}^{*}/\beta^{3}$ |  $t_{FSI}^{*}/(\beta^{5}\varepsilon Wi^{2})$
Dimensional front width ($x_{f}^{*}$) |  $l^{*}\sqrt{t^{*}/t_{f}^{*}}$ |  $l^{*}\left(t^{*}/t_{f}^{*}\right)^{1/4}$
Table 5: Scaling relationships for front width and peeling time for strong
conduit deformations $\beta\gg 1$. The definitions of the FSI parameter
$\beta$, Weissenberg number $Wi$, and non-dimensional time $t$ are given by
Eq.(27) and Eq. (28). The Newtonian FSI timescale in the small deformation
limit is $t_{FSI}^{*}=\frac{\eta^{*}}{\delta^{2}h_{0}^{*}\kappa^{*}}$

####  3.4.1 Moderate conduit deformation ($\beta\sim O(1)$ or smaller)

When the conduit’s deformations are moderate ($\beta\sim O(1)$), the relative
contribution of the two terms on the right hand side of Eq.(29a) depend on the
quantity $\varepsilon Wi^{2}$. The conduit’s FSI is dominated by Newtonian
fluid rheology when $\varepsilon Wi^{2}\ll 1$, while the FSI dominated by non-
Newtonian rheology when $\varepsilon Wi^{2}\gg 1$. Below describes similarity
solutions under these regimes.

Newtonian FSI: Let us inspect the case of Newtonian FSI – i.e., $\varepsilon
Wi^{2}\ll 1$. Upon taking this limit, the PDE Eq.(29a) becomes:

$\beta\frac{\partial p}{\partial t}=\frac{\partial}{\partial
x}\left(\frac{1}{3}h^{3}\frac{\partial p}{\partial x}\right);\qquad h=1+\beta
p\\\ $ (30)

with the same boundary conditions as before (Eq.(29b)). For early times
$t\ll\beta$, the height profile admits a similarity solution. The profile
looks like a front that propagates from the conduit’s inlet to the exit, with
the front’s width scaling as $x_{f}\sim\sqrt{t/\beta}$. The front reaches the
end when $x_{f}\sim O(1)$, which gives an expression for the peeling time as
$t_{f}\sim\beta$. Using these ideas, we write the form of the solution as
follows:

Figure 6:  Similarity solution for Newtonian FSI ($\varepsilon Wi^{2}\ll 1$)
at moderate sheet deformations ($\beta\sim O(1))$ or smaller). The similarity
solution is plotted against the numerical finite difference solution at
different time instances. For $\beta=0.01$, the time instances plotted are
$t=5\times 10^{-5},5\times 10^{-4},$ and $10^{-3}$. For $\beta=0.5$, the time
instances are $t=5\times 10^{-5},5\times 10^{-4},10^{-3}$, and $2\times
10^{-2}$. For both plots, $Wi=0.01$ and $\varepsilon=0.01$. For $\beta=0.01$,
the leading order analytical solution for $\beta\ll 1$ is also plotted
($p=\text{erfc}(\eta/\sqrt{2})$). The definitions of the FSI parameter
$\beta$, Weissenberg number $Wi$, and non-dimensional time $t$ are given by
Eq.(27) and Eq. (28).

$p=f(\eta);\qquad\eta\equiv\frac{x\sqrt{3\beta}}{\sqrt{2t}}$ (31)

The above transformation converts the PDE Eq.(33) into an ODE:

$\displaystyle\frac{d}{d\eta}\left(h^{3}\frac{df}{d\eta}\right)+\eta\frac{df}{d\eta}=0;\qquad
h=1+\beta f$ (32a) $\displaystyle f(0)=1;\quad f(\infty)=0$ (32b)

This ODE has an analytical solution for $\beta\ll 1$:
$f=\text{erfc}(\eta/\sqrt{2})$. For $\beta\sim O(1)$, one can numerically
solve the ODE using a shooting method, where one applies the boundary
condition at infinity at a sufficiently large value of $\eta$ (typically
$\eta>10$). We increase $\eta$ up to a point where the results become
insensitive to the value of $\eta$. At this value, the similarity plot matches
the numerical solution fairly well. Fig.6 plots the similarity solution versus
the fully resolved finite difference solution from the previous section. For
the conditions stated – i.e., Newtonian fluid ($\varepsilon Wi^{2}\ll 1$) and
early times ($t\ll\beta$), the similarity solution captures the deformation
profile well.

Non-Newtonian FSI: When the FSI is dominated by the non-Newtonian rheology,
$\varepsilon Wi^{2}\gg 1$. In this limit, the PDE Eq.(29a) for the pressure
becomes:

$\beta\frac{\partial p}{\partial t}=\frac{\partial}{\partial
x}\left(\frac{2\varepsilon Wi^{2}}{5}h^{5}\left(\frac{\partial p}{\partial
x}\right)^{3}\right);\qquad h=1+\beta p\\\ $ (33)

with the same boundary conditions as before (Eq.(29b)). For early times
$t\ll\beta/(\varepsilon Wi^{2})$, the height profile admits a similarity
solution. The profile propagates as a front with thickness scaling as
$x_{f}\sim\left(\varepsilon Wi^{2}t/\beta\right)^{1/4}$. The front reaches the
end when $x_{f}\sim O(1)$, which gives the peeling time as
$t_{f}\sim\beta/(\varepsilon Wi^{2})$. Using these ideas, we write the form of
the solution:

$p=f(\eta);\qquad\eta\equiv\frac{x(5\beta)^{1/4}}{(8t\varepsilon
Wi^{2})^{1/4}}$ (34)

The above transformation converts the PDE (33) into an ODE:

$\displaystyle\frac{d}{d\eta}\left(h^{5}\left(\frac{df}{d\eta}\right)^{3}\right)+\eta\frac{df}{d\eta}=0;\qquad
h=1+\beta f$ (35a) $\displaystyle f(0)=1;\quad f(\infty)=0$ (35b)

This ODE is different than the Newtonian ODE in that it is not well defined
for the entire domain $\eta\in[0,\infty)$. At a particular point $\eta_{max}$,
the pressure gradient will switch sign:

$\frac{df}{d\eta}(\eta_{max})=0$

Clearly, this behavior is aphysical so the similarity ODE can only be defined
in the region $\eta\in[0,\eta_{max}]$. What is physically happening is that
around $\eta=\eta_{max}$, the pressure gradients are so small that the
Newtonian contributions to the FSI start to matter. Because the pressure is so
small in this region and the region $\eta>\eta_{\max}$, the details are
unimportant for all practical purposes. Thus, one can reformulate the above
ODE as follows:

(a)

(b)

(c)

Figure 7: Similarity solution for viscoelastic FSI ($\varepsilon Wi^{2}\gg 1$)
at moderate sheet deformations ($\beta\sim O(1))$ or smaller). (a)
$\beta=0.005$. The similarity solution is plotted against the numerical finite
difference solution at different time instances $t=5\times 10^{-7},10^{-6}$,
and $10^{-5}$. The red curve is the leading order analytical solution for
$\beta\ll 1$ (Eq. (37)). (b) $\beta=0.1$. The similarity solution is plotted
against the numerical solution at times $t=1\times 10^{-6},5\times 10^{-6}$,
and $10^{-4}$. (c) $\beta=0.2$. The similarity solution is plotted against the
numerical solution at times $t=1\times 10^{-6},5\times 10^{-6}$, $10^{-4}$,
and $10^{-3}$. For all plots, $Wi=20$ and $\varepsilon=0.01$. The definitions
of the FSI parameter $\beta$, Weissenberg number $Wi$, and non-dimensional
time $t$ are given by Eq.(27) and Eq. (28)

$\displaystyle\frac{d}{d\eta}\left(h^{5}\left(\frac{df}{d\eta}\right)^{3}\right)+\eta\frac{df}{d\eta}=0;\qquad
h=1+\beta f$ (36a) $\displaystyle f(0)=1;\quad f^{\prime}(\eta_{max})=0;\qquad
f(\eta_{max})=0$ (36b)

where we note that the solution will have a small error in a thin region
around the endpoint. For the case of small deformations ($\beta\ll 1$), the
above ODE has an analytical solution:

$f=1-\frac{\eta_{max}^{2}}{\sqrt{3}}\int_{0}^{\eta/\eta_{max}}\sqrt{1-z^{2}}dz\qquad\eta_{max}=\left(\frac{4\sqrt{3}}{\pi}\right)^{1/2}\qquad\beta\ll
1$ (37)

For values of $\beta\sim O(1)$, one will have to numerically solve the above
ODE. The way we approach this task is to perform a nested shooting method. We
first define a function that takes an input $\eta_{max}$ and uses a shooting
method to compute the solution to the above ODE with boundary conditions
$f(0)=1$ and $f^{\prime}(\eta_{max})=0$. The function outputs the value at the
endpoint $f(\eta_{max})$. One performs root finding on the output of this
function until one obtains the value of $\eta_{max}$ at which
$f(\eta_{max})=0$. Fig.7 shows the solution to the similarity ODE and compares
it to the fully resolved finite difference solution described previously. This
figure encapsulates the similarity scaling discussed in Table 4.

Fig.8 summarizes the scaling behavior observed in the numerical simulations
for the peeling time at different Weissenberg numbers at moderate conduit
deformations ($\beta\sim O(1)$ or smaller). The figure clearly demarcates the
two regimes of a strongly Newtonian fluid (where $t_{f}\sim\beta$) and a
strongly viscoelastic fluid (where $t_{f}\sim\beta/(\varepsilon Wi^{2})$). At
moderate viscoelasticity ($\varepsilon Wi^{2}\sim O(0.1)$), neither scaling
regime holds and one must resort to a full numerical solution. The next
section will discuss the regime of strong conduit deformations ($\beta\gg 1$).

(a)  Figure 8:  Summary of peeling time for moderate conduit deformations
$(\beta\sim O(1)$ or smaller). The front width ($x_{f}$) is the coordinate at
which the scaled pressure $p$ drops to under $10^{-3}$. The peeling time is
defined as the time at which $x_{f}=1$.The definitions of the FSI parameter
$\beta$, Weissenberg number $Wi$, and non-dimensional time $t$ are given by
Eq.(27) and Eq. (28)

####  3.4.2 Strong conduit deformations ($\beta\gg 1$)

When the conduit’s deformation is strong, the conduit height scales as
$h\sim\beta\gg 1$ rather than $h\sim O(1)$ as stated in the previous section.
This fact alters the relative magnitude of the Newtonian and viscoelastic
terms in the pressure PDE eq (29a). Newtonian FSI now occurs when $\varepsilon
Wi^{2}\beta^{2}\ll 1$, while viscoelsatic FSI occurs in the opposite limit
$\varepsilon Wi^{2}\beta^{2}\gg 1$. Below is a brief summary of the similarity
solutions in both limits.

Newtonian FSI: Under the Newtonian limit $\varepsilon Wi^{2}\beta^{2}\ll 1$,
the limiting PDE for the pressure profile is the same as eqn (33). However,
since the conduit’s height scales as $h\sim\beta$, the peeling front width and
the peeling time are different. The front’s width scales as
$x_{f}\sim\sqrt{t\beta^{2}}$. The peeling time is determined when $x_{f}\sim
O(1)$, which occurs at $t_{f}\sim\beta^{-2}$.

We note that the peeling time in this strong deformation limit is markedly
different than the moderate deformation case discussed previously, where
$t_{f}\sim\beta$. Fig.9 shows numerical simulations of the peeling time for
different values of $\beta$ that clearly demarcates the two regimes. The non-
monotonic trend in peeling time on the conduit deformation $\beta$ is a non-
intuitive finding that arises due to the competition between lubrication
pressure ($p_{lub}\sim\frac{1}{h^{3}}\frac{\partial h}{\partial t}$) and the
sheet’s elastic stresses ($p_{elas}\sim(h-1)\beta$). At large conduit
deformations, the lubrication pressure is greatly attenuated, which gives rise
to the trends above.

(a)  Figure 9: Peeling time versus FSI parameter $\beta$ for the Newtonian FSI
regime ($\varepsilon Wi^{2}\ll 1$ for $\beta\sim O(1)$ or smaller, and
$\varepsilon Wi^{2}\beta^{2}\ll 1$ for $\beta\gg 1$). In the plot, we fix
$Wi=1$ and $\varepsilon=0.01$ and vary $\beta$. The definition of $\beta$ for
this problem is given in Eq.(27).

To obtain a similarity solution in the regime of strong conduit deformation
$(\beta\gg 1)$, we make the following transformation:

$p=f(\eta);\qquad\eta\equiv\frac{x\sqrt{3}}{\sqrt{2t\beta^{2}}}$ (38)

which converts the PDE Eq.(33) into the ODE:

$\displaystyle\frac{d}{d\eta}\left(\tilde{h}^{3}\frac{df}{d\eta}\right)+\eta\frac{df}{d\eta}=0;\qquad\tilde{h}=f+\beta^{-1}$
(39a) $\displaystyle f(0)=1;\quad f(\infty)=0$ (39b)

This equation is not well-defined on the entire domain $\eta\in[0,\infty)$.
There is a specific value of $\eta$ at which the pressure drops to zero:

$f(\eta_{max})=0$

We thus replace the boundary condition at infinity with the one above. The
value of $\eta_{max}$ is obtained from simulation.

Fig.10 plots the similarity solution in the high $\beta$ limit compared to the
finite difference simulations. Overall, we see a reasonable collapse of the
numerical data to the similarity curve, although there is some scatter. We
also note that the similarity solution is qualitatively different than the
moderate deformation case ($\beta\sim O(1)$ in Fig.7). Here, the shape of the
front is concave down, whereas the shape for the moderate deformation case is
concave up.

Figure 10: Similarity solution for Newtonian FSI ($\varepsilon
Wi^{2}\beta^{2}\ll 1$) for strong sheet deformations ($\beta\gg 1)$). The
similarity solution is plotted against the numerical finite difference
solution at different time instances $t=10^{-5},10^{-4},10^{-3}$. The
definitions of the FSI parameter $\beta$, Weissenberg number $Wi$, and non-
dimensional time $t$ are given by Eq.(27) and Eq. (28).

Non-Newtonian FSI: Table 4 lists the scaling relationships for non-Newtonian
FSI in the regime of strong conduit deformation ($\beta\gg 1$). Although one
can develop a similarity solution in this regime, we believe that this
situation might not be realizable using the current model. At high conduit
deformations, the value of $\delta Wi$ will likely no longer be small (likely
$O(0.1)$ or larger), making the upper convected derivative term in Eq. (6)
important. In this situation, the relaxation of the polymer chain will be
comparable to the average transit time in the conduit, giving a non-quasi
steady flow profile [69]. We will not discuss this situation in this
manuscript.

##  4 Lubrication analysis - Power law fluid

###  4.1 Lubrication Equations

The setup is exactly the same as described in section 2.1 except that the non-
Newtonian fluid is now a power law fluid with power law index $n$ and
consistency index $K^{*}$. For an initial gap thickness $h_{0}^{*}$ and
average fluid velocity $U^{*}$, the effective viscosity scale is:

$\eta_{eff}^{*}=K^{*}\left(\frac{U^{*}}{h_{0}^{*}}\right)^{n-1}$ (40)

When performing the lubrication analysis discussed before, we will keep the
same definitions for the non-dimensional lengths, velocities, time, pressure,
and stresses as in Eq.(12a)-Eq.(12c). The relationship between the lubrication
pressure scale $P^{*}_{lub}$ and the velocity scale $U^{*}$ will have a
similar form as Eq.(12d), except that one replaces the viscosity of the fluid
by the effective viscosity $\eta^{*}_{eff}$:

$P^{*}_{lub}=\frac{\eta_{eff}^{*}U^{*}}{\delta h_{0}^{*}}$ (41)

Depending on the problem at hand (pressure controlled or flow controlled), one
of the two scales $P^{*}_{lub}$ or $U^{*}$ will be specified with the other
determined by the above two relationships (Eq. (40), Eq.(41)). The lubrication
approximation will hold if the conduit is long and slender ($\delta\ll 1$) and
the effective Reynolds number $Re_{eff}\ll\delta^{-1}$. Here, the effective
Reynolds number is the Reynolds number using the effective viscosity defined
above: $Re_{eff}=\rho^{*}U^{*}h_{0}^{*}/\eta_{eff}^{*}$.

In non-dimensional form, the lubrication equations for the momentum balance
will take the form:

$\frac{\partial\tau_{xy}}{\partial y}=\frac{\partial p}{\partial
x};\qquad\frac{\partial p}{\partial y}=0$ (42)

where

$\tau_{xy}=\left|\frac{\partial u}{\partial y}\right|^{n-1}\frac{\partial
u}{\partial y}$ (43)

We substitute the above equation Eq. (43) into the first equation of the
momentum balance and integrate to find the velocity field $u$. In doing this
procedure, we apply the boundary conditions that there is no shear stress at
the center of the conduit ($y=0$) and no velocity at the top ($y=h$). This
procedure yields:

$u=\frac{n}{1+n}\left(-\frac{\partial p}{\partial
x}\right)^{1/n}\left(h^{\frac{1+n}{n}}-\left|y\right|^{\frac{1+n}{n}}\right)$
(44)

To derive the equation for the height of the conduit, we examine the
continuity equation $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}=0$ and integrate it over the height of the conduit from $y=0$ to $y=h$.
Noting that $v=\frac{\partial h}{\partial t}$ at the top wall, we obtain:

$\frac{\partial h}{\partial t}=-\frac{n}{2n+1}\frac{\partial}{\partial
x}\left(h^{\frac{2n+1}{n}}\left(-\frac{\partial p}{\partial
x}\right)^{1/n}\right)$ (45)

This equation is the PDE that describes the height distribution in the
conduit. It is coupled with the deformation relationship

$h=1+\beta p$ (46)

where $\beta$ is the channel’s FSI parameter.

Eq.(45) governs the evolution of the elastic deformation of the sheet. The
flow rate controlled case at steady state has been discussed in detail by
Anand et.al.[37]. We will proceed to solve the transient problem for the
pressure controlled dynamics.

###  4.2 Results – Pressure controlled case

####  4.2.1 Definitions

Similar to Section 3.3, we will examine a problem where a stagnant conduit at
zero pressure is subject to a pressure impulse at its inlet at time $t^{*}=0$.
The inlet of the conduit is clamped at pressure $\Delta P^{*}$ while the
outlet is clamped at zero pressure. If we specify $\Delta P^{*}$, this will
set the pressure scale, velocity scale, and effective viscosity scale for the
non-dimensional variables in eqns Eq. (12a)-Eq.(12b):

$P^{*}_{lub}=\Delta P^{*};\qquad U^{*}=h_{0}^{*}\left(\frac{\Delta
P^{*}\delta}{K^{*}}\right)^{1/n}\qquad\eta_{eff}^{*}=K^{*}\left(\frac{\Delta
P^{*}{\delta}}{K^{*}}\right)^{\frac{n-1}{n}}$ (47)

The FSI parameter for this problem is:

$\beta=\frac{\Delta P^{*}}{h_{0}^{*}\kappa^{*}}$ (48)

and the dimensionless time is:

$t=t^{*}\frac{U^{*}}{l^{*}}={\delta}t^{*}\left(\frac{\Delta
P^{*}{\delta}}{K^{*}}\right)^{1/n}$ (49)

In non-dimensional form, the PDE that describes the height evolution of the
channel is given by Eq.(45) subject to the condition that $h=1+\beta p$. The
boundary conditions are that the pressure is $p(x=0,t)=1$ and $p(x=1,t)=0$,
with the initial condition that $p(x,t=0)=0$.

####  4.2.2 Numerical results

In this section, we solve the PDE given by Eq.(45). Our numerical procedure is
the same as discussed in Section 3.3.2 except that we use an strongly implicit
Crank-Nicolson time-stepping scheme which is thoroughly described in reference
[67].

Fig.11 plots the pressure profile at different time instances for a Newtonian
fluid and a moderately shear thinning fluid ($n=0.4$). The results indicate
that shear thinning causes the front to propagate more quickly. The next
section goes through scaling analyses and similarity solutions for this
problem.

(a)

(b)

Figure 11: Numerical solution of transient pressure profile for (a) Newtonian
fluid ($n=1$) and (b) moderately shear thinning fluid ($n=0.4$). The time
instances sampled correspond to $t_{1}=0.001\beta$, $t_{2}=0.01\beta$,
$t_{3}=0.05\beta$, $t_{4}=0.1\beta$, and $t_{5}=2.5\beta$. The FSI parameter
is $\beta=0.4$. Note – the definition of FSI parameter $\beta$ and non-
dimensional time $t$ are given by eqns (48-49).

####  4.2.3 Similarity solution

Let us look at PDE Eq. (45) with the substitution $h=1+\beta p$:

$\displaystyle\beta\frac{\partial p}{\partial
t}=-\frac{n}{2n+1}\frac{\partial}{\partial
x}\left(h^{\frac{2n+1}{n}}\left(-\frac{\partial p}{\partial
x}\right)^{1/n}\right);\qquad h=1+\beta p$ (50a) $\displaystyle
p(x,t=0)=0;\quad p(x=0,t)=1;\quad p(x=1,t)=0$ (50b)

We will follow a similar procedure as Section 3.4 to determine a similarity
solution for moderate conduit deformations ($\beta\sim O(1)$ or smaller) and
large conduit deformations ($\beta\gg 1)$. Table 6 summarizes the scaling
results.

|  $\beta\sim O(1)$ |  $\beta\gg 1$
---|---|---
Non-dimensional peeling time ($t_{f}$) |  $\beta$ |  $\beta^{-\frac{(1+n)}{n}}$
Non-dimensional front width ($x_{f}$) |  $\left(t/t_{f}\right)^{\frac{n}{n+1}}$ |  $\left(t/t_{f}\right)^{\frac{n}{n+1}}$
Dimensional peeling time ($t_{f}^{*}$) |  $\frac{t_{f}}{{\delta}}\left(\frac{\Delta P^{*}{\delta}}{K^{*}}\right)^{-\frac{1}{n}}$ |  $\frac{t_{f}}{{\delta}}\left(\frac{\Delta P^{*}{\delta}}{K^{*}}\right)^{-\frac{1}{n}}$
Dimensional front width ($x_{f}^{*}$) |  $l^{*}\left(t^{*}/t_{f}^{*}\right)^{\frac{n}{n+1}}$ |  $l^{*}\left(t^{*}/t_{f}^{*}\right)^{\frac{n}{n+1}}$
Table 6: Scaling relationships for front width and peeling time for a power
law fluid with moderate channel deformations ($\beta\sim O(1)$ or smaller) and
strong channel deformations ($\beta\gg 1$). The definitions of the FSI
parameter $\beta$ and non-dimensional time are given by Eq. (48) and Eq. (49)

Moderate deformations ($\beta\sim O(1)$ or smaller): For moderate conduit
deformations, the above PDE admits a similarity solution at early times
$t\ll\beta$. The profile propagates as a front with thickness scaling as
$x_{f}\sim\left(t/\beta\right)^{\frac{n}{n+1}}$. The front reaches the end
when $x_{f}\sim O(1)$, which gives the peeling time as $t_{f}\sim\beta$. The
form of the similarity solution is as follows:

$p=f(\eta);\qquad\eta\equiv
x\left(\frac{(2n+1)\beta}{(n+1)t}\right)^{\frac{n}{n+1}}$ (51)

The above transformation converts the PDE (50a) into an ODE:

$\displaystyle-\frac{d}{d\eta}\left(h^{\frac{2n+1}{n}}\left(-\frac{df}{d\eta}\right)^{\frac{1}{n}}\right)+\eta\frac{df}{d\eta}=0;\qquad
h=1+\beta f$ (52a) $\displaystyle f(0)=1;\quad f(\infty)=0$ (52b)

(a)

(b)

Figure 12: Comparison between similarity solution and numerical solution for
power law fluid. (a) Moderate conduit deformations ($\beta\sim O(1)$ or
smaller). For $\beta=0.005$, the time instances plotted for the numerical
solution are $t=10^{-7},10^{-5}$, and $5\times 10^{-4}$. The analytical
solution is the leading order solution for $\beta\ll 1$ as stated in eqn (54).
For $\beta=0.5$, the time instances plotted for the numerical solution are
$t=10^{-7},10^{-5}$, and $1.5\times 10^{-4}$. (b) Strong conduit deformation
($\beta\gg 1$). The plot is for $\beta=10$ and the numerical solutions are
plotted at time instants $t=5\times 10^{-7},5\times 10^{-6},10^{-5}$, and
$10^{-4}$. For all plots, the shear-thinning index is $n=0.5$. The definitions
of the FSI parameter $\beta$ and non-dimensional time are given by Eq. (48)
and Eq. (49).

Again, just like in Section 3.4.1, the above ODE is well defined only in a
finite region $\eta\in[0,\eta_{max}]$. At a particular location
$\eta=\eta_{max}$, the pressure gradient will equal zero and hence other
physics will become important. In actuality, the region around $\eta_{max}$
has a pressure so small that the power-law constitutive relationship likely no
longer holds. Thus, we consider replacing the boundary conditions of the ODE
with the following:

$f(0)=1;\quad f(\eta_{max})=0;\qquad f^{\prime}(\eta_{max})=0$ (53)

where $\eta_{max}$ needs to be found.

For $\beta\ll 1$, the above similarity ODE admits an analytical solution:

$\displaystyle
f=1-\left(\frac{1-n}{2}\eta_{max}^{2}\right)^{\frac{n}{1-n}}\eta_{max}\int_{0}^{\eta/\eta_{max}}\left[1-z^{2}\right]^{\frac{n}{1-n}}dz$
(54a)
$\displaystyle\eta_{max}^{\frac{n+1}{1-n}}=\frac{\left(\frac{2}{1-n}\right)^{\frac{n}{1-n}}}{\int_{0}^{1}\left(1-z^{2}\right)^{\frac{n}{1-n}}dz}$
(54b)

For $\beta\sim O(1)$, we solve the ODE using the nested shooting method
discussed before. The results have been shown in Fig.12.

From Table 6, we note that the front propagation scaling - $x_{f}\sim
t^{\frac{n}{n+1}}$ matches with the scaling law provided by Boyko et.al. [70]
thereby suggesting that we could expect the analysis for power law fluids to
also extend to long, thin walled, elastic shells. This can be owed to the fact
that although the axisymmetric governing equations in both the problems differ
by the addition of circumferential stresses in the tube case, the radially
averaged model predicts a spring like behavior for the pressure drop-
deformation relation.

Large deformation ($\beta\gg 1$): When the conduit deformation is large
$(\beta\gg 1)$, the conduit height behaves as $h\sim\beta\gg 1$, which gives
rise to different scaling for the front width and peeling time. For times
$t\ll\beta^{-(1+n)/n}$ one obtains the scalings in Table 5. The similarity
transformation takes the following form:

$p=f(\eta);\qquad\eta\equiv\frac{x}{\beta}\left(\frac{(2n+1)}{(n+1)t}\right)^{\frac{n}{n+1}}$
(55)

The above transformation converts the PDE (50a) into an ODE:

$\displaystyle-\frac{d}{d\eta}\left(\tilde{h}^{\frac{2n+1}{n}}\left(-\frac{df}{d\eta}\right)^{\frac{1}{n}}\right)+\eta\frac{df}{d\eta}=0;\qquad\tilde{h}=f+\beta^{-1}$
(56a) $\displaystyle f(0)=1;\quad f(\infty)=0$ (56b)

Like the previous section, the equations are not well defined in the entire
domain. We will replace the boundary conditions with the following:

$f(1)=0;\qquad f(\eta_{max})=0$ (57)

where $\eta_{max}$ is determined from simulation.

Fig.12 shows the similarity solution for $\beta=10$ and compares it with the
numerical solution at different time instances. Overall, the simulation data
strongly suggests that similarity holds in the small time limit, with the
solution given by the above ODE.

##  5 Conclusion

In this study, we explore a particular class of interactions between
deformable slender sheets and viscoelastic fluids. When two elastic sheets are
separated by a lubricating layer of fluid, the flow of the fluid can cause the
sheets to peel. This situation is common in many applications in nature and
industry, for example exfoliation of graphene sheets[71], hydraulic fracturing
[72], and cell sorting [73]. This study investigates how non-Newtonian liquids
can alter peeling behavior of the sheets. Specifically, the study quantifies
the startup, transient behavior of the peeling process, characterizing how the
peeling front propagates as well as the peeling time. Two classes of complex
fluids are studied – a simplified Phan-Thien Tanner (sPTT) fluid and a
generalized power law fluid. The sPTT model captures many of the qualitative
features of polymeric solutions such as shear thinning, extensional
thickening, viscoelasticity, and normal stress differences. The second fluid
only exhibits shear thinning, but can demonstrate a much wider range of shear
thinning behavior than the sPTT model.

Our study finds that shear thinning plays a major role in modifying the
peeling characteristics of the elastic sheet. Other aspects of non Newtonian
rheology including nonzero normal stress difference, elasticity, and
nonlinearity are subdued for the system under consideration. To explain this
observation for the s-PTT model we note that $\tau_{yy}=0$ (see Eq. (14b) due
to the lubrication approximation, which makes the hydrodynamic normal load on
the structure arise only from pressure and not due to the deviatoric stress.
On the other hand, because of nonlinearity of the rheological model, we have a
finite $\tau_{xx}$ (see Eq. (16), even though the corresponding strain rate
$\frac{\partial u}{\partial x}=0$. This normal stress works to augment the
shear thinning property of the fluid, as shown in Eq.(17), and ultimately
alters the lubrication pressure. Finally, the assumptions pertaining to small
$\epsilon Wi$ allow us to neglect the elasticity/relaxation time of the fluid
completely.

When compared to a Newtonian fluid with equivalent zero-shear viscosity, a
shear thinning liquid will have its peeling front propagate with a smaller
power law exponent than a Newtonian fluid (e.g., $x_{f}\sim t^{-n/(n+1)}$ for
a power law index $n$, compared to $x_{f}\sim t^{1/2}$ for Newtonian), but
with an order-of-magnitude larger prefactor. The overall consequence of this
result is that the peeling time for the shear thinning fluid is much smaller
than the equivalent Newtonian fluid. For both constitutive relationships
studied, scaling relationships are provided for the peeling front and peeling
time in the Newtonian and non-Newtonian (e.g., power law) regions of the
viscosity curves, under the specific limits of moderate sheet deformation and
strong sheet deformation (see Tables 3-5). Similarity solutions are provided
to give a quantitative description in these limiting regimes, as well as full
numerical solutions in the regimes where the similarity solutions no longer
hold (e.g., when ${\varepsilon}Wi^{2}\sim O(1)$ for the sPTT fluid, where $Wi$
is the Weissenberg number and ${\varepsilon}$ is the elongation parameter). We
hope that these results will give much needed insight for those interested in
manipulating slender, deformable sheet structures with non-Newtonian liquids.

##  6 Acknowledgements

The authors acknowledge funding from the American Chemical Society Petroleum
Research Fund (Grant No. ACS PRF 61266-DNI9).

##  Appendix A Derivation of flow controlled case solution at steady state

From section 3, we know that the height evolution equation is given as
follows:

${q}=-\left(\frac{d{p}}{d{x}}\right)\left(\frac{2h^{3}}{3}\right)-\left(\frac{d{p}}{d{x}}\right)^{3}\left(\frac{{\varepsilon}Wi^{2}}{5}\right)h^{5}$
(58)

If we set the non-dimensional flow rate as $1$, we get the following governing
equation:

$1=-\left(\frac{d{p}}{d{x}}\right)\left(\frac{2h^{3}}{3}\right)-\left(\frac{d{p}}{d{x}}\right)^{3}\left(\frac{{\varepsilon}Wi^{2}}{5}\right)h^{5}$
(59)

From the pressure-deformation relation, we have

${h}=1+\beta{p}$ (60)

which when differentiated, yields

$\frac{d{h}}{d{x}}=\beta\frac{d{p}}{d{x}}\implies\frac{1}{\beta}\frac{d{h}}{d{x}}=\frac{d{p}}{d{x}}$
(61)

On substituting (61) into (59), we get

$1+\frac{1}{\beta}\left(\frac{dh}{dx}\right)\left(\frac{2h^{3}}{3}\right)+\frac{1}{\beta^{3}}\left(\frac{dh}{dx}\right)^{3}\left(\frac{{\varepsilon}Wi^{2}}{5}\right)h^{5}=0$
(62)

We make the following transformation:

${h^{2}}={t_{1}}\implies{2h}\frac{d{h}}{d{x}}=\frac{d{t_{1}}}{d{x}}$ (63)

This transformation converts (62) into

$\left(\frac{d{t_{1}}}{d{x}}\right)^{3}\left(\frac{t_{1}kWi^{2}}{40\beta^{3}}\right)+\left(\frac{d{t_{1}}}{d{x}}\right)\left(\frac{t_{1}}{3\beta}\right)+1=0$
(64)

$\frac{d{t_{1}}}{d{x}}=\theta$ (65)

On choosing the variables ${a}=\frac{t_{1}{\varepsilon}Wi^{2}}{40\beta^{3}}$;
${b}=\frac{t_{1}}{3\beta}$; ${c}=1$

We get the following depressed cubic equation:

$a\theta^{3}+b\theta+c=0$ (66)

If we write the terms as $\frac{b}{a}={p_{1}}$ ; $\frac{c}{a}={q_{1}}$

This equation has a well known solution given by the Cardano’s formula:

$\theta=\left(-\frac{q_{1}}{2}+\left(\frac{q_{1}^{2}}{4}+\frac{p_{1}^{3}}{27}\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}+\left(-\frac{q_{1}}{2}-\left(\frac{q_{1}^{2}}{4}+\frac{p_{1}^{3}}{27}\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}$
(67)

This creates the ordinary differential

$\frac{d{t_{1}}}{d{x}}=\left(-\frac{q}{2}+\left(\frac{q^{2}}{4}+\frac{p^{3}}{27}\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}+\left(-\frac{q}{2}-\left(\frac{q^{2}}{4}+\frac{p^{3}}{27}\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}$
(68)

We perform a simple transformation $p=m$ and $q=\frac{n}{t_{1}}$

$\frac{d{t_{1}}}{d{x}}=\left(-\frac{n}{2t_{1}}+\left(\frac{n^{2}}{4t_{1}^{2}}+\frac{m^{3}}{27}\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}+\left(-\frac{n}{2t_{1}}-\left(\frac{n^{2}}{4t_{1}^{2}}+\frac{m^{3}}{27}\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}$
(69)

We now take
$M=\frac{n}{2m^{3/2}}=\left(\frac{27{\color[rgb]{0,0,1}\varepsilon}Wi^{2}}{160}\right)^{1/2}$

$\frac{d{t_{1}}}{d{x}}=m^{1/2}\Big{(}\left(-\frac{M}{t_{1}}+\left(\frac{M^{2}}{t_{1}^{2}}+\frac{1}{27}\right)^{1/2}\right)^{1/3}+\\\
\left(-\frac{M}{t_{1}}-\left(\frac{M^{2}}{t_{1}^{2}}+\frac{1}{27}\right)^{1/2}\right)^{1/3}\Big{)}$
(70)

On substituting $t_{1}=3\sqrt{3}M\tan\theta$, we get the following equation

$\frac{9M}{\sqrt{m}}\frac{\sec^{2}\theta\tan^{1/3}\theta}{(\sec\theta-1)^{1/3}-(\sec\theta+1)^{1/3}}d\theta=dx$
(71)

This can be further simplified to give us :

$\frac{9M}{\sqrt{m}}\frac{\sec^{2}\theta\tan^{1/3}\frac{\theta}{2}}{\tan^{2/3}\frac{\theta}{2}-1}d\theta=dx$
(72)

We then take $v=(\tan\frac{\theta}{2})^{1/3}$

$\frac{54M}{\sqrt{m}}\frac{(v^{6}+1)(v^{3})dv}{(v^{2}-1)(v^{6}-1)^{2}}=dx$
(73)

Using partial fractions, we can get the final form of the integrated equation

$I(v)=x+C_{1}$ (74)

where $I(v)$ is given by :

$I(v)=\frac{54M}{\sqrt{m}}\left[\frac{log|v^{2}-1|}{27}-\frac{log|v^{4}+v^{2}+1|}{54}+\frac{v^{6}-5v^{4}+v^{2}}{18(v^{2}-1)(v^{6}-1)}\right]$
(75)

and $C_{1}$ is an integration constant. Furthermore, on substituting the BC
$h(x=1)=1$ we get an implicit equation for $h(x)$ and $x$ given by :

$I(v)=x+I(v_{1})$ (76)

where $v_{1}=(\tan\frac{\theta_{1}}{2})^{1/3}$ and
$\theta_{1}=\arctan{\frac{1}{3\sqrt{3}M}}$

##  Appendix B  Equations for finite difference

We use a central difference scheme for the spatial derivatives given by :

$\frac{\partial{p}}{\partial{x}}\equiv\frac{p^{n}_{i+1}-2p^{n}_{i}+p^{n}_{i-1}}{(\Delta
x)^{2}}$ (77)

and temporal derivatives using the backward difference scheme :

$\frac{\partial p}{\partial t}\equiv\frac{p^{n}_{i}-p^{n-1}_{i}}{\Delta t}$
(78)

The value of pressure at the $(n-1)^{th}$ step is taken to be the same as the
$k^{th}$ Newton iteration of the $(n-1)^{th}$ step. The discretized nonlinear
equation, whose Jacobian is calculated for the Newton iterations, is given as
follows:

$\frac{p^{n}_{i}-p^{n-1}_{i}}{\Delta
t}-\frac{2}{3\beta}\left(\frac{p^{n}_{i+1}-2p^{n}_{i}+p^{n}_{i-1}}{(\Delta
x)^{2}}\right)-2(1+\beta
p^{n}_{i})^{2}\left(\frac{p^{n}_{i+1}-p^{n}_{i-1}}{2\Delta x}\right)^{2}\\\
-\left(\frac{6{\varepsilon}Wi^{2}}{5\beta}\right)\left(\frac{p^{n}_{i+1}-p^{n}_{i-1}}{2\Delta
x}\right)^{2}\left(\frac{p^{n}_{i+1}-2p^{n}_{i}+p^{n}_{i-1}}{(\Delta
x)^{2}}\right)(1+\beta p^{n}_{i})^{5}-\\\
\left(\frac{2{\varepsilon}Wi^{2}}{5}\right)\left(\frac{p^{n}_{i+1}-p^{n}_{i-1}}{2\Delta
x}\right)^{4}(1+\beta p^{n}_{i})^{4}\equiv
f(p^{n}_{i+1},p^{n}_{i},p^{n}_{i-1})$ (79)

We then proceed to calculate the Jacobian of $f$ and solve the Newton’s method
problem until an error of $1e-7$ is attained based on the absolute value of
the difference between $p^{k}$ and $p^{k+1}$.

## References

  * [1] G. Whitesides, The origins and the future of microfluidics, Nature 442 (2006) 368–373.  doi:10.1038/nature05058.
  * [2] H. Bruus, Theoretical Microfluidics, Oxford Master Series in Physics, OUP Oxford, 2007.
URL https://books.google.com/books?id=wQxREAAAQBAJ

  * [3] J. Berthier, Microfluidics for biotechnology, 2nd Edition, Integrated microsystems series, Artech House, Boston, 2009.
  * [4] A. W. Martinez, S. T. Phillips, G. M. Whitesides, E. Carrilho, Diagnostics for the developing world: Microfluidic paper-based analytical devices, Analytical Chemistry 82 (1) (2010) 3–10, pMID: 20000334.  arXiv:https://doi.org/10.1021/ac9013989, doi:10.1021/ac9013989.
URL https://doi.org/10.1021/ac9013989

  * [5] O. Strohmeier, M. Keller, F. Schwemmer, S. Zehnle, D. Mark, F. von Stetten, R. Zengerle, N. Paust, Centrifugal microfluidic platforms: advanced unit operations and applications, Chemical Society reviews 44 (17) (2015) 6187–6229.  doi:https://doi.org/10.1039/C4CS00371C.
  * [6] N.-T. Nguyen, S. Wereley, S. A. M. Shaegh, Fundamentals and applications of microfluidics, 3rd Edition, Artech house Integrated Microsystems Series, Artech House, Norwood, Massachusetts, 2019.
  * [7] J.-J. Shu, J. Bin Melvin Teo, W. Kong Chan, Fluid velocity slip and temperature jump at a solid surface, Applied Mechanics Reviews 69 (2).
  * [8] S. Chakraborty, Electroosmotically driven capillary transport of typical non-newtonian biofluids in rectangular microchannels, Analytica Chimica Acta 605 (2) (2007) 175–184.  doi:https://doi.org/10.1016/j.aca.2007.10.049.
URL https://www.sciencedirect.com/science/article/pii/S0003267007017850

  * [9] S. Zeng, C.-H. Chen, J. C. Mikkelsen, J. G. Santiago, Fabrication and characterization of electroosmotic micropumps, Sensors and Actuators B: Chemical 79 (2) (2001) 107–114.  doi:https://doi.org/10.1016/S0925-4005(01)00855-3.
URL https://www.sciencedirect.com/science/article/pii/S0925400501008553

  * [10] V. Anand, I. C. Christov, Transient compressible flow in a compliant viscoelastic tube, Physics of Fluids 32 (11) (2020) 112014.  doi:https://doi.org/10.1063/5.0022406.
  * [11] A. Gat, I. Frankel, D. Weihs, Gas flows through constricted shallow micro-channels, Journal of Fluid Mechanics 602 (2008) 427–442.  doi:10.1017/S0022112008001055.
  * [12] I. C. Christov, Soft hydraulics: from newtonian to complex fluid flows through compliant conduits, Journal of Physics:Condensed Matter 34 (6).  doi:https://doi.org/10.1088/1361-648X/ac327d.
  * [13] N. L. Jeon, D. T. Chiu, C. J. Wargo, H. Wu, I. S. Choi, J. R. Anderson, G. M. Whitesides, Design and fabrication of integrated passive valves and pumps for flexible polymer 3-dimensional microfluidic systems, Biomedical Microdevices 4 (2002) 117–121, 775.  doi:https://doi.org/10.1023/A:1014683114796.
  * [14] M. A. Unger, H.-P. Chou, T. Thorsen, A. Scherer, S. R. Quake, Monolithic microfabricated valves and pumps by multilayer soft lithography, Science 288 (5463) (2000) 113–116.  doi:https://doi.org/10.1126/science.288.5463.113.
  * [15] J. C. Lötters, W. Olthuis, P. H. Veltink, P. Bergveld, The mechanical properties of the rubber elastic polymer polydimethylsiloxane for sensor applications, Journal of Micromechanics and Microengineering 7 (3) (1997) 145–147.  doi:10.1088/0960-1317/7/3/017.
URL https://doi.org/10.1088/0960-1317/7/3/017

  * [16] S. Elbaz, H. Jacob, A. Gat, Transient gas flow in elastic microchannels, Journal of Fluid Mechanics 846 (2018) 460–481.  doi:10.1017/jfm.2018.287.
  * [17] C. Duprat, H. Stone (Eds.), Fluid-Structure Interactions in Low-Reynolds-Number Flows, Soft Matter Series, The Royal Society of Chemistry, 2016.  doi:10.1039/9781782628491.
URL http://dx.doi.org/10.1039/9781782628491

  * [18] T. Gervais, J. El-Ali, A. Günther, K. F. Jensen, Flow-induced deformation of shallow microfluidic channels, Lab Chip 6 (2006) 500–507.  doi:10.1039/B513524A.
URL http://dx.doi.org/10.1039/B513524A

  * [19] M. K. Raj, S. Dasgupta, S. Chakraborty, Hydrodynamics in deformable microchannels, Microfluidics and Nanofluidics 21 (2017) 1–12.  doi:https://doi.org/10.1007/s10404-017-1908-5.
  * [20] I. C. Christov, V. Cognet, T. C. Shidhore, H. A. Stone, Flow rate–pressure drop relation for deformable shallow microfluidic channels, Journal of Fluid Mechanics 841 (2018) 267–286.  doi:10.1017/jfm.2018.30.
  * [21] B. S. Hardy, K. Uechi, J. Zhen, H. Pirouz Kavehpour, The deformation of flexible pdms microchannels under a pressure driven flow, Lab Chip 9 (2009) 935–938.  doi:10.1039/B813061B.
URL http://dx.doi.org/10.1039/B813061B

  * [22] E. Seker, D. C. Leslie, H. Haj-Hariri, J. P. Landers, M. Utz, M. R. Begley, Nonlinear pressure-flow relationships for passive microfluidic valves, Lab Chip 9 (2009) 2691–2697.  doi:10.1039/B903960K.
URL http://dx.doi.org/10.1039/B903960K

  * [23] P. Cheung, K. Toda-Peters, A. Q. Shen, In situ pressure measurement within deformable rectangular polydimethylsiloxane microfluidic devices, Biomicrofluidics 6 (2) (2012) 026501.  arXiv:https://doi.org/10.1063/1.4720394, doi:10.1063/1.4720394.
URL https://doi.org/10.1063/1.4720394

  * [24] A. Raj, A. K. Sen, Flow-induced deformation of compliant microchannels and its effect on pressure–flow characteristics, Microfluidics and Nanofluidics 20 (2) (2016) 1–13.  doi:10.1007/s10404-016-1702-9.
  * [25] A. Martínez-Calvo, A. Sevilla, G. G. Peng, H. A. Stone, Start-up flow in shallow deformable microchannels, Journal of Fluid Mechanics , Vol. 885(2020)doi:https://doi.org/10.1017/jfm.2019.994.
  * [26] A. Hosoi, L. Mahadevan, Peeling, healing, and bursting in a lubricated elastic sheet, Physical review letters 93 (13) (2004) 137802.
  * [27] T. Ball, J. Neufeld, Static and dynamic fluid-driven fracturing of adhered elastica, Physical Review Fluids 3.  doi:10.1103/PhysRevFluids.3.074101.
  * [28] V. Ciriello, A. Lenci, S. Longo, V. Di Federico, Relaxation-induced flow in a smooth fracture for ellis rheology, Advances in Water Resources 152 (2021) 103914.  doi:https://doi.org/10.1016/j.advwatres.2021.103914.
URL https://www.sciencedirect.com/science/article/pii/S0309170821000695

  * [29] Y. Matia, A. D. Gat, Dynamics of elastic beams with embedded fluid-filled parallel-channel networks, Soft robotics 2 (1) (2015) 42–47.
  * [30] J. B. Grotberg, Pulmonary flow and transport phenomena, Annual review of fluid mechanics 26 (1) (1994) 529–571.
  * [31] J. B. Grotberg, O. E. Jensen, Biofluid mechanics in flexible tubes, Annual Review of Fluid Mechanics 36 (1) (2004) 121–147.  doi:10.1146/annurev.fluid.36.050802.121918.
  * [32] M. HEIL, Stokes flow in collapsible tubes: computation and experiment, Journal of Fluid Mechanics 353 (1997) 285–312.  doi:10.1017/S0022112097007490.
  * [33] J. R. Womersley, Oscillatory flow in arteries: the constrained elastic tube as a model of arterial flow and pulse transmission, Physics in Medicine and Biology 2 (2) (1957) 178–187.  doi:10.1088/0031-9155/2/2/305.
URL https://doi.org/10.1088/0031-9155/2/2/305

  * [34] K. A. S. Boster, T. C. Shidhore, A. A. Cohen-Gadol, I. C. Christov, V. L. Rayz, Challenges in modeling hemodynamics in cerebral aneurysms related to arteriovenous malformations, Cardiovascular Engineering and Technology (2022) 1–12doi:https://doi.org/10.1007/s13239-022-00609-3.
  * [35] R. Bird, R. Bird, R. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Volume 1: Fluid Mechanics, Dynamics of Polymeric Liquids, Wiley, 1987\.
URL https://books.google.com/books?id=posvAQAAIAAJ

  * [36] R. Chhabra, J. Richardson, Non-Newtonian Flow and Applied Rheology: Engineering Applications, Elsevier Science, 2011.
URL https://books.google.com/books?id=_6nnoh9PtF0C

  * [37] V. Anand, J. David Jr, I. C. Christov, Non-newtonian fluid–structure interactions: Static response of a microchannel due to internal flow of a power-law fluid, Journal of Non-Newtonian Fluid Mechanics 264 (2019) 62–72.  doi:https://doi.org/10.1016/j.jnnfm.2018.12.008.
  * [38] L. Ramos Arzola, O. Bautista, Fluid structure-interaction in a deformable microchannel conveying a viscoelastic fluid, Journal of Non-Newtonian Fluid Mechanics 296.  doi:10.1016/j.jnnfm.2021.104634.
  * [39] V. Kumaran, G. H. Fredrickson, P. Pincus, Flow induced instability of the interface between a fluid and a gel at low reynolds number, J. Phys. II France 4 (6) (1994) 893–911.  doi:10.1051/jp2:1994173.
URL https://doi.org/10.1051/jp2:1994173

  * [40] V. Kumaran, Stability of the viscous flow of a fluid through a flexible tube, Journal of Fluid Mechanics 294 (1995) 259–281.  doi:10.1017/S0022112095002886.
  * [41] S. A. Roberts, S. Kumar, Stability of creeping couette flow of a power-law fluid past a deformable solid, Journal of Non-Newtonian Fluid Mechanics 139 (1) (2006) 93–102.  doi:https://doi.org/10.1016/j.jnnfm.2006.07.006.
URL https://www.sciencedirect.com/science/article/pii/S0377025706001704

  * [42] V. Shankar, S. Kumar, Instability of viscoelastic plane couette flow past a deformable wall, Journal of Non-Newtonian Fluid Mechanics 116 (2) (2004) 371–393.  doi:https://doi.org/10.1016/j.jnnfm.2003.10.003.
URL https://www.sciencedirect.com/science/article/pii/S0377025703002490

  * [43] X. Wang, I. C. Christov, Reduced modeling and global instability of finite-reynolds-number flow in compliant rectangular channelsdoi:10.48550/ARXIV.2202.11704.
URL https://arxiv.org/abs/2202.11704

  * [44] V. Anand, S. C. Muchandimath, I. C. Christov, Hydrodynamic Bulge Testing: Materials Characterization Without Measuring Deformation, Journal of Applied Mechanics 87 (5), 051012.  arXiv:https://asmedigitalcollection.asme.org/appliedmechanics/article-pdf/87/5/051012/6512238/jam\\_87\\_5\\_051012.pdf, doi:10.1115/1.4046297.
URL https://doi.org/10.1115/1.4046297

  * [45] D. Coyle, Forward roll coating with deformable rolls: A simple one-dimensional elastohydrodynamic model, Chemical Engineering Science 43 (10) (1988) 2673–2684.  doi:https://doi.org/10.1016/0009-2509(88)80011-3.
URL https://www.sciencedirect.com/science/article/pii/0009250988800113

  * [46] J. M. Skotheim, L. Mahadevan, Soft lubrication: The elastohydrodynamics of nonconforming and conforming contacts, Physics of Fluids 17 (9) (2005) 092101.  arXiv:https://doi.org/10.1063/1.1985467, doi:10.1063/1.1985467.
URL https://doi.org/10.1063/1.1985467

  * [47] J. M. Skotheim, L. Mahadevan, Soft lubrication, Phys. Rev. Lett. 92 (2004) 245509.  doi:10.1103/PhysRevLett.92.245509.
URL https://link.aps.org/doi/10.1103/PhysRevLett.92.245509

  * [48] X. Yin, S. Kumar, Lubrication flow between a cavity and a flexible wall, Physics of Fluids 17 (6) (2005) 063101\.  arXiv:https://doi.org/10.1063/1.1914819, doi:10.1063/1.1914819.
URL https://doi.org/10.1063/1.1914819

  * [49] X. Yin, S. Kumar, Two-dimensional simulations of flow near a cavity and a flexible solid boundary, Physics of Fluids 18 (6) (2006) 063103.  arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.2204061, doi:10.1063/1.2204061.
URL https://aip.scitation.org/doi/abs/10.1063/1.2204061

  * [50] X. Wang, I. Christov, Theory of the flow-induced deformation of shallow compliant microchannels with thick walls, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475 (2019) 20190513.  doi:10.1098/rspa.2019.0513.
  * [51] N. P. Thien, R. I. Tanner, A new constitutive equation derived from network theory, Journal of Non-Newtonian Fluid Mechanics 2 (4) (1977) 353–365.  doi:https://doi.org/10.1016/0377-0257(77)80021-9.
  * [52] L. Ferrás, M. Morgado, M. Rebelo, G. H. McKinley, A. Afonso, A generalised phan–thien—tanner model, Journal of Non-Newtonian Fluid Mechanics 269 (2019) 88–99.  doi:https://doi.org/10.1016/j.jnnfm.2019.06.001.
URL https://www.sciencedirect.com/science/article/pii/S0377025718301502

  * [53] A. Mehboudi, J. Yeom, Experimental and theoretical investigation of a low-reynolds-number flow through deformable shallow microchannels with ultra-low height-to-width aspect ratios, Microfluidics and Nanofluidics 23 (5) (2019) 1–14.  doi:https://doi.org/10.1007/s10404-019-2235-9.
  * [54] M. Kumar, J. S. Guasto, A. M. Ardekani, Transport of complex and active fluids in porous media, Journal of Rheology 66 (2) (2022) 375–397.  doi:10.1122/8.0000389.
  * [55] R. G. Larson, E. S. G. Shaqfeh, S. J. Muller, A purely elastic instability in taylor–couette flow, Journal of Fluid Mechanics 218 (1990) 573–600.  doi:10.1017/S0022112090001124.
  * [56] L. Zhong, M. Oostrom, M. Truex, V. Vermeul, J. Szecsody, Rheological behavior of xanthan gum solution related to shear thinning fluid delivery for subsurface remediation, Journal of Hazardous Materials 244-245 (2013) 160–170.  doi:https://doi.org/10.1016/j.jhazmat.2012.11.028.
  * [57] S. Ansari, M. A. I. Rashid, P. R. Waghmare, D. S. Nobes, Measurement of the flow behavior index of newtonian and shear-thinning fluids via analysis of the flow velocity characteristics in a mini-channel, SN Applied Sciences 2 (11) (2020) 1787.
  * [58] R. Dangla, F. Gallaire, C. N. Baroud, Microchannel deformations due to solvent-induced pdms swelling, Lab Chip 10 (2010) 2972–2978.  doi:10.1039/C003504A.
URL http://dx.doi.org/10.1039/C003504A

  * [59] G. Romeo, G. D’Avino, F. Greco, P. A. Netti, P. L. Maffettone, Viscoelastic flow-focusing in microchannels: scaling properties of the particle radial distributions, Lab Chip 13 (2013) 2802–2807.  doi:10.1039/C3LC50257K.
URL http://dx.doi.org/10.1039/C3LC50257K

  * [60] F. Del Giudice, S. J. Haward, A. Q. Shen, Relaxation time of dilute polymer solutions: A microfluidic approach, Journal of Rheology 61 (2) (2017) 327–337.  doi:10.1122/1.4975933.
  * [61] K. W. Ebagninin, A. Benchabane, K. Bekkour, Rheological characterization of poly(ethylene oxide) solutions of different molecular weights, Journal of Colloid and Interface Science 336 (1) (2009) 360–367.  doi:https://doi.org/10.1016/j.jcis.2009.03.014.
URL https://www.sciencedirect.com/science/article/pii/S002197970900304X

  * [62] S. Longo, V. Di Federico, L. Chiapponi, A dipole solution for power-law gravity currents in porous formations, Journal of Fluid Mechanics 778 (2015) 534–551.  doi:10.1017/jfm.2015.405.
  * [63] P. J. Oliviera, F. T. Pinho, Analytical solution for fully developed channel and pipe flow of phan-thien–tanner fluids, Journal of Fluid Mechanics 387 (1999) 271–280.  doi:10.1017/S002211209900453X.
  * [64] H. Ahmed, L. Biancofiore, A new approach for modeling viscoelastic thin film lubrication, Journal of Non-Newtonian Fluid Mechanics 292 (2021) 104524.  doi:https://doi.org/10.1016/j.jnnfm.2021.104524.
URL https://www.sciencedirect.com/science/article/pii/S037702572100046X

  * [65] E. Boyko, H. A. Stone, Reciprocal theorem for calculating the flow rate–pressure drop relation for complex fluids in narrow geometries, Phys. Rev. Fluids 6 (2021) L081301.  doi:10.1103/PhysRevFluids.6.L081301.
URL https://link.aps.org/doi/10.1103/PhysRevFluids.6.L081301

  * [66] A. Martínez-Calvo, A. Sevilla, G. G. Peng, H. A. Stone, Start-up flow in shallow deformable microchannels, J. Fluid Mech. 885 (2020) A25.  doi:10.1017/jfm.2019.994.
  * [67] A. A. Ghodgaonkar, I. C. Christov, Solving Nonlinear Parabolic Equations by a Strongly Implicit Finite Difference Scheme, Springer International Publishing, Cham, 2019, pp. 305–342.  doi:10.1007/978-3-030-29951-4_14.
URL https://doi.org/10.1007/978-3-030-29951-4_14

  * [68] H. Langtangen, S. Linge, Finite Difference Computing with PDEs, Vol. 16, 2017.  doi:10.1007/978-3-319-55456-3.
  * [69] E. Boyko, H. A. Stone, Pressure-driven flow of the viscoelastic oldroyd-b fluid in narrow non-uniform geometries: analytical results and comparison with simulations, Journal of Fluid Mechanics 936 (2022) A23.  doi:10.1017/jfm.2022.67.
  * [70] E. Boyko, M. Bercovici, A. D. Gat, Viscous-elastic dynamics of power-law fluids within an elastic cylinder, Phys. Rev. Fluids 2 (2017) 073301.  doi:10.1103/PhysRevFluids.2.073301.
URL https://link.aps.org/doi/10.1103/PhysRevFluids.2.073301

  * [71] M. Yi, Z. Shen, Fluid dynamics: an emerging route for the scalable production of graphene in the last five years, RSC Adv. 6 (2016) 72525–72536.  doi:10.1039/C6RA15269D.
URL http://dx.doi.org/10.1039/C6RA15269D

  * [72] G. G. Peng, J. R. Lister, Viscous flow under an elastic sheet, Journal of Fluid Mechanics 905 (2020) A30.  doi:10.1017/jfm.2020.745.
  * [73] A. Y. Fu, H.-P. Chou, C. Spence, F. H. Arnold, S. R. Quake, An integrated microfabricated cell sorter, Analytical Chemistry 74 (11) (2002) 2451–2457, pMID: 12069222.  doi:10.1021/ac0255330.