# Low Reynolds number pulsatile flow of a viscoelastic fluid through a
channel: Effects of fluid rheology and pulsation parameters

Subhasisa Rath Corresponding Author<EMAIL_ADDRESS>ORCID:0000-0002-4202-7434 [Address all correspondence related to ASME style
format and figures to this author]  School of Energy Science & Engineering
Indian Institute of Technology Kharagpur
Kharagpur 721302, West Bengal, India  Bimalendu Mahapatra Corresponding
Author<EMAIL_ADDRESS>ORCID:0000-0003-0845-9775 [Address
all correspondence for other issues to this author]  Department of Mechanical
Engineering
Indian Institute of Technology Kharagpur
Kharagpur 721302, West Bengal, India

###### Abstract

As the first endeavour, we have analyzed the pulsatile flow of Oldroyd-B
viscoelastic fluid where the combined effects of fluid elasticity and
pulsation parameters on the flow characteristics are numerically studied at a
low Reynolds number. Computations are performed using a finite-volume based
open-source solver OpenFOAM® by appending the log-conformation tensor approach
to stabilize the numerical solution at high Deborah number. Significant flow
velocity enhancement is achieved by increasing the viscoelastic behaviour of
the fluid. High-velocity gradient zones and high polymeric stress regions are
observed near the channel wall. The magnitude of axial velocity attenuates
with increasing pulsation amplitude or pulsation frequency, and the extent of
this attenuation is highly dependent on the Deborah number or the retardation
ratio. This work finds application in the transport of polymeric solutions,
extrusion, and injection moulding of polymer melts in several process
industries.

##  1 Introduction

Pulsatile flow of viscoelastic fluids plays an indispensable role in various
process industries as well as in biomedical and biochemical applications [1,
2, 3, 4, 5]. It promotes inherent advantages in mixing and separations of
chemical species, transport of polymeric solutions or suspensions, injection
moulding of polymer melts, blood flow through arteries and veins, and drug
delivery [6, 7, 8, 9, 10, 11, 12]. Viscoelastic fluids are the class of non-
Newtonian fluid where some polymers are blend with a solvent. Thus,
viscoelastic fluids exhibit both viscous and elastic behaviours, i.e., elastic
recovery after deformation and stress relaxation [13]. Owing to the time-
dependent behaviour or relaxation phenomena and memory effects, the underlying
physics of viscoelastic fluids is quite complex. The interplay of complex
rheology of the viscoelastic fluid and the pulsatile nature of the flow
further accentuate the complexity of the hydrodynamic problem. Hence, proper
understanding and characterization of viscoelastic fluids are vital for their
applicability in real-life practical scenarios.

Over the last few decades, a noticeable development in research have been seen
for non-Newtonian viscoelastic fluids. Several linear, non-linear, and quasi-
linear viscoelastic fluid models have been developed to apprehend the
intricate flow physics associated with the elastic behaviour of the complex
fluids. For instance, the Giesekus model[14], Maxwell model[15], Jeffreys
model[16], Oldroyd-B model[17], Phan–Thien–Tanner (PTT) model[18], Finitely
Extensible Nonlinear Elastic (FENE) model[19], and Pom–Pom model[20] are few
well-known viscoelastic fluid models in literature. In numerical computations,
an elementary problem associated with the viscoelastic fluid model is the
instability or convergence issue at high Weisenberg number ($Wi$), or Deborah
number ($De$)[21]. Hence, a verity of approaches have been proposed to
overcome such problem and stabilize the numerical solution such as: log-
conformation tensor approach [22, 23, 24], both-sides-diffusion technique[25,
26, 27], elastic-viscous stress splitting methods[28, 29], and solvent-polymer
decomposition method[13, 30].

Aiming to better understand the complex flow dynamics of pulsatile flow of
biological fluids, several studies were conducted using Newtonian and
viscoelastic fluid models [31, 32, 33, 34, 1]. Cable and Boger [35]
experimentally investigated the unstable abrupt entry flow of viscoelastic
fluids in a tube using a periodic spiralling frequency method. The periodic
entry flow disturbances, criteria for the onset of periodic disturbances, and
random entry flow disturbances were elaborately explained in their study. The
oscillating and pulsating flow of viscoelastic fluids were experimentally
studied by Khabakhpasheva $et~{}al.$ [36] using a sinusoidal time-varying
pressure gradient. Axial velocity profiles as a function of the pulsation
phase are plotted and compared with the pressure drop profiles. An enhancement
in flow rate has been observed compared to a steady flow, and the flow
enhancement with amplitude of pulsation for different ranges of $Re$, $De$,
and $Wi$ were discussed in their study. Letelier $et~{}al.$ [37] and Siginer
and Letelier [38] analytically solved the pulsatile flow of viscoelastic
fluids in different cross-section of straight tubes to delineate the
longitudinal field and secondary flows, respectively. Duarte $et~{}al.$ [30]
numerically and analytically investigated both start-up and pulsating test
case problem for unsteady viscoelastic flows. A periodic pressure gradient has
been superimposed on the constant Poiseuille flow to demonstrate the pulsatile
flow problem using various viscoelastic fluid models. Rachid and Ouazzani [39]
analytically investigated the effect of pulsatile flow on low $Re$ peristaltic
transport of a viscoelastic fluid in a cylindrical tube using the Maxwell
model.

Pulsatile flow of viscoelastic fluids plays a vital role in biochemical and
biomedical applications to mimic the blood flow in arteries, veins, and
vessels. Javadzadegan [40] numerically studied the pulsatile flow of both
viscous and viscoelastic fluids in a constricted stenosed artery using
oscillatory pressure gradient with Oldroyd-B and Cross models. Flow rate,
velocity profiles, and wall shear-stress have been clearly elaborated in their
study. Bakhti $et~{}al.$ [41] numerically investigated the pulsatile blood
flow through a constricted tapered artery using the generalized Oldroyd-B
model with variable-order fractional derivative technique. Moyers $et~{}al.$
[42] analytically solved the high-frequency oscillatory blood flow in a tube
using the non-homogeneous hemorheological model. The amplitude relationship
between the pressure gradient and volume flow rate has been investigated in
their study. Sajid $et~{}al.$ [43] numerically demonstrated the effect of body
acceleration on the pulsatile flow of blood in a vessel using the Oldroyd-B
viscoelastic model. To elucidate the effects of stress relaxation and
retardation in unsteady pulsatile flow was the primary aim of their work.

The previous studies on the pulsatile flow of viscoelastic fluids have
revealed many physical insights into the flow physics where a time-varying
pressure gradient was imposed to generate the pulsating flow. The time-varying
pressure gradient gives rise to a time-varying velocity pulsation in the flow,
which may or may not be in the same phase depending upon the flow conditions
[36]. As the first endeavour in the field of viscoelastic fluids, a sinusoidal
axial velocity pulsation has been superimposed on the mean flow to trigger a
pulsatile flow in the channel, albeit research has been done in Newtonian and
power-law viscoplastic fluids [44, 45]. In real life practical experiments,
one can convert the velocity to flow rate and use the time-varying flow rate
to mimic the pulsatile flow in a duct or channel. In the present study, a
pulsatile flow of viscoelastic fluid in a parallel plate channel has been
studied numerically to elucidate the combined effects of pulsation parameters
and complex rheology of the viscoelastic fluid on flow characteristics. This
work finds application in the transport of polymeric solutions, extrusion and
injection moulding of the polymer melt in several process industries, and
transport of biofluids in microfluidic devices.

##  2 Problem description

In the present study, we have considered a parallel-plate channel of height
$2H$ and length $60H$ to simulate the pulsatile flow of viscoelastic fluid, as
shown in Fig. 1. The origin of the Cartesian coordinate system is located
exactly at the centre of the inlet such that the top and bottom walls are
placed at $y=\pm H$, respectively. A sinusoidal pulse, $a\sin(\omega t)$, is
superimposed on the non-zero mean or base flow $(U_{0})$ to generate a time-
varying axial pulsatile velocity $U_{x}$. Where $a$ is the amplitude of
pulsation, $\omega$ is the angular frequency of pulsation, and $t$ is the
time. The viscoelastic fluid described by the Oldroyd-B constitutive equations
flows inside the channel due to the imposition of velocity pulsation. We have
illustrated the governing equations, non-dimensionalization of equations, and
relevant boundary conditions for the present flow configuration in the
forthcoming sections.

###  2.1 Governing equations

The dimensional governing equations for conservation of mass and linear-
momentum in the absence of any external body force for an incompressible,
isothermal, laminar, and two-dimensional flow are given by:

$\nabla\cdot\mathbf{U}=0$ (1) $\rho\frac{\partial\mathbf{U}}{\partial
t}+\rho\nabla\cdot(\mathbf{U}\cdot\mathbf{U})=-\nabla
p+\nabla\cdot\boldsymbol{\tau}$ (2)

Where, $\rho$ denotes the mass density, $\mathbf{U}$ is the velocity vector
having its components $U_{x}$ and $U_{y}$ in axial and transverse directions,
respectively. The symbol ‘$p$’ represents the fluid pressure and
$\boldsymbol{\tau}$ is the stress tensor of the viscoelastic fluid which shows
a complex relationship with the strain rate [13]. In this study, we implement
a quasi-linear viscoelastic fluid more precisely Oldroyd-B model [17] to
capture the rheological behavior of the fluid. The extra stress tensor
$\boldsymbol{\tau}$ for the Oldroyd-B constitutive model is described by:

$\boldsymbol{\tau}+\lambda\overset{\nabla}{\boldsymbol{\tau}}=2\eta\Big{(}\mathbf{D}+\lambda(1-\kappa)\overset{\nabla}{\mathbf{D}}\Big{)}$
(3)

In Eq. 3, $\lambda$ is the fluid relaxation time [46] and $\eta$ is the total
viscosity of the fluid, where $\eta=\eta_{s}+\eta_{p}$. Here, $\eta_{s}$ and
$\eta_{p}$ are the solvent and polymeric viscosity contributions at zero
shear-rate, respectively. The symbol ‘$\kappa$’ is defined as the ratio of
polymeric viscosity to the total viscosity and the fluid retardation time is
of the form $\lambda_{r}=\lambda(1-\kappa)$. The rate of deformation tensor is
denoted as $\mathbf{D}$, where
$\mathbf{D}$=$\frac{1}{2}[\nabla\mathbf{U}+\nabla\mathbf{U}^{\mathrm{T}}]$.
Finally, $\overset{\nabla}{\boldsymbol{\tau}}$ and
$\overset{\nabla}{\mathbf{D}}$ are the upper-convected derivative (UCD) terms
of $\boldsymbol{\tau}$ and $\mathbf{D}$. The UCD of any second-order tensor
$\mathbf{M}$ is defined as

$\overset{\nabla}{\mathbf{M}}=\frac{\partial\mathbf{M}}{\partial
t}+\mathbf{U}\cdot\nabla\mathbf{M}-\mathbf{M}\cdot\nabla\mathbf{U}-\nabla\mathbf{U}^{\mathrm{T}}\cdot\mathbf{M}$
(4)

###  2.2 Non-dimensionalization

The dimensional governing equations employed to solve the present problem are
non-dimensionalized using proper scaling of variables to obtain a
dimensionless set of equations and reduced number of variables to define the
flow physics. Here, the superscript ‘$\sim$’ denotes the non-dimensional
parameters which are scaled using the characteristic values as defined below.

$\displaystyle\widetilde{t}=\frac{tU_{0}}{H}~{},~{}~{}\widetilde{x}=\frac{x}{H}~{},~{}~{}\widetilde{y}=\frac{y}{H}~{},~{}~{}\widetilde{U}_{x}=\frac{u_{x}}{U_{0}}~{},~{}~{}\widetilde{U}_{y}=\frac{u_{y}}{U_{0}}~{},{}{}$
(5) $\displaystyle\widetilde{p}=\frac{pH}{\eta
U_{0}}~{},~{}~{}\widetilde{\tau}_{xx}=\frac{\tau_{xx}H}{\eta
U_{0}}~{},~{}~{}\widetilde{\tau}_{xy}=\frac{\tau_{xy}H}{\eta
U_{0}}~{},~{}~{}\widetilde{\tau}_{yy}=\frac{\tau_{yy}H}{\eta U_{0}}$

Here, we have considered the advection time to be the characteristic time
scale, the half-height of the channel as the characteristic length scale, the
mean flow velocity as the characteristic velocity scale and viscous stress to
scale the fluid pressure and stress components.

###  2.3 Flow governing parameters

It is worthwhile to define all the governing dimensionless parameters and
there range of values associated with the present work. We have defined all
the dimensionless parameters and their range of values in Table 1. Here the
Reynolds number is defined as the ratio of inertial forces to the viscous
forces based on the mean flow velocity. The ratio of the characteristic time
scale, which is the relaxation time ($\lambda$) of the fluid in this study, to
the time scale of deformation is termed as the Deborah number [47]. The ratio
of the retardation time ($\lambda_{r}$) to the relaxation time ($\lambda$) of
the viscoelastic fluid, which is equal to the ratio of solvent viscosity to
the total viscosity of the fluid, is defined as the retardation ratio or
viscosity ratio ($\beta$) [30, 47]. The dimensionless pulsatile flow frequency
scaled with the viscous effects is known as the Womersley number ($Wo$)[48].
The dimensionless pulsation amplitude ($A$) is defined as the amplitude of
pulsatile velocity ($a$) scaled with the flow velocity ($U_{0}$). After the
scaling analysis, the dimensional governing equations, Eqs. 1,2, and 3, are
transformed to the their dimensionless form as Eqs. 6,7, and 8, respectively.

$\widetilde{\nabla}\cdot\widetilde{\mathbf{U}}=0$ (6)
$\frac{\partial\widetilde{\mathbf{U}}}{\partial\widetilde{t}}+\widetilde{\nabla}\cdot(\widetilde{\mathbf{U}}\widetilde{\mathbf{U}})=-\widetilde{\nabla}\widetilde{p}+\frac{1}{Re}\widetilde{\nabla}\widetilde{\boldsymbol{\tau}}$
(7)
$\widetilde{\bm{\tau}}+De~{}\overset{\nabla}{\widetilde{\boldsymbol{\tau}}}=2\Bigg{(}\widetilde{\mathbf{D}}+\beta
De~{}\overset{\nabla}{\widetilde{\mathbf{D}}}\Bigg{)}$ (8)

##  3 Computational details

In the present study, we have used an open-source toolbox ‘$RheoTool$’ [49]
for the numerical computations of viscoelastic fluid flows in the presence of
inlet velocity pulsation. The governing conservation equations, Eqs. 1 and 2,
along with the constitutive relation, Eq. 3, are solved by fully implicit
finite-volume method (FVM) based on a time-marching pressure-correction
algorithm [50, 51] in OpenFOAM® framework. We have adopted the Semi-Implicit
Method for Pressure-Linked Equations-Consistent (SIMPLEC) algorithm [52, 53]
for pressure-velocity coupling. To increase the stability of the solver, a new
momentum-stress coupling algorithm (log-conformation tensor approach) is
implemented in the flow solver [52]. The detailed explanation for the log-
conformation tensor approach is presented in the Suplementary material.

###  3.1 Boundary conditions:

In order to solve the aforementioned dimensionless governing equation
numerically, we impose a time varying pulsating axial flow velocity variation
at the inlet boundary and both pressure and polymeric extra-stress components
are set to zero-gradient. At the outlet, pressure is set to atmospheric
whereas, all other variables are assigned as zero-gradient. No-slip and no-
penetration stationary wall boundary condition is imposed to both the top and
bottom walls. At these walls, pressure is set to zero normal gradient whereas,
polymeric extra-stress components are linearly extrapolated.

###  3.2 Grid independence analysis

In order to ensure that the numerical results are independent of the size of
the user-generated grids, a grid sensitivity test has been demonstrated by
generating three different mesh structures, namely M1, M2, and M3, as given in
Table 2. In this study, the $blockMesh$ utility tool is employed to generate
mesh for the simulations, which is available in ‘OpenFOAM®. For distributing
mapped cells in the domain, the computational domain has been split into two
zones, as shown in Fig. 2(a). The subplot (b) of Fig. 2 shows the non-uniform
grid used to fully resolve the flow field in the vicinity of the channel wall.
In Table 2, $n_{x}$ and $n_{y}$ represent the number of cells in $x$ and $y$
directions, respectively in each zone and $N_{c}$ represents the total number
of cells in the computational domain. Figure 2(c) shows the axial velocity
profiles along the transverse direction at $\widetilde{x}=40$ for three
different meshes at $De=10$. From Table 2 and Fig. 2(c), the relative
variation in peak axial velocity is seen to be less than $1\%$ between the
meshes M2 and M3. Hence, the mesh M2 is seemed as grid-independent and taken
for the present computations. To cross verify the grid sensitivity test, the
“Grid Convergence Index” (GCI) method [54, 55, 56] is also employed in this
study. The details about the GCI algorithm and its results are provided in the
Supplementary material.

###  3.3 Model validation

In order to check the accuracy and reliability of the present numerical model,
the present numerical result is validated with the existing analytical
solution of Na and Yoo [57] for a uniform flow of viscoelastic fluid in a
parallel plate channel. Figure 3 represents the variation of axial velocity,
axial normal stress, and shear stress profiles along the transverse direction
at $\widetilde{x}=40$ for $Re=1$, $\beta=0.01$, and $De=0.1$. The present
numerical result shows a close agreement with Na and Yoo [57] within $1\%$
deviation. Hence, it can be believed that the numerical model implemented here
would be able to predict much accurate solution for the present study.

##  4 Results and discussion

In this study, the time-averaged velocity and polymeric stress component
profiles along the transverse and axial directions are presented to elucidate
the effects of $De$, $\beta$, $Wo$, and $A$ on the flow characteristics
associated with the interplay of pulsatile flow and viscoelastic nature of the
fluid. However, before going to discuss the results, it is worthwhile to
understand the nature of the pulsation parameters, i.e., $Wo$ and $A$. Thus,
in Fig. 4(a,b), the pulsation of axial velocity is denoted at a probe location
of $\widetilde{x}=40,\widetilde{y}=+0.5$ in the flow domain to explain the
effect of $Wo$ and $A$, respectively for $\beta=0.1$ and $De=1$. In Fig. 4(a),
pulsation shows for different frequencies having same amplitude of $A=0.1$ for
four seconds which denotes two complete cycles of $Wo=3.545$ or $\omega=4\pi$
, four complete cycles of $Wo=2.507$ or $\omega=2\pi$, and eight complete
cycles of $Wo=1.772$ or $\omega=\pi$. Similarly, in Fig. 4(b), pulsation shows
for different amplitudes having the same frequency of $Wo=2.507$ or
$\omega=2\pi$ for two seconds which denotes two complete cycles. From Fig.4,
it is clear that with an increase in $Wo$, the frequency of the axial velocity
pulsation increases and a higher $A$ enhances the magnitude of
$\widetilde{U}_{x}$. In the forthcoming sections, we will discuss the effects
of changing fluid properties and pulsation parameters on velocity and stress
distributions in the channel.

###  4.1 Effect of fluid elasticity

Unlike the Newtonian fluid where the stress is a linear function of the rate
of strain, the viscoleastic fluid exhibit a non-linear time dependent relation
between the fluid stress and the rate of strain. The viscoelastic fluid model
constitute of two-time scales, namely relaxation time $\lambda$ and
retardation time $\lambda_{r}$, which control the elastic behaviour of the
fluid. In the presence of flow pulsation, it is expected that these fluid
properties may significantly alter the flow dynamics. In Fig. 5(a,b), we have
highlighted the axial velocity variation for varying $De$ (non-dimensional
representation of relaxation time $\lambda$) at $\widetilde{x}=40$ for two
different $\beta$ values. From Fig. 5(a), we observe that for a constant
$\beta$ of 0.1, the centerline velocity increases with an increase in $De$ due
to the increased elasticity of the fluid. Further, keeping the $De=50$, the
variation of axial velocity with changing $\beta$ is highlighted in Fig. 5(b).
This figure shows that for higher $\beta$, we obtain a low centerline
velocity, as compared to the low $\beta$ value. This is because as the $\beta$
increases, the elastic behaviour of the fluid diminishes, and at $\beta=1$,
the viscoelastic fluid becomes Newtonian fluid. The polymeric shear stress
$\widetilde{\tau}_{xy}$ variation along the $\widetilde{y}$ direction at
$\widetilde{x}=40$ for changing $De$ and $\beta$ is illustrated in Fig. 5(c,d)
implies that for a lower $\beta$ the magnitude of $\widetilde{\tau}_{xy}$ at
the channel wall is higher and zero at the centerline. The value of
$|\widetilde{\tau}_{xy}|$ reduces with increasing $\beta$ and becomes zero for
$\beta=1$, as the viscoelastic fluid behaves as a Newtonian fluid. The smooth
variation of $\widetilde{\tau}_{xy}$ along $\widetilde{y}$ is altered with
increased fluid elasticity which is evident from the comparison of
$\widetilde{\tau}_{xy}$ vs $\widetilde{y}$ plots for $De=1$ and $De=50$. With
increasing $De$, the fluid relaxation time increases, and with reducing
$\beta$, the polymeric component of the viscosity increases. Thus, the elastic
behavior of the viscoelastic fluid dominates over the viscous behavior near
the walls at higher $De$ and lower $\beta$. This makes the flow
characteristics more complex and as a result, the polymeric stress profiles
become more bump near the walls.

###  4.2 Effect of amplitude of pulsation

In the pulsatile flow of viscoelastic fluid, the amplitude of pulsation ($A$)
is a major flow governing parameter, which can significantly alter the flow
dynamics. Hence, in Fig. 6, we have presented the variation of velocity and
stress components for changing $A$. It is evident from Fig. 6(a,b) that with
the increase in pulsation amplitude, the centerline velocity decreases because
of the interplay of continuous fluctuation in axial velocity and the elastic
nature of the viscoelastic fluid. It is interesting to note that for
$\beta=0.9$, which denotes less elastic behaviour of fluid than $\beta=0.1$,
application of equal amplitude of axial velocity fluctuation gives minimal
deviations in the velocity distributions as evident from Fig. 8(a). Thus, it
is clear that for high $De$ or low $\beta$ the alterations in fluctuation
amplitude affects the velocity distribution considerably. The decrease in the
centerline velocity for higher $A$ values is also evident from the variation
of $\widetilde{U}_{x}$ along the axial direction, as illustrated in Fig. 6(b).
Figure 6(c-e) presents the variations in the time-averaged polymeric stress
components at $\widetilde{x}=40$ for changing pulsation amplitude. From this
figure, we observe that the magnitude of polymeric stresses is high at the
near-wall region, which is generated because of the large velocity gradients
at that region. We also observe that the magnitude of $\widetilde{\tau}_{xx}$
and $\widetilde{\tau}_{xy}$ are of $O(1)$, whereas the magnitude of
$\widetilde{\tau}_{yy}$ is of the order $\sim O(10^{-4})$ at the channel
walls.

###  4.3 Effect of pulsation frequency or effect of Womersley number

In this section, we highlight the effect of pulsation frequency, the
dimensionless representation of which is the Womersley number $Wo$, on the
overall flow dynamics. The alterations in the velocity and stress components
with varying $Wo$ is illustrated in Fig. 7. In the subplot (a) of Fig. 7, we
represent the axial velocity profiles along the transverse direction at
$\widetilde{x}=40$. The centerline axial velocity shows a maximum value for
$Wo=0$ or uniform flow condition, and it decreases marginally with the
increase in $Wo$. The reduction in centerline axial velocity can also be seen
in Fig. 7(b), which highlights the variation of $\widetilde{U}_{x}$ as a
function of $Wo$ along the axial direction. With an increase in $Wo$, the
frequency of pulsation increases, and the interplay of non-zero fluid
elasticity with the frequency of pulsation manifests retardation in the flow
velocity. This can be further clarified from Fig. 8(b), which shows the
velocity distribution for exact variations in pulsation frequency at a higher
$\beta=0.9$. The comparison between the two figures, i.e. Fig. 7(a) and Fig.
8(b) shows that for high $\beta$ or low $De$ the alterations in pulsation
frequency does not affect the velocity distribution considerably as the
elastic behaviour of the fluid diminishes at higher $\beta$ or lower $De$
values. From Fig. 7(b), it can be observed that at the entry of the channel,
the boundary layer develops, and the flow becomes fully developed after the
developing region. Thus, the axial velocity increases along the axial
direction up to the entry length and beyond which, it becomes constant. This
phenomenon is instantaneous due to the pulsatile nature of the flow. The time-
averaged polymeric stress components at $\widetilde{x}=40$ for changing
pulsation frequency can be seen in Fig. 7(c-e). From this figure, it is
evident that the near-wall region of high-velocity gradients is the regions
for high polymeric stresses and the magnitude of polymeric stresses increases
with increasing pulsation frequency or $Wo$. We also observe that the
magnitude of $\widetilde{\tau}_{xx}$ and $\widetilde{\tau}_{xy}$ are of
$O(1)$, whereas the magnitude of $\widetilde{\tau}_{yy}$ is of the order $\sim
O(10^{-7})$ at the channel walls. Detailed discussion on the effect of both
pulsation frequency and pulsation amplitude at high retardation ratio and high
Deborah number can be found form Appendix A and Appendix B, respectively.

##  5 CONCLUSIONS

To conclude, we have analyzed the pulsatile flow of viscoelastic fluid in a
parallel plate channel numerically to elucidate the combined effects of
pulsation parameters and complex rheology of the fluid on the flow
characteristics. The numerical solution is obtained by a fully implicit
finite-volume method (FVM) based on a time-marching pressure-correction
algorithm in OpenFOAM® framework. The grid independence study and the solver
validations were carried out to check the accuracy of obtained numerical
results. With increasing viscoelastic behaviour of the fluid, we observe a
significant enhancement in the magnitude of axial velocity and polymeric shear
stress at the near-wall region. The elasticity of the viscoelastic fluid can
be enhanced either by increasing the Deborah number $De$ or by reducing the
retardation ratio $\beta$. For a time-averaged pulsatile flow, we observe a
reduction in the axial velocity with increasing pulsation amplitude or
frequency. The extent of reduction in velocity and polymeric stress magnitude
with changing pulsation parameters is highly dependent on the elastic
properties of the viscoelastic fluid. The results of the present work may find
its application in the transport of polymeric solutions, extrusion, and
injection moulding of polymer melts in several process industries, as well as
transport of biofluids in microfluidic devices.

## Supplementary Material

The supplementary material contains a detailed explanation of the log-
confirmation tensor approach and the solution steps for grid independence
study using the GCI method for the present problem.

## Data availability statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

$a$ | amplitude of pulsation, (m/s)
---|---
$A$ | dimensionless amplitude
$\mathbf{D}$ | deformation rate tensor, (s-1)
$De$ | Deborah number
$H$ | half height of the channel, (m)
$L$ | length of the channel, (m)
$\mathbf{M}$ | arbitary second-order tensor
$n$ | number of cells
$N$ | total number of cells
$p$ | fluid pressure, (Pa)
$Re$ | Reynolds number
$t$ | time, (s)
$u,v$ | velocity components, (m/s)
$\mathbf{U}$ | fluid velocity, (m/s)
$Wi$ | Weisenberg number
$Wo$ | Womersley number
$x,y$ | space coordinates, (m)

### Greek symbols

$\beta$ | retardation ratio
---|---
$\eta$ | total viscosity of the fluid, (Pa.s)
$\eta_{p}$ | polymeric viscosity, (Pa.s)
$\eta_{s}$ | solvent viscosity, (Pa.s)
$\lambda$ | relaxation time, (s)
$\lambda_{r}$ | retardation time, (s)
$\omega$ | angular frequency, (s-1)
$\rho$ | fluid density, (kg/m3)
$\mathbf{\tau}$ | stress tensor, (Pa)

### Subscripts and superscripts

$0$ | mean or base
---|---
$max$ | maximum
$p$ | polymeric
$s$ | solvent
$T$ | transpose
$x,y$ | coordinate directions
$\sim$ | non-dimensional quantities

## Appendix

##  Appendix A Axial velocity distribution for changing pulsation parameters
at high retardation ratio

In the present flow configuration, the interplay of fluid elasticity and the
pulsation parameters gives rise to alterations in the velocity and stress
distributions. The variations in the pulsation parameters alone for a
Newtonian fluid will not show any deviation in the velocity or stress profiles
in a time-averaged flow field. Here, for a viscoelastic fluid having high
$\beta$ and low $De$ values, we have shown the axial velocity variation by
changing the amplitude and frequency of pulsation in Fig. 8. For $\beta=0.9$
and $De=0.1$, the amplitude of pulsation is varied from $A=0\to 0.5$ and we
observe small deviations in the centerline velocity only (refer to Fig. 8(a)),
which is expected because of the non-zero $De$. We have also altered the
pulsation frequency, i.e. $Wo=0\to 3.545$ keeping constant fluid elastic
properties and the deviations in the velocity profile is found to be
negligible (refer to Fig. 8(b)). From the above discussion, it is clear that
the existence of fluid elasticity enhances the effect of pulsation parameters
in a pulsatile flow of viscoelastic fluid.

##  Appendix B Axial velocity distribution for changing pulsation parameters
at high Deborah number

The viscoelastic fluid flow is highly dependent on the elastic behaviour of
the fluid which are non-dimensionally denoted by Deborah number $De$ and
retardation ratio $\beta$. The effect of changing $\beta$ on the axial
velocity distribution is previously discussed. Here, we focus on highlighting
the effect of high $De$ on the axial velocity distribution for changing $A$
and $Wo$. We have shown the axial velocity variation by changing the amplitude
and frequency of pulsation in subplots ($a$) and ($b$) of Fig. 9,
respectively. For a higher Deborah number i.e. $De=100$, we observe a
reduction in the centerline velocity with increasing the value of $A$ (refer
to Fig. 9(a)). However it is observed that the axial velocity is enhanced at
the near wall regions to satisfy the continuity. With alterations in the
pulsation frequency, i.e. $Wo=0\to 3.545$ for a higher $De$, it is found that
the higher $Wo$ is responsible for increasing the axial velocity at the near
wall region and decreases the centerline velocity as seen from Fig. 9(b). It
is clear from the above discussion that the effect of pulsation parameters in
a pulsatile flow of viscoelastic fluid is enhanced due to the increased
elastic property of the working fluid.

## References

  * [1] Jia, J.-h., and Hua, H.-x., 2008.  “Study of Oscillating Flow of Viscoelastic Fluid With the Fractional Maxwell Model”.  Journal of Fluids Engineering, 130(4), 04.  041201\.
  * [2] Shah, S. D., 1996.  Pulsating gas-assisted injection molding method and apparatus, Jan. 9.  US Patent 5,482,669.
  * [3] Wang, Q., and Jin-Ping, Q., 2007.  “Analysis on pulsating flow of melt for reciprocating extrusion in the screw channel”.  Polymer-Plastics Technology and Engineering, 46(11), pp. 1055–1061.
  * [4] Phan-Thien, N., 1978.  “On pulsating flow of polymeric fluids”.  Journal of Non-Newtonian Fluid Mechanics, 4(3), pp. 167–176.
  * [5] Sankar, D., and Lee, U., 2009.  “Mathematical modeling of pulsatile flow of non-newtonian fluid in stenosed arteries”.  Communications in Nonlinear Science and Numerical Simulation, 14(7), pp. 2971–2981.
  * [6] Phan-Thien, N., and Dudek, J., 1982.  “Pulsating flow of a plastic fluid”.  Nature, 296(5860), pp. 843–844.
  * [7] Jinping, Q., Dantao, C., Yanhong, F., Xiaochun, Y., and Jianhua, D., 2008.  “Influence of pulsating-pressure-induced injection molding on rheologic property of polymer melt in spiral cavity”.  Engineering Plastics Application, 4.
  * [8] Razavi, A., Shirani, E., and Sadeghi, M., 2011.  “Numerical simulation of blood pulsatile flow in a stenosed carotid artery using different rheological models”.  Journal of Biomechanics, 44(11), pp. 2021–2030.
  * [9] Tu, C., and Deville, M., 1996.  “Pulsatile flow of non-newtonian fluids through arterial stenoses”.  Journal of Biomechanics, 29(7), pp. 899–908.
  * [10] Hettiarachchi, H., Hsu, Y., Harris, T. J., and Linninger, A. A., 2011.  “The effect of pulsatile flow on intrathecal drug delivery in the spinal canal”.  Annals of Biomedical Engineering, 39(10), p. 2592.
  * [11] Dhananchezhiyan, P., and Hiremath, S. S., 2018.  “Improving the performance of micro pumps by reduction of flow pulsation for drug delivery application”.  Advances in Materials and Processing Technologies, 4(1), pp. 24–38.
  * [12] Mahapatra, B., and Bandopadhyay, A., 2020.  “Electroosmosis of a viscoelastic fluid over non-uniformly charged surfaces: Effect of fluid relaxation and retardation time”.  Physics of Fluids, 32(3), p. 032005.
  * [13] Bird, R., Armstrong, R., and Hassager, O., 1987.  Dynamics of Polymeric Liquids: vol.1. Fluid mechanics.  Wiley, New York.
  * [14] Giesekus, H., 1982.  “A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility”.  Journal of Non-Newtonian Fluid Mechanics, 11(1-2), pp. 69–109.
  * [15] Phan-Thien, N., and Tanner, R., 1983.  “Viscoelastic squeeze-film flows–maxwell fluids”.  Journal of Fluid Mechanics, 129, pp. 265–281.
  * [16] Khadrawi, A., Al-Nimr, M., and Othman, A., 2005.  “Basic viscoelastic fluid flow problems using the jeffreys model”.  Chemical engineering science, 60(24), pp. 7131–7136.
  * [17] Oldroyd, J., 1958.  “Non-newtonian effects in steady motion of some idealized elastico-viscous liquids”.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 245(1241), pp. 278–297.
  * [18] Phan-Thien, N., 1978.  “A nonlinear network viscoelastic model”.  Journal of Rheology, 22(3), pp. 259–283.
  * [19] Bird, R., Dotson, P., and Johnson, N., 1980.  “Polymer solution rheology based on a finitely extensible bead—spring chain model”.  Journal of Non-Newtonian Fluid Mechanics, 7(2-3), pp. 213–235.
  * [20] McLeish, T., and Larson, R., 1998.  “Molecular constitutive equations for a class of branched polymers: The pom-pom polymer”.  Journal of Rheology, 42(1), pp. 81–110.
  * [21] Amoreira, L., and Oliveira, P., 2010.  “Comparison of different formulations for the numerical calculation of unsteady incompressible viscoelastic fluid flow”.  Adv. Appl. Math. Mech, 4, pp. 483–502.
  * [22] Fattal, R., and Kupferman, R., 2004.  “Constitutive laws for the matrix-logarithm of the conformation tensor”.  Journal of Non-Newtonian Fluid Mechanics, 123(2-3), pp. 281–285.
  * [23] Fattal, R., and Kupferman, R., 2005.  “Time-dependent simulation of viscoelastic flows at high weissenberg number using the log-conformation representation”.  Journal of Non-Newtonian Fluid Mechanics, 126(1), pp. 23–37.
  * [24] Afonso, A., Oliveira, P. J., Pinho, F., and Alves, M., 2009.  “The log-conformation tensor approach in the finite-volume method framework”.  Journal of Non-Newtonian Fluid Mechanics, 157(1-2), pp. 55–65.
  * [25] Chen, X., Marschall, H., Schäfer, M., and Bothe, D., 2013.  “A comparison of stabilisation approaches for finite-volume simulation of viscoelastic fluid flow”.  International Journal of Computational Fluid Dynamics, 27(6-7), pp. 229–250.
  * [26] Mahapatra, B., and Bandopadhyay, A., 2021.  “Numerical analysis of combined electroosmotic-pressure driven flow of a viscoelastic fluid over high zeta potential modulated surfaces”.  Physics of Fluids, 33(1), p. 012001.
  * [27] Mahapatra, B., and Bandopadhyay, A., 2021.  “Microconfined electroosmotic flow of a complex fluid with asymmetric charges: Interplay of fluid rheology and physicochemical heterogeneity”.  Journal of Non-Newtonian Fluid Mechanics, 289, p. 104479.
  * [28] Perera, M., and Walters, K., 1977.  “Long-range memory effects in flows involving abrupt changes in geometry: Part i: flows associated with i-shaped and t-shaped geometries”.  Journal of Non-Newtonian Fluid Mechanics, 2(1), pp. 49–81.
  * [29] Sun, J., Phan-Thien, N., and Tanner, R. I., 1996.  “An adaptive viscoelastic stress splitting scheme and its applications: $avss/si$ and $avss/supg$”.  Journal of non-newtonian fluid mechanics, 65(1), pp. 75–91.
  * [30] Duarte, A., Miranda, A. I., and Oliveira, P. J., 2008.  “Numerical and analytical modeling of unsteady viscoelastic flows: The start-up and pulsating test case problems”.  Journal of non-newtonian fluid mechanics, 154(2-3), pp. 153–169.
  * [31] Rosenfeld, M., and Einav, S., 1995.  “The Effect of Constriction Size on the Pulsatile Flow in a Channel”.  Journal of Fluids Engineering, 117(4), 12, pp. 571–576.
  * [32] Daidzic, N. E., 2014.  “Application of Womersley Model to Reconstruct Pulsatile Flow From Doppler Ultrasound Measurements”.  Journal of Fluids Engineering, 136(4), 02.  041102\.
  * [33] Budwig, R., Egelhoff, C. J., and Tavoularis, S., 1997.  “Laminar Pulsatile Flow Through an Axisymmetric Sudden Expansion”.  Journal of Fluids Engineering, 119(1), 03, pp. 208–211.
  * [34] Yapici, K., Karasozen, B., and Uludag, Y., 2012.  “Numerical Analysis of Viscoelastic Fluids in Steady Pressure-Driven Channel Flow”.  Journal of Fluids Engineering, 134(5), 05.  051206\.
  * [35] Cable, P., and Boger, D., 1979.  “A comprehensive experimental investigation of tubular entry flow of viscoelastic fluids: Part iii. unstable flow”.  AIChE Journal, 25(1), pp. 152–159.
  * [36] Khabakhpasheva, E., Popov, V., Kekalov, A., and Mikhailova, E., 1989.  “Pulsating flow of viscoelastic fluids in tubes”.  Journal of non-newtonian fluid mechanics, 33(3), pp. 289–304.
  * [37] Letelier, M. F., Siginer, D. A., and Cáceres, C., 2002.  “Pulsating flow of viscoelastic fluids in straight tubes of arbitrary cross-section—part i: longitudinal field”.  International Journal of Non-Linear Mechanics, 37(2), pp. 369–393.
  * [38] Siginer, D. A., and Letelier, M. F., 2002.  “Pulsating flow of viscoelastic fluids in straight tubes of arbitrary cross-section—part ii: secondary flows”.  International journal of non-linear mechanics, 37(2), pp. 395–407.
  * [39] Rachid, H., and Ouazzani, M., 2014.  “Interaction of pulsatile flow with peristaltic transport of a viscoelastic fluid: case of a maxwell fluid”.  International Journal of Applied Mechanics, 6(05), p. 1450061.
  * [40] Javadzadegan, A., Esmaeili, M., Majidi, S., and Fakhimghanbarzadeh, B., 2009.  “Pulsatile flow of viscous and viscoelastic fluids in constricted tubes”.  Journal of mechanical science and technology, 23(9), pp. 2456–2467.
  * [41] Bakhti, H., Azrar, L., and Dumitru, B., 2017.  “Pulsatile blood flow in constricted tapered artery using a variable-order fractional oldroyd-b model”.  Thermal Science, 21(1 Part A), pp. 29–40.
  * [42] Moyers-Gonzalez, M. A., Owens, R. G., and Fang, J., 2009.  “On the high frequency oscillatory tube flow of healthy human blood”.  Journal of non-newtonian fluid mechanics, 163(1-3), pp. 45–61.
  * [43] Sajid, M., Zaman, A., Ali, N., and Siddiqui, A. M., 2015.  “Pulsatile flow of blood in a vessel using an oldroyd-b fluid”.  International Journal of Nonlinear Sciences and Numerical Simulation, 16(5), pp. 197–206.
  * [44] Qamar, A., Seda, R., and Bull, J. L., 2011.  “Pulsatile flow past an oscillating cylinder”.  Physics of fluids, 23(4), p. 041903.
  * [45] Gupta, S., Patel, S., and Chhabra, R., 2020.  “Pulsatile flow of power-law fluids over a heated cylinder: Flow and heat transfer characteristics”.  International Journal of Thermal Sciences, 152, p. 106330.
  * [46] Vamerzani, B., Norouzi, M., and Firoozabadi, B., 2014.  “Analytical solution for creeping motion of a viscoelastic drop falling through a newtonian fluid”.  Korea-Australia Rheology Journal, 26(1), pp. 91–104.
  * [47] Alves, M. A., Oliveira, P. J., and Pinho, F. T., 2003.  “Benchmark solutions for the flow of oldroyd-b and ptt fluids in planar contractions”.  Journal of Non-Newtonian Fluid Mechanics, 110(1), pp. 45–75.
  * [48] Kiran Raj, M., DasGupta, S., and Chakraborty, S., 2019.  “Biomimetic pulsatile flows through flexible microfluidic conduits”.  Biomicrofluidics, 13(1).
  * [49] Pimenta, F., and Alves, M., 2016.  rheotool.  https://github.com/fppimenta/rheoTool.
  * [50] Oliveira, P. J., Pinho, F. d., and Pinto, G., 1998.  “Numerical simulation of non-linear elastic flows with a general collocated finite-volume method”.  Journal of Non-Newtonian Fluid Mechanics, 79(1), pp. 1–43.
  * [51] Alves, M., Oliveira, P., and Pinho, F., 2003.  “A convergent and universally bounded interpolation scheme for the treatment of advection”.  International journal for numerical methods in fluids, 41(1), pp. 47–75.
  * [52] Pimenta, F., and Alves, M., 2017.  “Stabilization of an open-source finite-volume solver for viscoelastic fluid flows”.  Journal of Non-Newtonian Fluid Mechanics, 239, pp. 85–104.
  * [53] Van Doormaal, J., and Raithby, G., 1984.  “Enhancements of the simple method for predicting incompressible fluid flows”.  Numerical heat transfer, 7(2), pp. 147–163.
  * [54] Celik, I. B., Ghia, U., Roache, P. J., Freitas, C. J., Coleman, H., Raad, P. E., Celik, Ì., Freitas, C., and Coleman, H., 2008.  “Procedure for estimation and reporting of uncertainty due to discretization in $\\{$CFD$\\}$ applications”.  Journal of Fluids Engineering, 130(078001).
  * [55] Rath, S., and Dash, S. K., 2019.  “Effect of horizontal spacing and tilt angle on thermo-buoyant natural convection from two horizontally aligned square cylinders”.  International Journal of Thermal Sciences, 146, p. 106113.
  * [56] Rath, S., 2021.  “Complex interplay of power-law rheology and non-oberbeck-boussinesq effects on natural convection heat transfer in a confined domain”.  International Journal of Heat and Mass Transfer, 176, p. 121462.
  * [57] Na, Y., and Yoo, J., 1991.  “A finite volume technique to simulate the flow of a viscoelastic fluid”.  Computational mechanics, 8(1), pp. 43–55.

###### List of Figures

  1. 1 (Color Online) Schematic of the problem; (a) axial pulsation representation of the flow for a particular case of unity frequency; (b) computational domain of the present flow configuration where the point $p(x,y)$ denots the probe location
  2. 2 (Color Online) (a) Grid resolution of the computational domain which is divided into two zones, (b) magnified view of the near wall grid region, (c) variation of axial velocity with changing mesh resolution along the transverse direction at $\widetilde{x}=40$ for $De=10$, $\beta=0.1$, $Wo=1.772$, and $A=0.5$
  3. 3 (Color Online) Validation with the analytical work of Na and Yoo [57] for a uniform flow having $Re=1$, $\beta=0.01$, and $De=0.1$, representing axial velocity and polymeric stress component profiles along the transverse direction at $\widetilde{x}=40$
  4. 4 (Color Online) Variation of axial velocity at a probe location ($\widetilde{x}=40,\widetilde{y}=+0.5$) in the channel with time for (a) varying $Wo$ at $A=0.1$, (b) varying $A$ at $Wo=2.507$; other simulation parameters are $\beta=0.1$ and $De=1$
  5. 5 (Color Online) (a) Variation of the axial velocity along the transverse for varying $De$ at $\widetilde{x}=40$ for $Wo=2.507$, $\beta=0.1$, and $A=0.1$; (b) Variation of the axial velocity along the transverse for varying $\beta$ at $\widetilde{x}=40$ for $Wo=2.507$, $De=50$, and $A=0.1$; (c) $\widetilde{\tau}_{xy}$ vs $\widetilde{y}$ for varying $\beta$ at $\widetilde{x}=40$, $A=0.25$, $De=1$, and $Wo=2.507$; (d) $\widetilde{\tau}_{xy}$ vs $\widetilde{y}$ for varying $\beta$ at $\widetilde{x}=40$, $A=0.25$, $De=50$, and $Wo=2.507$
  6. 6 (Color Online) Variation of velocity and stress componenets with changing $A$: (a) $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $Wo=2.507$; (b) $\widetilde{U}_{x}$ vs $\widetilde{x}$ at $\widetilde{y}=0$, $De=0.1$, $\beta=0.1$, and $Wo=2.507$; (c) variation of axial normal stress along the transverse direction, i.e., $\widetilde{\tau}_{xx}$ vs $\widetilde{y}$; (d) variation of shear stress along the transverse direction, i.e., $\widetilde{\tau}_{xy}$ vs $\widetilde{y}$; (e) variation of transverse normal stress along the transverse direction, i.e., $\widetilde{\tau}_{yy}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $Wo=2.507$
  7. 7 (Color Online) Variation of velocity and stress componenets with changing $Wo$: (a) $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $A=0.1$; (b) $\widetilde{U}_{x}$ vs $\widetilde{x}$ at $\widetilde{y}=0$, $De=0.1$, $\beta=0.1$, and $A=0.1$; (c) variation of axial normal stress along the transverse direction, i.e., $\widetilde{\tau}_{xx}$ vs $\widetilde{y}$; (d) variation of shear stress along the transverse direction, i.e., $\widetilde{\tau}_{xy}$ vs $\widetilde{y}$; (e) variation of transverse normal stress along the transverse direction, i.e., $\widetilde{\tau}_{yy}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $A=0.1$
  8. 8 (Color Online) (a) Variation of axial velocity with changing $A$, i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.9$, and $Wo=2.507$; (b) Variation of axial velocity with changing $Wo$, i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.9$, and $A=0.1$
  9. 9 (Color Online) (a) Variation of axial velocity with changing $A$, i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=100$, $\beta=0.1$, and $Wo=2.507$; (b) Variation of axial velocity with changing $Wo$, i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=100$, $\beta=0.1$, and $A=0.1$

###### List of Tables

  1. 1 Range of values of the non-dimensional quantities used in the numerical solution
  2. 2 Characteristics of blockMesh and grid independence test showing in terms of maximum peak velocity (Note: mesh in z-direction is having unity depth i.e., $n_{z}=1$ and grading value= 1)

## FIGURES

Figure 1: (Color Online) Schematic of the problem; (a) axial pulsation
representation of the flow for a particular case of unity frequency; (b)
computational domain of the present flow configuration where the point
$p(x,y)$ denots the probe location

Figure 2: (Color Online) (a) Grid resolution of the computational domain which
is divided into two zones, (b) magnified view of the near wall grid region,
(c) variation of axial velocity with changing mesh resolution along the
transverse direction at $\widetilde{x}=40$ for $De=10$, $\beta=0.1$,
$Wo=1.772$, and $A=0.5$

Figure 3: (Color Online) Validation with the analytical work of Na and Yoo
[57] for a uniform flow having $Re=1$, $\beta=0.01$, and $De=0.1$,
representing axial velocity and polymeric stress component profiles along the
transverse direction at $\widetilde{x}=40$

Figure 4: (Color Online) Variation of axial velocity at a probe location
($\widetilde{x}=40,\widetilde{y}=+0.5$) in the channel with time for (a)
varying $Wo$ at $A=0.1$, (b) varying $A$ at $Wo=2.507$; other simulation
parameters are $\beta=0.1$ and $De=1$

Figure 5: (Color Online) (a) Variation of the axial velocity along the
transverse for varying $De$ at $\widetilde{x}=40$ for $Wo=2.507$, $\beta=0.1$,
and $A=0.1$; (b) Variation of the axial velocity along the transverse for
varying $\beta$ at $\widetilde{x}=40$ for $Wo=2.507$, $De=50$, and $A=0.1$;
(c) $\widetilde{\tau}_{xy}$ vs $\widetilde{y}$ for varying $\beta$ at
$\widetilde{x}=40$, $A=0.25$, $De=1$, and $Wo=2.507$; (d)
$\widetilde{\tau}_{xy}$ vs $\widetilde{y}$ for varying $\beta$ at
$\widetilde{x}=40$, $A=0.25$, $De=50$, and $Wo=2.507$

Figure 6: (Color Online) Variation of velocity and stress componenets with
changing $A$: (a) $\widetilde{U}_{x}$ vs $\widetilde{y}$ at
$\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $Wo=2.507$; (b)
$\widetilde{U}_{x}$ vs $\widetilde{x}$ at $\widetilde{y}=0$, $De=0.1$,
$\beta=0.1$, and $Wo=2.507$; (c) variation of axial normal stress along the
transverse direction, i.e., $\widetilde{\tau}_{xx}$ vs $\widetilde{y}$; (d)
variation of shear stress along the transverse direction, i.e.,
$\widetilde{\tau}_{xy}$ vs $\widetilde{y}$; (e) variation of transverse normal
stress along the transverse direction, i.e., $\widetilde{\tau}_{yy}$ vs
$\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $Wo=2.507$

Figure 7: (Color Online) Variation of velocity and stress componenets with
changing $Wo$: (a) $\widetilde{U}_{x}$ vs $\widetilde{y}$ at
$\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $A=0.1$; (b)
$\widetilde{U}_{x}$ vs $\widetilde{x}$ at $\widetilde{y}=0$, $De=0.1$,
$\beta=0.1$, and $A=0.1$; (c) variation of axial normal stress along the
transverse direction, i.e., $\widetilde{\tau}_{xx}$ vs $\widetilde{y}$; (d)
variation of shear stress along the transverse direction, i.e.,
$\widetilde{\tau}_{xy}$ vs $\widetilde{y}$; (e) variation of transverse normal
stress along the transverse direction, i.e., $\widetilde{\tau}_{yy}$ vs
$\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$, $\beta=0.1$, and $A=0.1$

Figure 8: (Color Online) (a) Variation of axial velocity with changing $A$,
i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=0.1$,
$\beta=0.9$, and $Wo=2.507$; (b) Variation of axial velocity with changing
$Wo$, i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$,
$De=0.1$, $\beta=0.9$, and $A=0.1$

Figure 9: (Color Online) (a) Variation of axial velocity with changing $A$,
i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$, $De=100$,
$\beta=0.1$, and $Wo=2.507$; (b) Variation of axial velocity with changing
$Wo$, i.e., $\widetilde{U}_{x}$ vs $\widetilde{y}$ at $\widetilde{x}=40$,
$De=100$, $\beta=0.1$, and $A=0.1$

## TABLES

Table 1: Range of values of the non-dimensional quantities used in the numerical solution Symbols |  Name |  Definition |  Range of values
---|---|---|---
$Re$ |  Reynolds number |  $\frac{\rho U_{0}H}{\eta}$ |  $1$
$De$ |  Deborah number |  $\frac{\lambda U_{0}}{H}$ |  $0-100$
$\beta$ |  Retardation ratio |  $\frac{\lambda_{r}}{\lambda}=\frac{\eta_{s}}{\eta_{s}+\eta_{p}}$ |  $0-1$
$Wo$ |  Womersley number |  $H\sqrt{\frac{\omega\rho}{\eta}}$ |  $0-3.545$
$A$ |  Pulsation amplitude |  $\frac{a}{U_{0}}$ |  $0-0.5$

Table 2: Characteristics of blockMesh and grid independence test showing in terms of maximum peak velocity (Note: mesh in z-direction is having unity depth i.e., $n_{z}=1$ and grading value= 1) Mesh |  Zone 1 |  |  Zone 2 |  | $N_{c}$ | $\widetilde{U}_{x,max}$ | $\%~{}error$
---|---|---|---|---|---|---|---
| $n_{x}\times n_{y}$ | Grading | $n_{x}\times n_{y}$ | Grading |  |  |
M1 | $50\times 20$ | $[10~{}~{}0.5]$ | $50\times 20$ | $[10~{}~{}2]$ | 2000 | 1.3956 | -
M2 | $100\times 40$ | $[10~{}~{}0.5]$ | $100\times 40$ | $[10~{}~{}2]$ | 8000 | 1.4111 | 1.11
M3 | $150\times 60$ | $[10~{}~{}0.5]$ | $150\times 60$ | $[10~{}~{}2]$ | 18000 | 1.4148 | 0.26


2021

[1,2]Mrityunjay Kothari

[1,3]Tal Cohen

1]Department of Civil and Environmental Engineering, Massachusetts Institute
of Technology, 77 Massachusetts Ave, Cambridge, 02139, MA, USA

2]Department of Mechanical Engineering, University of New Hampshire, Durham,
03824, NH, USA

3]Department of Mechanical Engineering, Massachusetts Institute of Technology,
77 Massachusetts Ave, Cambridge, 02139, MA, USA

# The crucial role of elasticity in regulating liquid-liquid phase separation
in cells

<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [

###### Abstract

Liquid-liquid phase separation has emerged as a fundamental mechanism
underlying intracellular organization, with evidence for it being reported in
numerous different systems. However, there is a growing concern regarding the
lack of quantitative rigor in the techniques employed to study phase
separation, and their ability to account for the complex nature of the
cellular milieu, which affects key experimentally observable measures, such as
the shape, size and transport dynamics of liquid droplets. Here we bridge this
gap by combining recent experimental data with theoretical predictions that
capture the subtleties of nonlinear elasticity and fluid transport. We show
that within a biologically accessible range of material parameters, phase
separation is highly sensitive to elastic properties and can thus be used as a
mechanical switch to rapidly transition between different states in cellular
systems. Furthermore, we show that this active mechanically mediated mechanism
can drive transport across cells at biologically relevant timescales and could
play a crucial role in promoting spatial localization of condensates; whether
cells exploit such mechanisms for transport of their constituents, remains an
open question.

###### keywords:

Liquid-Liquid Phase Separation, Ostwald ripening, Elasticity, Front
Propagation

##  1 Introduction

Liquid-liquid phase separation (LLPS) has emerged as a fundamental mechanism
by which eukaryotic cells organize themselves into membraneless compartments
called biomolecular condensates that carry out important cellular functions
(brangwynne2009germline, ; hyman2014liquid, ; alberti2019considerations, ).
Despite its recognition over the last decade as a means of self-organization
in biology, there is an emerging skepticism regarding the supporting evidence
for LLPS (mcswiggen2019evaluating, ; Leslie336, ). A major criticism levelled
at the current approaches to study LLPS is their reliance on qualitative
measures and the lack of concomitant quantitative modeling. For instance,
Fluorescence Recovery After Photobleaching (FRAP) – a technique used commonly
to ascertain the liquid-like nature of biomolecular condensates – measures the
time it takes for photo bleached molecules to diffuse out of these liquid
droplets. However, the observed diffusion times from FRAP studies
(mcswiggen2019evaluating, ) show a huge spread (from a few seconds to minutes)
and are likely to be sensitive to the droplet’s complex environment (such as
the elasticity of biopolymer networks in the medium surrounding the droplet),
yet this connection is poorly understood. Further, qualitative approaches are
not infallible to experimental limitations (like diffraction limit), post-
processing induced artifacts, and presence of alternate mechanisms that lead
to similar outcomes. Given the small size of biomolecular condensates and the
complex nature of the physiological environment inside cells, purely
qualitative observations can be misleading without accounting for the various
physical effects that are significant at cellular length scales. Therefore,
there is a need to supplement the current methods of investigation with
rigorous quantitative analyses that establish the crucial link between the
medium’s material properties and the outcome of LLPS, and can help guide and
interpret future studies.

Recent experiments have revealed that the sea of elastic network contained in
cells – which impart them an elastic-solid like mechanical resistance to
deformation – play a very important role in LLPS by providing elastic
resistance to the growth of biomolecular condensates (zhang2020mechanical, ;
lee2021chromatin, ; boddeker2021non, ; fernandez2021putting, ). This mechanism
may hold answers to why LLPS inside the cellular milieu leads to tightly
regulated condensate sizes and associated dynamics, unlike the classic example
of oil-water phase separation.

Synthetic polymer based systems have provided a simpler droplet-in-polymer
analog to condensate-in-cell systems to experimentally and theoretically study
the effects of elastic resistance (style2018liquid, ; rosowski2020elastic, ;
rosowski2020elastic2, ; ronceray2021liquid, ; vidal2021cavitation, ;
mukherjee2021statistical, ; biswas2021thermodynamics, ; paulin2021fluid, ;
wei2020modeling, ). These systems have raised some intriguing questions.
Existing studies frequently treat droplet growth as a cavitation process
(ball1982discontinuous, ), which assumes a constant, size-independent pressure
inside the droplet, and thus does not energetically distinguish between a
single large droplet and multiple small droplets (overall droplet volume being
the same). Then why do we see the latter and not the former? Furthermore, how
does elasticity affect the commonly measured experimental quantities, such as
size of the droplets, their coarsening dynamics, and their spatial
localization?

In our previous work on liquid-liquid phase separation inside synthetic
polymers (kothari2020effect, ), we showed that the elastic resistance imposed
by the medium delays phase separation, arrests droplet growth, and can even
drive flux of liquid from regions of high stiffness to low stiffness. For the
specific problem at hand, our analysis in (kothari2020effect, ) was restricted
to the assumption of infinite domains and employed a simple transport law
where liquid flux was assumed to be driven by concentration gradients.

To quantitatively understand the role of elasticity in intracellular phase
separation, where the material properties and length scales differ
substantially from synthetic systems studied so far in the literature, there
is a need to port the knowledge from synthetic systems to biological systems
and to address the above limitations. In particular, in recent years there has
been an accumulation of evidence that strain stiffening plays an important
role in various functions of the cell li2021nonlinear ; han2018cell ;
van2016strain ; hu2019high . Could this effect be crucial, even at the
cellular length scales, in explaining the size regulation of membraneless
organelles? Motivated by these questions, we develop a theoretical model that
builds on (kothari2020effect, ), and apply it to explain experimental
observations in relevant cellular systems. We then extend the analysis beyond
the experimental range to broaden our understanding of the constitutive
sensitivities.

##  2 Model

The model presented in this section closely follows the authors’ previous work
(kothari2020effect, ) with a key point of departure: in this work we employ a
generalized kinetic model by prescribing the flux as proportional to its
thermodynamic conjugate, i.e. gradients in chemical potential111The chemical
potential is defined here as the change in Gibbs free energy of a system when
a molecule is removed or added to it at constant temperature and pressure., as
opposed to prescribing the flux proportional to concentration gradients. This
kinetic law allows us to model a much broader range of material
heterogeneities and is capable of describing both downhill (from higher to
lower concentration) and uphill diffusion (from lower to higher
concentration), thereby widening the applicability of the theory to real-world
systems.

Our model system comprises an elastic network of cross-linked polymer,
permeated by a liquid mixture that is made of the following two components:
free liquid of the corresponding un-crosslinked or partially crosslinked
polymer chains (jensen2015wetting, ) (denoted by A), and another liquid of a
different species (denoted by B). The model system is a simplified
representation of cellular milieus that show phase separation. For example, in
cells, proteins and RNAs are mixed inside the cytoplasm, which contains
elastic network as well as many other liquids. Proteins and RNAs phase
separate to create P granules, that segregate but coexist with cytoplasm
(brangwynne2009germline, ; hyman2014liquid, ). Similarly, inside the nucleus,
the nucleolar protein fibrillarin (FIB-1) phase separates from the nucleoplasm
to create nucleoli (berry2015rna, ).

Typically, there are two main time scales associated with phase separation in
heterogeneous materials. Short timescale: as the mixture (supersaturated with
liquid B) is stimulated, droplets of liquid B nucleate and grow inside the
elastic matrix before achieving a quasi-equilibrium size, which is governed by
a balance between the local chemical and elastic properties of the mixture. If
the nucleation and growth process is rapid, the global heterogeneity does not
influence the local behavior of the droplets, which involves only short range
migration of liquid and thus only depends on local elastic properties. Long
timescale: heterogeneity can dominate transport and drive long-range migration
of liquid across the material at a longer timescale. In the current work, we
will restrict our attention to scenarios where the timescales are well
separated. This allows us to treat nucleation and growth separately from the
ensuing dynamics.

We make the following assumptions in our model. We consider situations in
which the matrix does not swell appreciably in the entire process, and thus
changes in elastic energy due to the swelling of the matrix can be neglected.
Accordingly, we model the process as mixing between liquids A and B, where the
cross-linked network only provides elastic resistance to droplet growth, and
only liquid B can migrate spatially. Restricting our attention to the dilute
limit, we also ignore the strain energy due to elastic interaction among
droplets, and we assume the droplet distribution to be monosdisperse.

The concentration of liquid B in the mixture is denoted by $\phi$, and is
defined with respect to the matrix volume (i.e. excluding the droplet volume).
Concentration of the liquid in the droplet phase is denoted by
$\phi_{D}=\frac{4\pi}{3}r^{3}n_{d}$, where $r$ is the radius of the droplets,
and $n_{d}$ is the number of droplets per unit volume; $\phi_{D}$ is defined
with respect to the total volume within a representative element.

As a supersaturated mixture of initial concentration $\phi^{sup}$ phase-
separates, the droplets nucleate and start to grow in the elastic matrix.
Phase separation is driven by the lowering of the mixing free energy
($\Delta\bar{G}_{mix}$). At later stages of phase separation, reduction of
surface energy ($\Delta\bar{G}_{sur}$) drives coarsening; however, the growth
and coarsening of droplets incurs substantial elastic energy
($\Delta\bar{G}_{el}$). This energy competition drives the system to choose a
droplet size that minimizes the total free energy. In terms of driving forces,
at equilibrium, the chemical potential of liquid B $(\mu)$ in the matrix
balances the chemical potential of the droplet $(\mu_{D})$.

The change in total free energy density of the system, calculated per unit
volume, can be expressed as the sum of contributions from mixing, surface, and
elastic energies

$\Delta\bar{G}(\phi,\phi_{D})=\Delta\bar{G}_{mix}(\phi,\phi_{D})+\Delta\bar{G}_{sur}(\phi_{D})\\\
+\Delta\bar{G}_{el}(\phi_{D}),$ (1)

respectively. In this work we will use the Flory–Huggins solution theory
(flory1953principles, ) to represent changes in mixing free energy,222We
assume that the number of lattice sites occupied by a single chain of the
liquid $N_{A}\gg 1$ for simplicity.

$\Delta\bar{G}_{mix}=\frac{kT}{\nu_{m}}(1-\phi_{D})(\phi\ln\phi+\chi(T)\phi(1-\phi)),$

where $k$ denotes the Boltzmann constant, $T$ denotes the temperature,
$\nu_{m}$ denotes the molecular volume of liquid B, and $\chi$ denotes the
Flory-Huggins parameter, which can be inferred by energy minimization for a
given saturation concentration of liquid B (see Appendix A). Changes in
surface energy are written as $\Delta\bar{G}_{sur}(\phi_{D})=4\pi r^{2}\Gamma
n_{d}$, where $\Gamma$ is the surface energy between the two liquids, and
$\Delta\bar{G}_{el}(\phi_{D})=\frac{4\pi}{3}r^{3}W(r)n_{d}$, where $W(r)$
denotes the elastic strain energy density due to growth of a single droplet
(see Appendix B).

In the short timescale, before long-range migration occurs, the local
concentration of the absorbed liquid, must be conserved, namely

$\phi^{sup}=\phi(1-\phi_{D})+\phi_{D},$ (2)

where $\phi^{sup}$ denotes the concentration of the supersaturated mixture
prior to phase separation. With the specific response functions chosen, the
short-timescale equilibrium droplet size can now be calculated by minimizing
the total free energy (1) subject to the mass conservation constraint (2).
Minimizing the free energy is equivalent to finding the concentration, $\phi$,
that equilibriates the chemical potential of the liquid in the matrix ($\mu$)
to the chemical potential of the liquid in the droplet ($\mu_{D}$). When the
droplets and matrix are in equilibrium (minimum energy state), exchange of a
molecule will not alter the total energy. In other words, their chemical
potentials are the same. The derivative ${\rm d}\Delta\bar{G}/{\rm d}\phi=0$
gives

$\displaystyle\begin{split}\underbrace{kT(\ln\phi+(1-\phi)+\chi(1-\phi)^{2})}_{\mu}=\quad\quad\\\
\underbrace{\nu_{m}\left(\frac{2\Gamma}{r}+W(r)+\frac{r}{3}W^{\prime}(r)\right)}_{\mu_{D}},\end{split}$
(3)

where we have used ${\rm d}(\phi(1-\phi_{D})+\phi_{D})=0$ from (2) (see
Appendix C for details). Here, $\mu$ is due to changes in free energy of the
mixed state with changes in the concentration, while $\mu_{D}$ arises due to
change in droplets’ surface energy and elastic energy upon change in its size.

The solution of equation (3) along with mass conservation constraint (2), for
each homogeneous region in the problem, fully describes the quasi-equilibrium
state of the system, and serves as the initial condition for the long
timescale dynamics of heterogeneous systems.333For a detailed treatment of the
short timescale behavior from nucleation to growth of droplets, we refer the
readers to (kothari2020effect, ).

At the longer timescale, spatial heterogeneities in material properties can
give rise to long-range migration in phase-separated systems. Rather than
modeling the individual droplets, we consider a 1D aggregate representation
with variations in droplet size and liquid concentration in the matrix
represented by the two field variables $\phi_{D}(x)$ and $\phi(x)$,
respectively. Such a 1D description is consistent with the symmetry of the
problems we will study in this work. Next, prescribing the flux $J$ of liquid
B as proportional to its thermodynamic conjugate, i.e. gradients in chemical
potential, we write

$J=-\frac{D\phi}{kT}\bf{\nabla}\mu,$ (4)

where $D$ is the diffusivity of the liquid (upon linearization (4) reduces to
the commonly used Fick’s law).

As the liquid starts to migrate across the matrix, the droplet-matrix
equilibrium is disturbed $(\mu_{D}\neq\mu)$ and a cascading effect follows
where droplet size may increase or decrease. This can be modeled by thinking
of droplets as a source/sink term which gives the following equations for the
dynamics of the system,

$\displaystyle\frac{\partial\phi(x,t)}{\partial t}+{\rm Div}J=s(x,t)$ (5)
$\displaystyle\frac{\partial\phi_{D}(x,t)}{\partial t}=-s(x,t).$ (6)

Here, $s(x,t)$ is the source term that captures the behavior of the droplets,
which can dissolve back into the matrix to replenish it or grow in size by
absorbing the excess liquid in the matrix. This process is driven by the
difference in chemical potential of liquid B; a thermodynamically consistent
form of the source term is thus chosen as,

$\displaystyle
s(x,t)=\frac{K}{\nu_{m}}\left(\mu_{D}-\mu\right)H(\phi_{D}),\quad\mathcal{K}=\frac{kTL^{2}K}{D\nu_{m}},$
(7)

where $H(\phi_{D})$ is the Heaviside function to ensure that the source is
exhausted at $\phi_{D}=0$, $L$ is a characteristic length scale of the system,
and $K$ is a material property, which in its non-dimensionalized form,
$\mathcal{K}$, is called the dissolution number — it quantifies the relative
eagerness of the droplet to give out the liquid to the matrix
(kothari2020effect, ). The specific value of $K$ is not known and estimating
it remains an area of future work. The kinetics of the process are governed
primarily by the dissolution number — a larger dissolution number makes the
source term bigger and makes the droplets respond faster to any changes in the
matrix.

Equations (4)-(7) complete the set of governing equations, which must be
supplemented with appropriate initial and boundary conditions to constitute a
well-defined problem. In the following, we will apply the model to two cases
of interest.

##  3 Results and Discussion

###  3.1 Case I - Size Regulation and Coarsening

The size of various cellular organelles must be tightly regulated for their
proper function. However, a distribution of droplets resulting from phase
separation typically tends to coarsen over time, where larger droplets grow at
the expense of smaller droplets — a process known as Ostwald ripening. How do
cells then achieve a stable distributions of multiple droplets? A few possible
explanations have been proposed. For instance, chemical reactions taking place
inside the droplets can potentially arrest the Ostwald ripening
(zwicker2015suppression, ). Alternatively, certain cellular components can act
as surfactants to stabilize the droplet size by reducing the surface energy of
the droplet-matrix interface (cuylen2016ki, ). In the current work, we explore
another mechanism that can stabilize the droplet distribution — the elastic
resistance to the growth and merger of droplets imposed by the medium. This
mechanism was considered for phase separation in synthetic polymers in
(style2018liquid, ; kothari2020effect, ). Here, we extend it to biopolymers
and show that this mechanism can still be important at the scales relevant for
the cells.

Biopolymers often show strain-stiffening behavior, where they progressively
stiffen as they are stretched. Using our model and specializing $W(r)$ for
strain-stiffening materials (see Appendix B), we study how elasticity
regulates the size of droplets. To illustrate the results, we select
representative material values for chromatin – a biopolymer found in nucleus –
from the literature (ronceray2021liquid, ; zhang2020mechanical, ;
berry2015rna, ). Since the values of strain-stiffening parameter $n$ and the
number density of droplets $n_{d}$ are not reported in these studies, we
select $n$ to represent moderate strain-stiffening (with more values reported
in the Appendix D), and choose $n_{d}$ to be proportional to the elastic
modulus, as motivated by the experimental finding from Style et al.
(style2018liquid, ). Accordingly, we set $n_{d}=\alpha E$, where the choice of
$\alpha=10^{13}$N-1m-1 is made to obtain droplet sizes in ballpark of those
reported in (zhang2020mechanical, ).

The droplet sizes at short timescales (or in absence of heterogeneity that can
lead to long timescale migration) are shown in Fig. 1 for two different
supersaturations and moderate strain-stiffening. In contrast to the commonly
assumed neo-Hookean response, in a strain-stiffening matrix the elastic energy
cost of expanding a single large droplet is higher than growing multiple small
droplets for the same total volume. Thus, strain-stiffening provides a clear
explanation for the observed multi-droplet state.

Figure 1: Typical droplet size as a function of chromatin stiffness at short
timescales. Material Properties: $n=0.95,\Gamma=5\times 10^{-7}$ Nm-1,
$\nu_{m}=1\times 10^{-23}$ m3, $\alpha=10^{13}$ N-1m-1, $\chi=2.14$. Inset:
optogenetically nucleated droplets in chromatin network show typical size of a
few microns (zhang2020mechanical, ). Figure 2: Elasticity-driven coarsening in
heterogeneous systems. Normalized droplet radius, $\tilde{R}$, for droplets at
the interface on the right side. Material Properties: $E=100$ Pa (Left);
$E=50$ Pa (Right); $\Gamma=5\times 10^{-7}$ Nm-1, $\nu_{m}=5\times 10^{-24}$
m3, $D=5\times 10^{-11}$m2s-1, $\alpha=10^{13}$ N-1m-1, $\chi=2.14$ (both
sides). Figure 3: Dynamics of coarsening in heterogeneous systems. Same system
as shown in the schematic in Fig. 2. Material Properties: $E=100$ Pa (Left);
$E=50$ Pa (Right); $\Gamma=5\times 10^{-7}$ Nm-1, $\nu_{m}=5\times 10^{-24}$
m3, $D=5\times 10^{-11}$m2s-1, $\alpha=10^{13}$ N-1m-1, $\chi=2.14$ (both
sides). (A) Rate of coarsening for different values of strain stiffening and
dissolution numbers; $\tilde{R}$ denotes the normalized droplet radius. Rate
of Ostwald ripening is shown for comparison in blue. (B) Spatial distribution
of droplet sizes for $t=120$ mins. Also see Appendix E for coarsening behavior
at longer times.

We find that even before onset of long timescale effects, the response of the
system is highly nonlinear in its dependence on the constitutive properties.
By choosing a reasonable value for $\alpha$, the droplet sizes obtained from
the model are in good agreement with the experimentally observed range of
droplet sizes zhang2020mechanical . This supports the hypothesis that elastic
resistance from the medium can play an important role in size regulation of
the droplets. However, experimental measurements of $n$ and $n_{d}$ are
required to provide accurate predictions. Beyond the intuitive trends that are
captured (i.e. that a stiffer medium results in smaller droplets), the rapid
decay of droplet radius, from a maximum value to a nearly unaffected size,
indicates a narrow zone of high sensitivity within a range of stiffnessses
that are accessible to active biological materials. Namely, biological systems
can use this process as a mechanical switch, whereby activating relatively
small changes in stiffness can induce large changes in droplet size. (Also see
Appendix D for strain-stiffening effects.)

While at the long timescale the surface tension driven coarsening is hindered
by the presence of elastic networks, these systems may show limited coarsening
driven by gradients in stiffness. To study this coarsening behavior in an
elastically heterogeneous medium, we take inspiration from
(rosowski2020elastic, ) and construct the simplest heterogeneous setup by
considering two homogeneous connected regions of differing stiffnesses, as
shown in Fig. 2. We require continuity of flux at the interface, and consider
and infinite domain (i.e. $x\in(-\infty,\infty)$). After the initial
nucleation and growth of droplets is complete, a longer timescale transport
emerges that causes the liquid to flow along the direction of decreasing
chemical potential.

Fig. 2 shows the results for the normalized droplet size at the interface,
denoted by $\tilde{R}$, evolving with time, where we use the initial droplet
radius – at the end of the short timescale – to normalize:
$\tilde{R}=r(x=0^{+},t)/r(x=0^{+},t=0)$. For strain-stiffening media, the
droplets eventually reach a final size (Figs. E3 and E4 in Appendix E show
coarsening plots up to 1200 minutes that show the size convergence). Figs.
3(A1-A3) show the rate of coarsening for a range of strain-stiffening
parameters and dissolution numbers, where the blue slopes show the Ostwald
ripening rate for comparison. At the initial times, the coarsening rate
follows the Ostwald ripening power law ($\tilde{R}^{3}-1\propto t$), and
diverges from it at longer times. Cases with $\mathcal{K}=10$ are the first to
show a significant deviation because of their faster dynamics, whereas the
cases with smaller $\mathcal{K}$ will take longer to diverge. In particular,
we find that even for moderate levels of strain-stiffening (i.e. Figs. 3(A2,
A3)), significant deviation from Ostwald ripening can occur and the coarsening
is eventually arrested. For a neo-Hookean material $(n=1)$, the coarsening
continues indefinitely. Nonetheless, it eventually becomes slower than Ostwald
ripening. Finally, Figs. 3(B1-B3) show the spatial distribution of droplet
sizes for all the different cases at $t=120$ minutes. Fig. 2 and Figs. E2-E4
(in Appendix E) also highlight the role of the dissolution number: since it
controls the dynamics of matrix-droplet transport, a higher dissolution number
leads to faster coarsening, and the droplets arrive at their final size faster
(in the case of a strain-stiffening medium).

###  3.2 Case II - Spatial Localization of Droplets in Finite Domains

Cells can spatially localize condensates by employing phase separation. For
instance, Brangwynne et al. (brangwynne2009germline, ) showed that in
Caenorhabditis elegans, P granules, starting from a uniformly distributed
state, localize to the posterior of the $\sim 50\mu m$-long cell over a period
of $\sim$ 10 minutes. Cells are also highly dynamic and regulate their
internal structure as well as mechanical properties both spatially and
temporally (heidemann2004towards, ). These two observations naturally raise
the question: how does the change of stiffness in the cellular medium impact
the localization of the droplets444While phase separation in
(brangwynne2009germline, ) was shown to be controlled by concentration of
polarity proteins, we use spatiotemporal localization of condensates as a
motivating example to explore the role of elasticity.? To answer this
question, we consider a scenario where the system develops a stiffness
gradient over time and study its impact on the localization of droplets.
Initially, the system is homogeneous and has a uniform distribution of
droplets, the size of which can be obtained by solving equations (2) and (3).
Over a timescale of minutes, the system develops a linear gradient in
stiffness, where the stiffness of the left end ($x=-L/2$) starts to increase,
while that of the right end ($x=L/2$) remains fixed: $E(x,t)=100-10(x/L-1/2)t$
Pa, where $t$ is the time in minutes and $x\in[-L/2,L/2].$ By establishing a
stiffness gradient, the system alters the chemical potential of the droplets,
which now exceeds that of the matrix, thus disturbing the equilibrium and
causing the liquid to migrate from left to right (along the decreasing
stiffness). We apply a no flux boundary condition on both the ends (i.e. $J=0$
at $x=-L/2,L/2$). Note that we take $n_{d}$ to be constant during this
process, corresponding to $E=100$Pa. In the initial process of nucleation and
growth, the regions near the droplet-matrix interface experience high strains
and, thus, inevitable damage. These locations also serve as preferential spots
for droplet growth. Furthermore, the elastic networks prevent the droplets
from coalescing. Thus, once the short timescale equilibrium has been achieved,
$n_{d}$ remains the same.

Figure 4: Spatial localization of droplets by stiffness gradient. A linear
stiffness gradient is established across a 50 $\mu$m region, which is
previously homogeneously stiff. Material Properties: $E(x,t)=100-10(x/L-1/2)t$
Pa where $t$ is in minutes, $n_{d}=10^{15}$ m-3, $\Gamma=5\times 10^{-7}$
Nm-1, $\nu_{m}=5\times 10^{-24}$ m3, $D=5\times 10^{-11}$m2s-1, $\chi=2.14$,
$n=0.95$.

Fig. 4 shows the results of this process over time for a $50\mu m$-long
domain. We find that the droplets localize over biologically relevant
timescales ($\sim$ 10 mins) to the right side. We also see that a higher
dissolution number promotes faster localization. When $\mathcal{K}$ is
sufficiently large, the droplet-matrix system can quickly adapt to changes in
the matrix stiffness and continuously maintain a quasi chemical equilibrium.
If $\mathcal{K}$ is too low, the droplet dynamics will lag behind the
stiffness changes in the matrix, and the system will always be out of chemical
equilibrium. This effect can be especially important in the tightly regulated
cellular environments. We speculate that the dissolution number lumps together
factors like local stress, swelling of the network, and diffused damage around
the droplet (kim2020extreme, ); a multi-scale approach can better resolve
these local dynamics and is beyond the scope of the current work.

##  4 Conclusion

In conclusion, these quantitative predictions in a simple and experimentally
realizable setting are aimed at uncovering the role of elastic driving forces
in intracellular phase separation. Even at the cellular length scales, the
elastic resistance from the cellular media can impact the growth and dynamics
of biomolecular condensates – both of which are critical to the functioning of
cells. While we have established that elasticity-driven mechanisms are
accessible to biological systems, the question of whether and to what extent
cells use these mechanisms remains an open question. Future work in this area
must focus on understanding if the changes in cellular stiffness are fast and
significant enough to affect the outcome of LLPS. The recent development of
advanced optogenetic techniques to selectively initiate phase separation
Spatio2017 , together with the capability to create spatially heterogeneous
crosslinking can be a promising technique to quantify the dynamics arising due
to differences in elastic properties, and to test the predictions of our
model. Finally, the model presented here is not without limitations. Future
work should extend beyond 1D geometries to reveal more complex phase
separation phenomena, and should account for additional physical mechanisms
that may become important in cellular systems, such as strain dependent
surface energy (jensen2017strain, ; xu2017direct, ).

Acknowledgments The authors acknowledge the support of the Office of Naval
Research, United States of America and Dr. Timothy B. Bentley, Program
Manager, under award number N00014-20-1-2561.

## Declarations

The authors declare no conflict of interests.

##  Appendix A Flory-Huggins Parameter

The Flory-Huggins parameter $\chi$ can be derived for a given solubility of
liquid B in the mixture, $\phi^{sat}$, as kothari2020effect ,

$\chi=-\frac{\log\phi^{sat}+(1-\phi^{sat})}{(1-\phi^{sat})^{2}},$ (8)

in the absence of any elastic and surface effects. Throughout this paper, we
use $\chi=2.14$ which translates to a solubility concentration of
$\phi^{sat}=0.06$.

##  Appendix B Strain Energy of Droplet Growing Inside a Strain-Stiffening
Elastic Matrix

To capture nonlinear material response at large strains, we use the
incompressible Mooney-Rivlin constitutive model for the elastic response of
the matrix (kothari2020effect, ). The work done in expanding a single droplet
from stress-free radius $r_{0}$ to an expanded radius $r$ is given as,

$\displaystyle\begin{split}W(r)=nE\bigg{(}\frac{5}{6}-\frac{r_{0}}{r}-\frac{r_{0}^{3}}{3r^{3}}+\frac{r_{0}^{4}}{2r^{4}}\bigg{)}+\\\
(1-n)E\bigg{(}\frac{r}{2r_{0}}-\frac{1}{3}-\frac{r_{0}^{2}}{r^{2}}+\frac{5}{6}\frac{r_{0}^{3}}{r^{3}}\bigg{)},\end{split}$
(9)

where $E$ is the stiffness of the crosslinked polymer, and $0\leq n\leq 1$ is
the strain-stiffening parameter; $n=1$ represents the neo-hookean material (no
strain stiffening) and the level of strain stiffening increases with
decreasing $n$. The stress-free radius is the length scale at which the
elastic strain energy dominates over the surface energy in the growth of the
droplet, and throughout this paper we choose $r_{0}=0.1\mu m$. As explained in
(kothari2020effect, ), $r_{0}$ may be different from the pore size of the
elastic network in the matrix. We also note that our choice of $r_{0}$ is
consistent with [6]: Fig. 3(b) in [6] shows that in a cross-linked matrix, the
pressure initially drops as the droplet grows, before increasing again when
the radius is $\sim 0.1\mu m$, which identifies the region of elastic
resistance from chromatin to the growth of droplets.

##  Appendix C Minimization of Free Energy

Following (1), the total free energy can be written as,

$\displaystyle\begin{split}\Delta\bar{G}(\phi,\phi_{D})=(1-\phi_{D})\frac{kT}{\nu_{m}}\\{\phi\ln\phi+\chi(T)\phi(1-\phi)\\}\\\
+4\pi r^{2}\Gamma{n_{d}}+\frac{4\pi}{3}r^{3}W(r){n_{d}}.\end{split}$ (10)

The system also satisfies the mass conservation constraint,

$\phi^{sup}=\phi(1-\phi_{D})+\phi_{D}.$ (11)

The energy minimization condition

$\frac{{\rm d}\Delta\bar{G}}{{\rm d}\phi}=0$ (12)

can then be evaluated and simplified using the following relations:
$\frac{\partial\left(1-\phi_{D}\right)}{\partial\phi}=\frac{(1-\phi_{D})}{(1-\phi)}$
and $\phi_{D}=\frac{4\pi}{3}r^{3}n_{d}$,

$\displaystyle\begin{split}&\frac{{\rm d}\Delta\bar{G}_{mix}}{{\rm
d}\phi}=\frac{kT}{\nu_{m}}\frac{\partial\left(1-\phi_{D}\right)}{\partial\phi}\\{\phi\ln\phi+\chi(T)\phi(1-\phi)\\}\\\
&+\frac{kT}{\nu_{m}}(1-\phi_{D})\\{1+\ln\phi+\chi(T)(1-2\phi)\\}\\\
=&\frac{kT}{\nu_{m}}\frac{\left(1-\phi_{D}\right)}{(1-\phi)}\\{\ln\phi+(1-\phi)+\chi(T)(1-\phi)^{2}\\}\\\
=&\frac{\left(1-\phi_{D}\right)}{(1-\phi)}\frac{\mu}{\nu_{m}}\end{split}$ (13)

$\displaystyle\begin{split}&\left(\frac{{\rm d}\Delta\bar{G}_{sur}}{{\rm
d}\phi}+\frac{{\rm d}\Delta\bar{G}_{el}}{{\rm d}\phi}\right)=\\\
&\frac{\partial}{\partial\phi_{D}}\left\\{\phi_{D}\left(\frac{3\Gamma}{r}+W(r)\right)\right\\}\frac{\partial\phi_{D}}{\partial\phi}\\\
=&-\frac{\left(1-\phi_{D}\right)}{(1-\phi)}\left\\{\frac{2\Gamma}{r}+W(r)+\frac{r}{3}W^{\prime}(r)\right\\}\\\
=&-\frac{\left(1-\phi_{D}\right)}{(1-\phi)}\frac{\mu_{D}}{\nu_{m}}\end{split}$
(14)

Finally, these equations yield

$\frac{{\rm d}\Delta\bar{G}}{{\rm
d}\phi}=\frac{\left(1-\phi_{D}\right)}{(1-\phi)}\left(\frac{\mu-\mu_{D}}{\nu_{m}}\right)=0$
(15)

as the equilibrium condition shown in (3).

##  Appendix D Equilibrium Droplet Size Variation with Strain-Stiffening

The equilibrium droplet size is determined by the solution of (3) together
with the mass conservation constraint (2). Using the form outlined in (9) for
$W(r)$, we study the sensitivity of the equilibrium droplet sizes to the level
of strain stiffening, as governed by the parameter $n$. Figure 5 shows that
the increasing strain-stiffening decreases the equilibrium droplet size.

Figure 5: Equilibrium droplet sizes for different levels of strain stiffening.

##  Appendix E Coarsening Plots for Longer Times

The following figures show normalized droplet radius, $\tilde{R}$, over time,
for droplets at the interface on the right side. The three figures are for
$n=1,0.95$ and $0.9$ respectively. Material Properties: $E=100$ Pa (Left);
$E=50$ Pa (Right); $\Gamma=5\times 10^{-7}$ Nm-1, $\nu_{m}=5\times 10^{-24}$
m3, $D=5\times 10^{-11}$ m2s-1, $\alpha=10^{13}$ N-1m-1, $\chi=2.14$ (both
sides) which translates to a solubility concentration of $\phi^{sat}=0.06$ in
the absence of elastic and surface effects. In strain-stiffening media, the
coarsening is arrested, leading to final size as is evident from Figs. E3 and
E4. The systems with smaller $\mathcal{K}$ take longer time to reach the final
state.

Figure 6: Coarsening dynamics for $n=1$. Figure 7: Coarsening dynamics for
$n=0.95$. Figure 8: Coarsening dynamics for $n=0.9$.

## References

  *   * (1) Brangwynne, C.P., Eckmann, C.R., Courson, D.S., Rybarska, A., Hoege, C., Gharakhani, J., Jülicher, F., Hyman, A.A.: Germline p granules are liquid droplets that localize by controlled dissolution/condensation. Science 324(5935), 1729–1732 (2009)
  * (2) Hyman, A.A., Weber, C.A., Jülicher, F.: Liquid-liquid phase separation in biology. Annual review of cell and developmental biology 30, 39–58 (2014)
  * (3) Alberti, S., Gladfelter, A., Mittag, T.: Considerations and challenges in studying liquid-liquid phase separation and biomolecular condensates. Cell 176(3), 419–434 (2019)
  * (4) McSwiggen, D.T., Mir, M., Darzacq, X., Tjian, R.: Evaluating phase separation in live cells: diagnosis, caveats, and functional consequences. Genes & development 33(23-24), 1619–1634 (2019)
  * (5) Leslie, M.: Separation anxiety. Science 371(6527), 336–338 (2021) https://science.sciencemag.org/content/371/6527/336.full.pdf. https://doi.org/10.1126/science.371.6527.336
  * (6) Zhang, Y., Lee, D.S., Meir, Y., Brangwynne, C.P., Wingreen, N.S.: Mechanical frustration of phase separation in the cell nucleus by chromatin. Physical review letters 126(25), 258102 (2021)
  * (7) Lee, D.S., Wingreen, N.S., Brangwynne, C.P.: Chromatin mechanics dictates subdiffusion and coarsening dynamics of embedded condensates. Nature Physics 17(4), 531–538 (2021)
  * (8) Böddeker, T.J., Rosowski, K.A., Berchtold, D., Emmanouilidis, L., Han, Y., Allain, F.H., Style, R.W., Pelkmans, L., Dufresne, E.R.: Non-specific adhesive forces between filaments and membraneless organelles. bioRxiv (2021)
  * (9) Fernández-Rico, C., Sai, T., Sicher, A., Style, R.W., Dufresne, E.R.: Putting the squeeze on phase separation. JACS Au (2021)
  * (10) Style, R.W., Sai, T., Fanelli, N., Ijavi, M., Smith-Mannschott, K., Xu, Q., Wilen, L.A., Dufresne, E.R.: Liquid-liquid phase separation in an elastic network. Physical Review X 8(1), 011028 (2018)
  * (11) Rosowski, K.A., Sai, T., Vidal-Henriquez, E., Zwicker, D., Style, R.W., Dufresne, E.R.: Elastic ripening and inhibition of liquid–liquid phase separation. Nature Physics, 1–4 (2020)
  * (12) Rosowski, K.A., Vidal-Henriquez, E., Zwicker, D., Style, R.W., Dufresne, E.R.: Elastic stresses reverse ostwald ripening. Soft Matter 16(25), 5892–5897 (2020)
  * (13) Ronceray, P., Mao, S., Košmrlj, A., Haataja, M.P.: Liquid demixing in elastic networks: Cavitation, permeation, or size selection? Europhysics Letters 137(6), 67001 (2022)
  * (14) Vidal-Henriquez, E., Zwicker, D.: Cavitation controls droplet sizes in elastic media. Proceedings of the National Academy of Sciences 118(40), 2102014118 (2021) 2102.02506. https://doi.org/10.1073/pnas.2102014118
  * (15) Mukherjee, S., Basu, A.: Statistical mechanics of phase transitions in elastic media. arXiv preprint arXiv:2103.06070 (2021)
  * (16) Biswas, S., Mukherjee, B., Chakrabarti, B.: Thermodynamics of droplets undergoing liquid-liquid phase separation. arXiv preprint arXiv:2104.00651 (2021)
  * (17) Paulin, O.W., Morrow, L.C., Hennessy, M.G., MacMinn, C.W.: Fluid–fluid phase separation in a soft porous medium. Journal of the Mechanics and Physics of Solids 164, 104892 (2022)
  * (18) Wei, X., Zhou, J., Wang, Y., Meng, F.: Modeling elastically mediated liquid-liquid phase separation. Physical Review Letters 125(26), 268001 (2020)
  * (19) Ball, J.M.: Discontinuous equilibrium solutions and cavitation in nonlinear elasticity. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 306(1496), 557–611 (1982)
  * (20) Kothari, M., Cohen, T.: Effect of elasticity on phase separation in heterogeneous systems. Journal of the Mechanics and Physics of Solids 145, 104153 (2020)
  * (21) Li, H., Zheng, Y., Han, Y.L., Cai, S., Guo, M.: Nonlinear elasticity of biological basement membrane revealed by rapid inflation and deflation. Proceedings of the National Academy of Sciences 118(11) (2021)
  * (22) Han, Y.L., Ronceray, P., Xu, G., Malandrino, A., Kamm, R.D., Lenz, M., Broedersz, C.P., Guo, M.: Cell contraction induces long-ranged stress stiffening in the extracellular matrix. Proceedings of the National Academy of Sciences 115(16), 4075–4080 (2018)
  * (23) Van Helvert, S., Friedl, P.: Strain stiffening of fibrillar collagen during individual and collective cell migration identified by afm nanoindentation. ACS applied materials & interfaces 8(34), 21946–21955 (2016)
  * (24) Hu, J., Li, Y., Hao, Y., Zheng, T., Gupta, S.K., Parada, G.A., Wu, H., Lin, S., Wang, S., Zhao, X., et al.: High stretchability, strength, and toughness of living cells enabled by hyperelastic vimentin intermediate filaments. Proceedings of the national Academy of Sciences 116(35), 17175–17180 (2019)
  * (25) Jensen, K.E., Sarfati, R., Style, R.W., Boltyanskiy, R., Chakrabarti, A., Chaudhury, M.K., Dufresne, E.R.: Wetting and phase separation in soft adhesion. Proceedings of the National Academy of Sciences 112(47), 14490–14494 (2015)
  * (26) Berry, J., Weber, S.C., Vaidya, N., Haataja, M., Brangwynne, C.P.: Rna transcription modulates phase transition-driven nuclear body assembly. Proceedings of the National Academy of Sciences 112(38), 5237–5245 (2015)
  * (27) Flory, P.J.: Principles of Polymer Chemistry. Cornell University Press, ??? (1953)
  * (28) Zwicker, D., Hyman, A.A., Juelicher, F.: Suppression of ostwald ripening in active emulsions. Physical Review E 92(1), 012317 (2015)
  * (29) Cuylen, S., Blaukopf, C., Politi, A.Z., Müller-Reichert, T., Neumann, B., Poser, I., Ellenberg, J., Hyman, A.A., Gerlich, D.W.: Ki-67 acts as a biological surfactant to disperse mitotic chromosomes. Nature 535(7611), 308–312 (2016)
  * (30) Heidemann, S.R., Wirtz, D.: Towards a regional approach to cell mechanics. Trends in cell biology 14(4), 160–166 (2004)
  * (31) Kim, J.Y., Liu, Z., Weon, B.M., Cohen, T., Hui, C.-Y., Dufresne, E.R., Style, R.W.: Extreme cavity expansion in soft solids: Damage without fracture. Science advances 6(13), 0418 (2020)
  * (32) Shin, Y., Berry, J., Pannucci, N., Haataja, M.P., Toettcher, J.E., Brangwynne, C.P.: Spatiotemporal Control of Intracellular Phase Transitions Using Light-Activated optoDroplets. Cell 168(1-2), 159–17114 (2017). https://doi.org/10.1016/j.cell.2016.11.054
  * (33) Jensen, K.E., Style, R.W., Xu, Q., Dufresne, E.R.: Strain-dependent solid surface stress and the stiffness of soft contacts. Physical Review X 7(4), 041031 (2017)
  * (34) Xu, Q., Jensen, K.E., Boltyanskiy, R., Sarfati, R., Style, R.W., Dufresne, E.R.: Direct measurement of strain-dependent solid surface stress. Nature communications 8(1), 1–6 (2017)