# Finite difference method in prolate spheroidal coordinates for freely
suspended spheroidal particles in linear flows of viscous and viscoelastic
fluids

Arjun Sharma  Sibley School of Mechanical and Aerospace Engineering, Cornell
University, Ithaca, NY, 14853, USA  Donald L. Koch  Robert Frederick Smith
School of Chemical and Biomolecular Engineering, Cornell University, Ithaca,
NY 14853, USA

###### Abstract

A finite difference scheme is used to develop a numerical method to solve the
flow of an unbounded viscoelastic fluid with zero to moderate inertia around a
prolate spheroidal particle. The equations are written in prolate spheroidal
coordinates, and the shape of the particle is exactly resolved as one of the
coordinate surfaces representing the inner boundary of the computational
domain. As the prolate spheroidal grid is naturally clustered near the
particle surface, good resolution is obtained in the regions where the
gradients of relevant flow variables are most significant. This coordinate
system also allows large domain sizes with a reasonable number of mesh points
to simulate unbounded fluid around a particle. Changing the aspect ratio of
the inner computational boundary enables simulations of different particle
shapes ranging from a sphere to a slender fiber. Numerical studies of the
latter particle shape allow testing of slender body theories. The mass and
momentum equations are solved with a Schur complement approach allowing us to
solve the zero inertia case necessary to isolate the viscoelastic effects. The
singularities associated with the coordinate system are overcome using
L’Hopital’s rule. A straightforward imposition of conditions representing a
time-varying combination of linear flows on the outer boundary allows us to
study various flows with the same computational domain geometry. For the
special but important case of zero fluid and particle inertia we obtain a
novel formulation that satisfies the force- and torque-free constraint in an
iteration-free manner. The numerical method is demonstrated for various flows
of Newtonian and viscoelastic fluids around spheres and spheroids (including
those with large aspect ratio). Good agreement is demonstrated with existing
theoretical and numerical results.

###### keywords:

Finite difference; prolate spheroidal coordinates; viscoelastic fluids;
moderate inertial effects; large aspect ratio fibers; spheres

††journal: Journal of Computational Physics

##  1 Introduction

The flow of viscoelastic or polymeric fluids around solid particles of various
shapes is important in many industrial processes. A particle shape is chosen
to achieve the manufactured product’s specific purpose or property. For
example, fibers allow the desired anisotropy in the roll-to-roll manufacturing
of high aspect ratio, low resistance films for flexible and transparent
electronics [1, 2]. In hydraulic fracturing [3], spheres may be used as
proppants to keep the pores of fractured rocks from closing. In extrusion
molding and fiber spinning [4, 5, 6, 7] spheres or fibers may be added to the
fluid to impart strength to the finished product. A spheroid is a convenient
shape to span a range of aspect ratios that can be synthesized by dispersing
polystyrene spheres in a solution of polyvinyl alcohol, followed by drying the
solution into thin sheets. The particles obtained after heating, stretching,
and cooling these sheets [8] may be used in experiments [9] or industrial
applications (such as the ones discussed above) involving the flow of
viscoelastic fluids around particles. Viscoelasticity arises in these fluids
from the underlying polymer molecules. Viscoelastic fluids exhibit several
properties such as shear thinning and first normal stress difference (rod-
climbing) due to the response of the polymers to the imposed flow field.
Therefore, numerical computations using specific polymer constitutive models
are very useful to isolate the origins of novel flow physics from the
interaction of the polymers with the suspended particles. Furthermore, due to
several parameters defining the properties of a viscoelastic fluid, such as
the polymer mobility, its maximum extensibility, and relaxation time [10],
computations complement laboratory experiments by exploring a wider range of
these parameters. Studies concerning moderate inertial effects such that the
flow is not turbulent but leads to mechanisms (absent in Stokes flow) such as
a significant Saffman lift force [11] on a particle are relevant to several
engineering applications and natural phenomena that involve particulate flows.
These include air and water pollution, pneumatic and slurry transport,
fluidized bed combustion, mineral separation, hemodynamics, and sedimentation
in rivers [12, 13].

The numerical challenges in studying the flow of particle suspensions in the
applications mentioned above are two-fold: fluid-particle interaction and
particle-particle interaction. To incorporate both these effects, numerical
computations involving more than one particle resort to the immersed boundary
method (IBM) [14, 15]. These are useful in studying dense particle
suspensions. However, the no-slip condition on the particle surface is not
directly imposed in IBM. Instead, the particle region is modeled using
fictitious forces required to enforce the necessary no-slip condition. The
diffuse particle-solid interface does not fully resolve the large polymer
stress gradients near the surface. In several of the scenarios above, the
particle concentration is dilute enough, so the particle-particle interaction
is rare. Thus, dilute particle suspensions can be modeled as an ensemble of
several isolated particles in an unbounded fluid. The fluid-particle
interaction is often analytical in Stokes flow (i.e. flow of inertia-less
Newtonian fluid). However, it is complex for a Newtonian fluid with moderate
inertia or a viscoelastic fluid because the particle-induced disturbance
affects the velocity field in a way that alters the interplay of viscous and
inertial or viscous and elastic forces. Therefore, valuable and more accurate
physical insight is obtained by studying the flow of a fluid around an
isolated particle where the no-slip on the particle surface is imposed. Such
numerical studies can also qualitatively complement IBM simulations (where the
particle-particle interaction is incorporated) for dense suspensions, such as
in [16]. Furthermore, accurate force and torque coefficients for spheroids of
different aspect ratios obtained from single particle simulations can be
useful in Lagrangian models. In this paper, we describe a numerical method
based on the finite-difference approximations to model the flow of
viscoelastic fluid around a prolate spheroidal particle. The equations are
solved in a particle-fixed reference frame which rotates and translates with
the spheroid.

Computational fluid mechanics of viscoelastic fluids, which first started in
the 1970s, is now a well-established research field [17, 18]. Viscoelasticity
of linear polymers is modeled through continuum equations governing second
moment of the polymer end-to-end distance averaged over the polymer
configuration [10]. Interesting and novel physical phenomena arise in the
industrially relevant parameter regimes when the polymer stretch is
significant, leading to significant polymer stress and its gradients. It is in
these parameter regimes where unique numerical challenges also arise that have
required ingenious solutions in the past, such as the log-conformation
formulation by Fattal and Kupferman (2004) [19] and algebraic numerical
treatment of the polymer stretch by Richter et al. (2010) [20] to ensure
finite polymer length and hence finite polymer stress. Accurately resolving
the polymeric flow around particles with large aspect ratios provides yet
another numerical challenge as it requires high spatial resolution to
accurately model the large polymer stress gradients along with the thin
particle surface. Therefore, this paper uses a prolate spheroidal coordinate
system to discretize the governing equations spatially. This exactly models
the particle surface as one of the coordinate surfaces and is well suited to
study the flow around a prolate spheroid, just as spherical [21] or
cylindrical [22, 23, 24] coordinate systems are beneficial in studying the
flow around a sphere and in a cylindrical pipe respectively. Furthermore, even
a uniform grid in our chosen coordinate system is naturally more clustered
(than a Cartesian grid) near the particle surface in the Euclidean sense,
allowing enhanced spatial resolution in the regions that require it the most.

The flow of viscoelastic fluids around particles in the aforementioned
industrial applications undergoes a series of local linear flows with time in
a Lagrangian reference frame. For example, in fiber spinning, the material is
first sheared within the spinneret and then pulled by the drawing mechanism
leading to a strong uni-axial extensional flow before solidifying to form a
fiber. The computational domain consists of the prolate spheroidal surface of
the solid particle as the inner boundary and a nearly spherical outer surface
where the imposed flow boundary conditions are applied. On the outer boundary
of our computational domain, we can apply any choice of imposed stationary,
time-varying, or alternating linear flow fields that mimic industrial
scenarios. This allows us greater flexibility in the choice of imposed flow
conditions as compared to previous numerical techniques where the
computational domain is problem specific such as using parallel, oppositely
moving walls to obtain simple shear flow [25, 26] or a cylindrical outer
surface for a uniform far-field flow [27]. In our method, boundary conditions
can be changed with time within a simulation.

Theoretical studies and numerical simulations are complementary. Theoretical
studies are however often done either for spherical [28, 29], or for slender
particles [30, 31]. For the latter, a matched asymptotic expansion in particle
aspect ratio, also known as slender body theory, [32, 33, 34] is used to
obtain useful physical insight by considering the particle aspect ratio to be
very large. Furthermore, flow around a particle is solved in an unbounded
fluid in such theoretical developments. The choice of a prolate spheroidal
coordinate system allows us to simulate the flow around a large aspect ratio
particle in a very large computational domain with good spatial resolution in
the regions near the particle surface where the velocity and polymer stress
gradients are expected to be large while maintaining a reasonable number of
mesh points. Therefore, our numerical method is a suitable testing ground for
slender body theories. As we demonstrate later, we can simulate the flow
around a sphere by allowing the particle aspect ratio to be 1+$\epsilon$,
where $\epsilon$ is a very small positive number.

We use a finite difference method to discretize the spatial gradients. Polymer
viscosity is large in polymer melts and concentrated polymer solutions. An
increase in the solvent viscosity in dilute polymer solutions leads to a large
polymer relaxation time, leading to interesting mechanisms that are
numerically challenging to resolve. Therefore, viscoelastic fluids are
generally highly viscous, and numerical studies with negligible fluid and
particle inertia are appropriate to study the effects of viscoelasticity [25,
35, 26]. Furthermore, the solutions of uniform and linear flows of an
unbounded inertia-less Newtonian fluid around spheres and spheroids [36] are
analytically known. Thus, numerical studies around such particles in the
absence of inertia are more relevant in the presence of viscoelasticity.

Some previous investigations of the flow of viscoelastic fluids, such as [25]
and [37], intending to ignore inertia, have used small but finite values of
Reynolds number, $Re$ (ratio of inertial to viscous forces). In such numerical
solvers, the momentum conservation and incompressibility (mass conservation)
equations are solved via a splitting method where momentum equations are
advanced in time. The incompressibility is imposed via a pressure Poisson
equation. Therefore, additional boundary conditions are required for the
pressure field, and Neumann boundary conditions are usually used at the solid
surface [38, 39]. The splitting method is more appropriate for large $Re$
flows where the splitting errors due to the introduction of artificial
boundary conditions are not dominant. However, when $Re$ is small, splitting
errors increase with the dominance of the viscous forces [40]. To avoid this
issue, we solve the coupled system of momentum and incompressibility equation
iteratively using a Schur complement method with GMRES [41]. In addition,
unlike the splitting method, we can solve a flow with $Re=0$ where the
momentum equation is quasi-steady, and the method used does not require us to
use time marching.

For a finite $Re$, we incorporate the inertial terms in the momentum equation
within the Schur complement method similar to [40]. Therefore, in addition to
solving the flow of inertia-less viscoelastic fluid, our method allows us to
access the effects of small to moderate inertia (with or without
viscoelasticity). Accurate simulations for small (but nonzero) inertia made
possible by the large computational domain enable us to test perturbation
theories for the limit $Re\ll 1$. Such simulations of a Newtonian fluid with
finite inertia for larger aspect ratio prolate spheroids can allow us to test
the slender body theories, such as in [42], describing the effect of inertia
on fibers.

In the next section, we present governing equations for the mass and the
momentum conservation in the fluid along with the polymer constitutive
equations in their original and log-conformal [19] formulation. We also
present the equations governing particle dynamics and treatment of the
boundary conditions. Section 3 deals with the temporal discretization of the
governing equations. This entails the description of the Schur complement
method [41] for solving mass and momentum equations, the technique of [20] to
treat finite extensibility of a polymer, quaternion formulation to account for
the particle orientation, and separate methodologies for zero and finite
particle inertia to account for the particle motion. In section 4 we
illustrate the spatial discretization of the equations using the finite
difference method. In this section, we also present the method to treat the
coordinate system dependent axis singularities that arise when the governing
equations are expressed in the prolate spheroidal coordinate system. This is
motivated by a similar method [23, 24] developed for the cylindrical
coordinate system. In section 5 we demonstrate the robustness and versatility
of our numerical solver through examples of a variety of flows past a particle
in (a) Stokes flow (inertia-less flows of a Newtonian fluid), (b) flow of a
Newtonian fluid with finite inertia, and (c) inertia-less flows of
viscoelastic fluids with different constitutive models. Comparison with
previous numerical studies or analytical results is provided for each case.
Our numerical solver is parallelized using the domain decomposition method
implemented in the Message Passing Interface (MPI). All the examples presented
are run on more than one processing unit. Finally, we present the conclusions
in section 6.

##  2 Governing equations

We consider the flow of an incompressible viscoelastic fluid around a prolate
spheroidal particle in a reference frame rotating and translating with the
particle (figure 1).

Fig. 1: Computational domain for the flow around a prolate spheroid: particle
surface, $\mathbf{r}_{p}$ is the inner boundary, and the exterior boundary is
$\mathbf{r}_{\infty}$ indicated with a dashed black curve.
$\mathbf{r}_{\infty}$ is at a large distance from the particle’s center and
represents a surface in the far-field where the velocity boundary conditions
are applied. A velocity field,
$\mathbf{r}\cdot\boldsymbol{\Gamma}+\mathbf{u}_{0}$ with time varying
$\boldsymbol{\Gamma}$ and $\mathbf{u}_{0}$ indicated by blue arrows and curves
can be imposed at $\mathbf{r}_{\infty}$. The gray shaded region is the
interior fluid region over which we solve the governing equations.

The equations governing the conservation of mass (incompressibility) and
momentum in this rotating reference frame are

$\nabla\cdot\mathbf{u}=0,$ (1)
$\rho_{f}\big{(}\frac{\partial\mathbf{u}}{\partial
t}+\textbf{ADV}(\mathbf{u},\mathbf{u}_{p},\boldsymbol{\omega}_{p};\mathbf{r})\big{)}=\nabla\cdot\boldsymbol{\sigma},$
(2)

where

$\textbf{ADV}(\mathbf{u},\mathbf{u}_{p},\boldsymbol{\omega}_{p};\mathbf{r})=\frac{d\mathbf{u}_{p}}{dt}+\mathbf{u}\cdot\nabla\mathbf{u}+2\boldsymbol{\omega}_{p}\times\mathbf{u}+\boldsymbol{\omega}_{p}\times\boldsymbol{\omega}_{p}\times\mathbf{r}+\frac{d\boldsymbol{\omega}_{p}}{dt}\times\mathbf{r}.$
(3)

$\boldsymbol{\omega}_{p}$ is the angular velocity of the particle,
$\mathbf{r}$ and $\mathbf{u}$ are the position and the fluid velocity vector
relative to the centroid of the particle, and $\boldsymbol{\sigma}$ is the
stress tensor field in the fluid. The stress in a viscoelastic fluid is the
sum of the Newtonian solvent stress, $\boldsymbol{\tau}$, (with viscosity
$\mu$) and the polymer, $\boldsymbol{\Pi}$, stress,

$\boldsymbol{\sigma}=\boldsymbol{\tau}+\boldsymbol{\Pi}=-p\boldsymbol{\delta}+2\mu\boldsymbol{e}+\boldsymbol{\Pi},$
(4)

where $p$ is the reduced pressure (the difference between the hydrodynamic and
the hydrostatic pressure),
$\boldsymbol{e}=(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})/2$ is the strain
rate tensor, and $\rho_{f}$ is the fluid density. A rotating frame fixed with
the particle avoids the need to introduce a mesh velocity but leads to the
first, third, fourth, and fifth terms on the RHS of equation (3). These arise
due to the non-inertial (rotating) reference frame and represent the linear
acceleration, Coriolis force, centrifugal force, and angular acceleration. As
shown in figure 1 the computational domain for these equations is bounded by
$\mathbf{r}_{p}$ (representing the particle surface) on the inside, and
$\mathbf{r}_{\infty}$ (representing the far-field) on the outside. The
boundary conditions on the velocity field are the imposed flow conditions in
the far-field and no-slip and no-penetration on the particle surface. In the
frame of reference rotating and translating with the particle, these are,

$\displaystyle\begin{split}\mathbf{u}=\mathbf{0},\text{ on particle
surface},\hskip
14.45377pt\mathbf{u}=\mathbf{u}_{\infty}(\mathbf{r})=\mathbf{r}\cdot\boldsymbol{\Gamma}+\mathbf{u}_{0}-\mathbf{u}_{p}-\boldsymbol{\omega}_{p}\times\mathbf{r},\text{
as }|\mathbf{r}|\rightarrow\infty\approx\mathbf{r}_{\infty},\end{split}$ (5)

where $\boldsymbol{\Gamma}$ and $\mathbf{u}_{0}$ are the imposed velocity
gradient and uniform velocity field. The angular and translational velocities
of the particle, $\boldsymbol{\omega}_{p}$ and $\mathbf{u}_{p}$, may either be
imposed or obtained from relevant equations governing the motion of the
particle due to the moments and forces acting on the particle. In the latter
case, Newton’s equations govern the particle motion,

$\rho_{p}V_{p}\frac{d\mathbf{u}_{p}}{dt}={\mathbf{f}_{\text{fluid}}+\mathbf{f}_{\text{ext.}}},$
(6)
$\mathbf{I}_{p}\cdot\frac{d\boldsymbol{\omega}_{p}}{dt}=\mathbf{q}_{\text{fluid}}+\mathbf{q}_{\text{ext.}},$
(7)

where $V_{p}$, $\rho_{p}$ and $\mathbf{I}_{p}$ are the volume, density and
moment of inertia tensor of the particle. $\mathbf{f}_{\text{ext.}}$ and
$\mathbf{q}_{\text{ext.}}$ are the external imposed force and torque that can
be prescribed. In case gravity is present, $\mathbf{f}_{\text{ext.}}$ includes
the buoyancy force acting on the particle due to the difference in particle
and fluid densities represented by $(\rho_{p}-\rho_{f})V_{p}{\mathbf{g}}$.
Here, ${\mathbf{g}}$ is the gravity vector in the particle fixed frame.
Although it is likely fixed in the inertial reference frame, the orientations
of $\mathbf{f}_{\text{ext.}}$, $\mathbf{q}_{\text{ext.}}$ and $\mathbf{g}$
vary in the particle reference frame with time for a rotating particle. For a
prolate spheroid with minor radius, $r_{\text{minor}}$ and aspect ratio,
$\kappa$,

$V_{p}=\frac{4}{3}\pi r_{\text{minor}}^{3}\kappa,\hskip
14.45377pt\mathbf{I}_{p}=\frac{1}{5}\rho_{p}V_{p}r_{\text{minor}}^{2}\begin{bmatrix}1+\kappa^{2}&0&0\\\
0&1+\kappa^{2}&0\\\ 0&0&2\end{bmatrix}.$ (8)

$\mathbf{f}_{\text{fluid}}=\mathbf{f}(\boldsymbol{\sigma})$ and
$\mathbf{q}_{\text{fluid}}=\mathbf{q}(\boldsymbol{\sigma})$ are the fluid
stress dependent hydrodynamic force and torque acting on the particle defined
as,

$\mathbf{f}(\boldsymbol{\sigma})=\int_{\mathbf{r}_{p}}dS\hskip
7.22743pt\boldsymbol{\sigma}\cdot\mathbf{n},\hskip
14.45377pt\mathbf{q}(\boldsymbol{\sigma})=\int_{\mathbf{r}_{p}}dS\hskip
7.22743pt\mathbf{r}\times(\boldsymbol{\sigma}\cdot\mathbf{n}),$ (9)

where $\mathbf{n}$ is the surface normal. In a viscoelastic fluid
$\mathbf{f}_{\text{fluid}}=\mathbf{f}(\boldsymbol{\sigma})$ and
$\mathbf{q}_{\text{fluid}}=\mathbf{q}(\boldsymbol{\sigma})$ are thus sums of
contributions from the Newtonian solvent and polymeric stress,

$\mathbf{f}_{\text{fluid}}=\mathbf{f}(\boldsymbol{\tau})+\mathbf{f}(\boldsymbol{\Pi}),\hskip
7.22743pt\mathbf{q}_{\text{fluid}}=\mathbf{q}(\boldsymbol{\tau})+\mathbf{q}(\boldsymbol{\Pi}),$
(10)

As considered in section 3.5.1 later, for finite particle inertia
($\rho_{p}\neq 0$) equations (6) and (7) are numerically integrated in time
after discretizing the time derivative on the left-hand side of these
equations. However, in the limit of zero particle inertia ($\rho_{p}=0$), the
appropriate conditions are zero net force and torque on the particle. In that
case, as discussed in section 3.5.2, either a secant iteration method (for a
massless particle in a fluid with finite density, $\rho_{f}$) or a novel
decomposition of the inertia-less non-Newtonian momentum equation combined
with a resistivity formulation is used to ensure that various components on
the right-hand side of these equations balance each other.

Finite difference discretization of the governing equations in the Cartesian
coordinate system is exactly satisfied by a uniform flow, but this is not the
case in curvilinear coordinates because the spatial gradients of the relevant
variables such as velocity, pressure, and polymer stress involve coordinate
system dependent metric derivatives [43]. Free-stream preservation of the
imposed linear flow field far away from the particle is particularly important
for the simulations of our interest since we use a large computational domain.
Various previously proposed techniques to treat this issue are in [43, 44] and
references therein. However, by simply simulating the deviation of the
velocity field from its far-field value, the discretization errors associated
with the violation of free-stream preservation are removed. In other words, we
do not need to numerically simulate the analytical value of the far-field
velocity, $\mathbf{u}_{\infty}$, and pressure, $p_{\infty}$.

Since the fluid is considered to be incompressible, in the imposed velocity
gradient tensor, $\boldsymbol{\Gamma}$, we require
$\text{trace}({\boldsymbol{\Gamma}})=0$ and hence we have
$\nabla\cdot\mathbf{u}_{\infty}=0$. The far-field flow (relative to particle
motion), $\mathbf{u}(\mathbf{r})_{\infty}=\mathbf{u}_{\infty}$ is linear in
position, the polymer stress generated by the linear imposed velocity field is
a spatially constant value, $\boldsymbol{\Pi}_{\infty}$. Therefore, the far-
field momentum equation is

$\frac{\partial\mathbf{u}_{\infty}}{\partial
t}+\textbf{ADV}(\mathbf{u}_{\infty},\mathbf{u}_{p},\boldsymbol{\omega}_{p};\mathbf{r})=-\frac{1}{\rho_{f}}\nabla
p_{\infty}.$ (11)

Hence, the governing mass and momentum equations for the deviation of velocity
field and the pressure field from the far-field flow,

$\widetilde{\mathbf{u}}=\mathbf{u}-\mathbf{u}_{\infty},\hskip
14.45377pt\widetilde{p}=p-p_{\infty}.$ (12)

are,

$\nabla\cdot\widetilde{\mathbf{u}}=0,$ (13)
${\rho_{f}}\big{(}\frac{\partial\widetilde{\mathbf{u}}}{\partial
t}+\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}},\mathbf{u}_{\infty},\boldsymbol{\omega}_{p};\mathbf{r})\big{)}=-\nabla\widetilde{p}+\nabla^{2}\widetilde{\mathbf{u}}+\nabla\cdot(\boldsymbol{\Pi}-\boldsymbol{\Pi}_{\infty}),$
(14)

where

$\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}},\mathbf{u}_{\infty},\boldsymbol{\omega}_{p};\mathbf{r})=\widetilde{\mathbf{u}}\cdot(\nabla\widetilde{\mathbf{u}}+\nabla\mathbf{u}_{\infty})+\mathbf{u}_{\infty}\cdot\nabla\widetilde{\mathbf{u}}+2\boldsymbol{\omega}_{p}\times\widetilde{\mathbf{u}},$
(15)

and the boundary conditions are

$\displaystyle\begin{split}\widetilde{\mathbf{u}}=-\mathbf{u}_{\infty}(\mathbf{r})=-\mathbf{r}\cdot\boldsymbol{\Gamma}-\mathbf{u}_{0}+\mathbf{u}_{p}+\boldsymbol{\omega}_{p}\times\mathbf{r},\text{
on particle surface },\hskip 14.45377pt\widetilde{\mathbf{u}}=0,\text{ as
}|\mathbf{r}|\rightarrow\infty\approx\mathbf{r}_{\infty}.\end{split}$ (16)

Simulating $\widetilde{\mathbf{u}}$ and $\widetilde{p}$ instead of
$\mathbf{u}$ and $p$ allows free-stream preservation trivially. Furthermore,
the momentum equation for $\widetilde{\mathbf{u}}$ does not include
centrifugal and angular acceleration terms (compare (3) with (15)) and hence
we do not need to numerically evaluate ${d\boldsymbol{\omega}_{p}}/{dt}$ at
different times during a simulation involving particle rotation.

To model the polymer stress, $\boldsymbol{\Pi}$, we consider various dumbbell
models that consider polymer molecules as a spring attached to Brownian beads
[10]. These dumbbell-based models define a constitutive equation for the
polymer configuration,
$\boldsymbol{\Lambda}=\langle\mathbf{qq}\rangle_{\text{polymer
configuration}}$, where $\mathbf{q}$ is the end-to-end vector of the dumbbell
and the angle brackets represent the average over polymer configurations.
$\boldsymbol{\Lambda}$ is non-dimensionalized with the square of the radius of
gyration of the polymer, and the dumbbell models have the following form of
the constitutive equation,

$\frac{\partial\boldsymbol{\Lambda}}{\partial
t}+\mathbf{u}\cdot\nabla(\boldsymbol{\Lambda}-\boldsymbol{\Lambda}_{\infty})=\nabla\mathbf{u}^{\text{T}}\cdot\boldsymbol{\Lambda}+\boldsymbol{\Lambda}\cdot\nabla\mathbf{u}{-\frac{1}{\lambda}\mathbf{R}(\boldsymbol{\Lambda})},$
(17)

where the polymer convection
$\mathbf{u}\cdot\nabla(\boldsymbol{\Lambda}-\boldsymbol{\Lambda}_{\infty})=\mathbf{u}\cdot\nabla\boldsymbol{\Lambda}$
is balanced by its stretching
($\nabla\mathbf{u}^{\text{T}}\cdot\boldsymbol{\Lambda}+\boldsymbol{\Lambda}\cdot\nabla\mathbf{u}$)
and relaxation (${\mathbf{R}(\boldsymbol{\Lambda})}/{\lambda}$) for a polymer
solution with relaxation time $\lambda$. Since any valid constitutive equation
is materially frame-invariant or objective [45], as long as the components of
the $\boldsymbol{\Lambda}$ tensor are expressed in the appropriate reference
frame, the constitutive equations retain their form in different reference
frames. In other words, unlike the momentum equation governing the fluid
velocity (a frame variant quantity), particle rotation does not introduce
additional terms in the constitutive equation governing $\boldsymbol{\Lambda}$
(a frame-invariant quantity) in a non-inertial or rotating reference frame. We
subtract the spatially constant far-field polymer configuration,
$\boldsymbol{\Lambda}_{\infty}$, before taking the gradient in the convective
term to ensure free-stream preservation. The polymeric stress,
$\boldsymbol{\Pi}$ is directly proportional to the polymer configuration,
$\boldsymbol{\Lambda}$ and can be written as

$\boldsymbol{\Pi}={\frac{c\mu}{\lambda}\mathbf{F}}(\boldsymbol{\Lambda}),$
(18)

where $c$ is the polymer concentration defined as the ratio of zero shear rate
polymer to solvent viscosity ($\mu$). The exact form of the relaxation tensor,
$\mathbf{R}(\boldsymbol{\Lambda})$, and the polymer force-configuration
relation, $\mathbf{F}(\boldsymbol{\Lambda})$, depends on the particular model
being used to represent the spring force. If the spring is considered to be
Hookean, one obtains the Oldroyd-B model; if the spring force is nonlinear and
has finite extensibility, one obtains a FENE (finite extensible nonlinear
elastic) model. A closure approximation is needed to obtain continuum-level
constitutive equations from the FENE model, and different choices lead to
FENE-P and FENE-CR models. The latter does not exhibit shear-thinning. The
Giesekus model is also a nonlinear model but has a quadratic nonlinearity. In
FENE models, the maximum length of the polymers is $L$ (non-dimensionalized
with the radius of gyration). The Giesekus model describes concentrated
polymer solutions and melts by including the anisotropic effect of nearby
dumbbells through the mobility parameter $\alpha$. See [10] for a review of
various constitutive models. $\mathbf{F}(\boldsymbol{\Lambda})$ and
$\mathbf{R}(\boldsymbol{\Lambda})$ for the different constitutive models we
consider are,

$\begin{matrix}\text{Model}:&\text{Oldroyd-B}&\text{FENE-P}&\text{FENE-
CR}&\text{Giesekus}\\\
{\mathbf{R}(\boldsymbol{\Lambda})}:&{\boldsymbol{\Lambda}-\boldsymbol{\delta}}&{f\boldsymbol{\Lambda}-b\boldsymbol{\delta}}&{f\boldsymbol{\Lambda}-f\boldsymbol{\delta}}&{(\boldsymbol{\Lambda}-\boldsymbol{\delta})-\alpha(\boldsymbol{\Lambda}-\boldsymbol{\delta})^{2}}\\\
\mathbf{F}(\boldsymbol{\Lambda}):&\boldsymbol{\Lambda}-\boldsymbol{\delta}&f\boldsymbol{\Lambda}-b\boldsymbol{\delta}&f\boldsymbol{\Lambda}-f\boldsymbol{\delta}&\boldsymbol{\Lambda}-\boldsymbol{\delta},\end{matrix}$
(19)

where ,

$f=1/(1-\text{tr}(\boldsymbol{\Lambda})/L^{2}),\text{ and
}b=1/(1-\text{tr}(\boldsymbol{\delta})/L^{2}).$ (20)

The polymer constitutive equation (17) (or equation (25) discussed later) is a
hyperbolic equation. Therefore, the boundary conditions are required at the
locations the streamlines of the velocity field, $\mathbf{u}$, enter the
computation domain. At these locations, the boundary condition is a time-
dependent, spatially constant polymer configuration tensor,
$\boldsymbol{\Lambda}=\boldsymbol{\Lambda}_{\infty}$, driven by the imposed
velocity field and is governed by,

$\frac{\partial\boldsymbol{\Lambda}_{\infty}}{\partial
t}=\nabla\mathbf{u}_{\infty}^{\text{T}}\cdot\boldsymbol{\Lambda}_{\infty}+\boldsymbol{\Lambda}_{\infty}\cdot\nabla\mathbf{u}_{\infty}{-\frac{1}{\lambda}\mathbf{R}(\boldsymbol{\Lambda}_{\infty})}.$
(21)

Using equation (18),
$\boldsymbol{\Pi}_{\infty}={c\mu/\lambda}\mathbf{F}(\boldsymbol{\Lambda}_{\infty})$
and $\mathbf{F}(\boldsymbol{\Lambda}_{\infty})$ is evaluated from equation
(19).

###  2.1 Log Conformal Form of Constitutive equations

An important property of $\boldsymbol{\Lambda}$ that is preserved by an
appropriate constitutive equation such as the ones introduced above is its
positive definiteness [19]. However, numerical discretization of constitutive
equations of the form in equations equation (17) and (18) with a model from
(19) may lead to violation of positive definiteness in the numerical solution
of $\boldsymbol{\Lambda}$ at high $De$. This manifests as a numerical
instability and is also termed the high Weissenberg number problem (HWNP)–
Weissenberg number is defined as the product of the characteristic flow
gradient magnitude and the polymer relaxation time; for steady linear flows
Deborah and Weissenberg numbers are equivalent. Fattal and Kupferman [19]
remedied the HWNP by introducing an equivalent constitutive equation for the
matrix logarithm of the conformation tensor, $\boldsymbol{\Lambda}$,

$\boldsymbol{\Psi}=\log(\boldsymbol{\Lambda}).$ (22)

Solving the governing equation for $\boldsymbol{\Psi}$ instead of
$\boldsymbol{\Lambda}$ provides a more stable numerical solution as found in
numerous numerical studies such as [25, 27, 46, 47, 48] after the seminal work
of [19]. Fattal and Kupferman [19] derived an equation for $\boldsymbol{\Psi}$
based on the eigen-decomposition of the velocity gradient. We use the
alternative derivation, provided by Hulsen and the previous authors [49],
based on the evolution of the principal axes of $\boldsymbol{\Lambda}$ (and
hence also $\boldsymbol{\Psi}$ since $\boldsymbol{\Lambda}$ and
$\boldsymbol{\Psi}$ have same eigenvectors). We find this form to be simpler
in treating cases such as biaxial extensional flow for which the two
eigenvalues are identical. We use the Jacobi algorithm provided in [50] to
obtain the eigen-decomposition of $\boldsymbol{\Lambda}$ and
$\boldsymbol{\Psi}$,

$\boldsymbol{\Lambda}=\mathbf{V}\cdot\mathbf{D}_{\boldsymbol{\Lambda}}\cdot\mathbf{V}^{T},\hskip
14.45377pt\boldsymbol{\Psi}=\mathbf{V}\cdot\mathbf{D}_{\boldsymbol{\Psi}}\cdot\mathbf{V}^{T},$
(23)

where $\mathbf{V}$ is a $3\times 3$ matrix with the eigenvectors,
$\mathbf{v}_{i},i\in[1,3]$ as its columns and
$\mathbf{D}_{\boldsymbol{\Lambda}}$ and

$\mathbf{D}_{\boldsymbol{\Psi}}=\log(\mathbf{D}_{\boldsymbol{\Lambda}}),$ (24)

are diagonal matrices with eigenvalues $\lambda_{i}$ and
$\psi_{i}=\log(\lambda_{i}),i\in[1,3]$ as their entries. The governing
equation for the matrix logarithm, $\boldsymbol{\Psi}$ is,

$\frac{\partial\boldsymbol{\Psi}}{\partial
t}+\mathbf{u}\cdot\nabla(\boldsymbol{\Psi}-\boldsymbol{\Psi}_{\infty})=\textbf{SR}(\boldsymbol{\Psi},\mathbf{u}),$
(25)

where,

$\textbf{SR}(\boldsymbol{\Psi},\mathbf{u})=2\Sigma_{i=1}^{3}L_{ii}\mathbf{v}_{i}\mathbf{v}_{i}+\Sigma_{i=1}^{3}\Sigma_{j=1,j\neq
i}^{3}\frac{\psi_{i}-\psi_{j}}{\lambda_{i}-\lambda_{j}}(\lambda_{j}L_{ij}+\lambda_{i}L_{ji})\mathbf{v}_{i}\mathbf{v}_{j}{-\exp(-\boldsymbol{\Psi})\cdot\mathbf{R}(\exp(\boldsymbol{\Psi}))},$
(26)

$\mathbf{L}=\nabla\mathbf{u}$ is the velocity gradient tensor,
$\exp(\boldsymbol{\Psi})=\mathbf{V}\cdot\mathbf{D}_{\boldsymbol{\Lambda}}\cdot\mathbf{V}^{T}$
and
$\exp(-\boldsymbol{\Psi})=\mathbf{V}\cdot\mathbf{D}_{\boldsymbol{\Lambda}}^{-1}\cdot\mathbf{V}^{T}$.
The three terms of $\textbf{SR}(\boldsymbol{\Psi},\mathbf{u})$ in equation
(26) represent the stretching of eigenvectors by $L_{ii}$, rotation of
eigenvectors by vorticity and their relaxation respectively [49]. When two
eigenvalues are identical [49] in the second term of equation (26),

$\lim\limits_{\lambda_{i}\rightarrow\lambda_{j}}\frac{\psi_{i}-\psi_{j}}{\lambda_{i}-\lambda_{j}}(\lambda_{j}L_{ij}+\lambda_{i}L_{ji})\rightarrow
L_{ij}+L_{ji}.$ (27)

The governing equation for the undisturbed matrix logarithm,
$\boldsymbol{\Psi}_{\infty}$, is obtained by setting
$\mathbf{L}=\nabla\mathbf{u}_{\infty}$ and
$\boldsymbol{\Psi}=\boldsymbol{\Psi}_{\infty}$ in equations (25) and (26).

##  3 Temporal discretization and coupling between equations and boundary
conditions

The constitutive equation (25) for the matrix logarithm, $\boldsymbol{\Psi}$
is driven by the velocity and velocity gradients. The polymer stress hence
generated (see equations (18), (19) and (22)) acts as a body force in the
momentum equation (2) and influences the velocity and pressure field. The
constraint of divergence-free velocity field defined by equation (1) ensures
mass conservation. In this section, we describe the temporal discretization
and methods that treat the coupling between mass (equation (1)) and momentum
(equation (2)) conservation and the polymer constitutive equation (25) . We
also describe the separate numerical treatments of the Newton’s equations
(equations (6) and (7)) governing the particle dynamics for the case when
particle inertia is negligible and when it is finite. The former case is of
particular interest when fluid inertia is also neglected.

We adopt a similar methodology as [40] to treat the time discretization of the
momentum and mass conservation equations. The mass conservation equation at
time step $n+1$ is

$\nabla\cdot\widetilde{\mathbf{u}}^{n+1}=0.$ (28)

Using a backward Euler temporal discretization scheme, at a time step $n+1$,
the momentum equation is written as,

$\displaystyle\begin{split}&\rho_{f}\big{(}\frac{\widetilde{\mathbf{u}}^{n+1}-\widetilde{\mathbf{u}}^{n}}{\Delta
t}+\frac{3}{2}\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}}^{n},\mathbf{u}_{\infty}^{n},\boldsymbol{\omega}_{p}^{n};\mathbf{r})-\frac{1}{2}\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}}^{n-1},\mathbf{u}_{\infty}^{n-1},\boldsymbol{\omega}_{p}^{n-1};\mathbf{r})\big{)}=\\\
&-\nabla\widetilde{p}^{n+1}+\nabla^{2}\widetilde{\mathbf{u}}^{n+1}+\nabla\cdot(\boldsymbol{\Pi}^{n+1}-\boldsymbol{\Pi}^{n+1}_{\infty}),\end{split}$
(29)

where the non-linear terms i.e.
$\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}},\mathbf{u}_{\infty},\boldsymbol{\omega}_{p};\mathbf{r})$
from equation (15) are treated explicitly using a second-order Adams–Bashforth
scheme after the first time step (in the first time step a first-order
explicit Euler scheme is used). Explicit treatment of these non-linear terms
dependent upon fluid inertia is a valid strategy because we aim to study the
effect of zero to moderate fluid inertia on flows of viscoelastic fluids.
After the first time step, we use a second-order implicit Crank-Nicholson
scheme where the polymer constitutive equation is temporally discretized as
following,

$\frac{\boldsymbol{\Psi}^{n+1}-\boldsymbol{\Psi}^{n}}{\Delta
t}+\frac{\mathbf{u}^{n+1}\cdot\nabla(\boldsymbol{\Psi}^{n+1}-\boldsymbol{\Psi}^{n+1}_{\infty})}{2}+\frac{\mathbf{u}^{n}\cdot\nabla(\boldsymbol{\Psi}^{n}-\boldsymbol{\Psi}^{n}_{\infty})}{2}=\frac{\textbf{SR}(\boldsymbol{\Psi}^{n+1},\mathbf{u}^{n+1})}{2}+\frac{\textbf{SR}(\boldsymbol{\Psi}^{n},\mathbf{u}^{n})}{2}.$
(30)

We consider both a fixed time step and a variable time step in our studies.
${\textbf{SR}(\boldsymbol{\Psi}^{n+1},\mathbf{u}^{n+1})}$ is dependent on the
unknown values at the current time step. We use a weighted Jacobi method to
iteratively solve equation (30) for $\boldsymbol{\Psi}^{n+1}$. In the first
time step we solve the polymer constitutive equation with a first order
implicit Euler method. To obtain $\boldsymbol{\Pi}^{n+1}$ from
$\boldsymbol{\Psi}^{n+1}$, we use the eigen-decomposition of
$\boldsymbol{\Psi}^{n+1}$ and then use equations (23) and (24) to obtain
$\boldsymbol{\Lambda}^{n+1}$ along with the relevant relation between
$\boldsymbol{\Lambda}^{n+1}$ and $\boldsymbol{\Pi}^{n+1}$ from equation (18)
and (19).

###  3.1 Decoupling mass-momentum system from polymer constitutive equations

The polymer stress, $\boldsymbol{\Pi}^{n+1}$, is considered constant while
solving the mass-momentum system of equations (28) and (29). Similarly, the
velocity field, $\widetilde{\mathbf{u}}^{n+1}$ is considered constant when
solving the constitutive equation (30). This allows us to numerically decouple
the mass and momentum equations from the polymer constitutive equations. To
properly account for the correct time location of $\boldsymbol{\Pi}^{n+1}$,
$\boldsymbol{\psi}^{n+1}$, $\widetilde{p}^{n+1}$ and
$\widetilde{\mathbf{u}}^{n+1}$, similar to [20], we consider $K$ inner
iterations. At $k^{th}$ iteration, first the polymer constitutive equations
are solved using velocity information from the previous,
$1^{st}<(k-1)^{th}\leq K^{th}$, step to update the values of
$\boldsymbol{\Pi}^{n+1}$ and $\boldsymbol{\psi}^{n+1}$ (at $k=1^{st}$ inner
iteration velocity information is taken from the previous time step). Then the
mass-momentum system of equations is solved using the latest polymer stress to
update the values of $\widetilde{\mathbf{u}}^{n+1}$ and $\widetilde{p}^{n+1}$.
The inner iterations are terminated when a user defined tolerance or $K$ is
reached. We obtain accurate results even with $K=1$.

###  3.2 Decoupled Schur complement approach to solve mass-momentum system

The coupled system of discrete mass and momentum equations (28) and (29) is
written in an operator form,

$\begin{bmatrix}\boldsymbol{\mathcal{L}}&\boldsymbol{\mathcal{G}}\\\
\boldsymbol{\mathcal{D}}&\boldsymbol{0}\end{bmatrix}\begin{bmatrix}\widetilde{\mathbf{u}}^{n+1}\\\
\widetilde{p}^{n+1}\end{bmatrix}=\begin{bmatrix}{\boldsymbol{m}}\\\
{0}\end{bmatrix},$ (31)

or, using a Schur complement reduction approach similar to Furuichi et al.
(2011), [41] as,

$\begin{bmatrix}\boldsymbol{\mathcal{L}}&\boldsymbol{\mathcal{G}}\\\
\boldsymbol{0}&\boldsymbol{\mathcal{S}}\end{bmatrix}\begin{bmatrix}\widetilde{\mathbf{u}}^{n+1}\\\
\widetilde{p}^{n+1}\end{bmatrix}=\begin{bmatrix}{\boldsymbol{m}}\\\
{h}\end{bmatrix},$ (32)

to decouple pressure from velocity. Here,
$\boldsymbol{\mathcal{S}}=-\boldsymbol{\mathcal{D}}\boldsymbol{\mathcal{L}}^{-1}\boldsymbol{\mathcal{G}}$
is the Schur complement of the matrix in equation (31) and
$h=-\boldsymbol{\mathcal{D}}\boldsymbol{\mathcal{L}}^{-1}{\boldsymbol{m}}$. In
these equations, $\boldsymbol{\mathcal{L}}=(\rho_{f}/\Delta
t)\boldsymbol{\delta}-\nabla^{2}$, $\boldsymbol{\mathcal{G}}=\nabla$, and,
$\boldsymbol{\mathcal{D}}=\nabla\cdot$ represent different spatial operators.
Here, $\nabla^{2}$ is a vector Laplacian operator. In the Cartesian basis, the
unit vectors are spatially constant. However, due to the spatial dependence of
unit vectors in a curvilinear basis, cross-terms exist that involve multiple
vector components in a particular component of the Laplacian of that vector.
This guides our choice of basis for vectors and tensors in section 4.2. The
discretization of these spatial operators will be described in more detail in
section 4.

${\boldsymbol{m}}=\rho_{f}\big{(}\frac{\widetilde{\mathbf{u}}^{n}}{\Delta
t}-\frac{3}{2}\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}}^{n},\mathbf{u}_{\infty}^{n},\boldsymbol{\omega}_{p}^{n};\mathbf{r})+\frac{1}{2}\widetilde{\textbf{ADV}}(\widetilde{\mathbf{u}}^{n-1},\mathbf{u}_{\infty}^{n-1},\boldsymbol{\omega}_{p}^{n-1};\mathbf{r})\big{)}+\nabla\cdot(\boldsymbol{\Pi}^{n+1}-\boldsymbol{\Pi}^{n+1}_{\infty})$
(33)

represents the sum of terms in the momentum equations that are considered
constant within a time-step (or within a $k^{th}$ iteration as described in
section 3.1), i.e., the explicit terms and the divergence of the polymer
stress. On the boundaries of the computational domain, i.e., the particle
surface and the outer boundary, the matrix-vector system of equation (31) is
appropriately changed to represent the velocity boundary conditions of
equation (16).

As highlighted in section 1, in the context of finite difference spatial
discretization, splitting methods (that advance the momentum equation in time
and solve a pressure-Poisson equation to obtain an appropriate pressure for
ensuring incompressibility) are not suitable for cases in which the fluid
inertia is small. This is related to artificial boundary conditions for
pressure (since a second-order partial differential equation now governs it)
required in the splitting methods [38, 39]. The splitting errors may be
ignored compared with more important momentum advection terms when fluid
inertia is larger than the viscous forces, making them suitable for studies of
turbulent flows. The dominance of viscous terms at lower fluid inertia values
in studies of our interest prevents us from using such methods. Fluid inertia
is quantified in the above equations by the fluid density, $\rho_{f}$.
Furthermore, due to the time advancement of the momentum equation in splitting
methods, studies with $\rho_{f}=0$ cannot be considered. $\rho_{f}=0$ studies
are useful to isolate viscoelasticity’s effect from fluid inertia completely.
This Schur complement reduction method was originally developed to study the
behavior of fluids with zero inertia and large viscosity variations, such as
in long time scale dynamics of the Earth’s convecting mantle [41]. We have
found it useful in studies of zero to moderate inertia viscoelastic fluids.

Similar to [41], we begin by solving the decoupled pressure equation

$\boldsymbol{\mathcal{S}}\widetilde{p}^{n+1}=h,$ (34)

using a Krylov subspace method that constructs a subspace
$K(\boldsymbol{\mathcal{S}},r_{0})=\text{Span}(r_{0},\boldsymbol{\mathcal{S}}r_{0},\boldsymbol{\mathcal{S}}^{2}r_{0},\cdots,\boldsymbol{\mathcal{S}}^{N-1}r_{0})$,
for $\boldsymbol{\mathcal{S}}\in\mathcal{R}^{N\times N}$, with

$\widetilde{\mathbf{u}}_{0}=\boldsymbol{\mathcal{L}}^{-1}(\boldsymbol{G}\widetilde{p}_{0}-{\boldsymbol{m}}),\hskip
14.45377ptr_{0}=h-\boldsymbol{\mathcal{S}}\widetilde{p}_{0}=\boldsymbol{\mathcal{D}}\widetilde{\mathbf{u}}_{0},$
(35)

for an initial guess $\widetilde{p}_{0}$ for $\widetilde{p}^{n+1}$. The
generalized minimum residual (GMRES) method [51] is our choice of Krylov
subspace method. The matrix vector product
$y_{i}=\boldsymbol{\mathcal{S}}x_{i}$ consists of three separate operations
defined as (similar to [41]),

$\displaystyle\mathbf{a}^{*}=\boldsymbol{\mathcal{G}}x_{i},\hskip
14.45377pt\widetilde{\mathbf{u}}^{*}_{i}=\boldsymbol{\mathcal{L}}^{-1}\mathbf{a}^{*},\hskip
14.45377pty_{i}=\boldsymbol{\mathcal{D}}\widetilde{\mathbf{u}}^{*}_{i}.$ (36)

The first and third steps are simple matrix-vector products of a known matrix
and vector that can be computed straightforwardly in an efficient manner. The
second step requires solution of a matrix equation
$\boldsymbol{\mathcal{L}}\widetilde{\mathbf{u}}*=\mathbf{a}^{*}$ and the
operator $\boldsymbol{\mathcal{L}}$ involves the sum of an identity and a
Laplacian operator. Therefore, we use the aggregation-based algebraic
multigrid (AGMG) method of [52, 53] to solve this elliptic equation
efficiently. We use the implementation provided in [54]. Once the GMRES method
terminates upon reaching a sufficient user-defined tolerance (see [51, 55] for
details), say in $M$ iterations, the solution for $\widetilde{p}^{n+1}$ is
formed as,

$\widetilde{p}^{n+1}=\widetilde{p}_{0}+\Sigma_{i=1}^{M}y_{i}/||y_{i}||_{2}l_{i},$
(37)

where for each GMRES iteration $i$, $y_{i}$ is a vector defined by equation
(36) and $l_{i}$ is a scalar defined within the GMRES procedure (see [51, 55]
for details). If we keep track of different $\widetilde{\mathbf{u}}^{*}_{i}$
in the GMRES iteration for solution of the pressure equation (34), by a simple
and fast operation we can construct,

$\delta\widetilde{\mathbf{u}}^{n+1}=\Sigma_{i=1}^{M}\widetilde{\mathbf{u}}^{*}_{i}l_{i}.$
(38)

The velocity field at step n+1 is thus obtained by

$\widetilde{\mathbf{u}}^{n+1}=\widetilde{\mathbf{u}}_{0}-\delta\widetilde{\mathbf{u}}^{n+1},$
(39)

where $\widetilde{\mathbf{u}}_{0}$ is already calculated and defined in
equation (35). Our methodology to solve the coupled system of mass and
momentum equation closely follows that of [41] for solving the decoupled
pressure equation (34). However, careful observation of the GMRES method
allows us to obtain the solution for velocity, $\widetilde{\mathbf{u}}^{n+1}$
as an auxiliary product of the same calculation. This avoids the need to solve
a velocity equation
$\widetilde{\mathbf{u}}^{n+1}=\boldsymbol{\mathcal{L}}^{-1}({\boldsymbol{m}}-\boldsymbol{\mathcal{G}}p^{n+1})$
and saves CPU time. In [41] authors note a similar point but still solve the
velocity equation after obtaining the pressure solution. This has so far
proven to be unnecessary for our studies.

###  3.3 Ensuring stretch limited by maximum polymer extensibility in FENE
models

In dumbbell based models for polymer configuration, as discussed in section 2,
the polymer configuration is
$\boldsymbol{\Lambda}=\langle\mathbf{qq}\rangle_{\text{polymer
configuration}}$, where $\mathbf{q}$ is end-to-end vector (non-dimensionalized
with polymer’s radius of gyration) of the dumbbell and the angle brackets
represent average over polymer configuration. Therefore, the polymer stretch
is $\sqrt{\text{tr}(\boldsymbol{\Lambda})}$. In FENE models such as the FENE-P
and FENE-CR (equation (19)) models, the polymers have a maximum extensibility
$L$. Following the technique introduced in [20] we numerically impose this
constraint by separately evolving the variable
$\gamma=1/f=1-\text{tr}(\boldsymbol{\Lambda})/L^{2}$ (equation (20)) used in
the FENE models. An evolution equation for $\gamma$ is obtained by taking the
trace of the polymer constitutive equation (17) with one of the FENE models
from equation (19),

$\frac{\partial\gamma}{\partial
t}+\mathbf{u}\cdot\nabla\gamma+\frac{2}{L^{2}}\boldsymbol{\Lambda}:\nabla\mathbf{u}+\frac{1}{De}(\frac{\gamma-1}{\gamma}+\beta)=0,\hskip
14.45377pt\beta=\begin{cases}\frac{3}{L^{2}-3},&\text{FENE-P}\\\
\frac{3}{L^{2}\gamma},&\text{FENE-CR}\end{cases}.$ (40)

Similar to [20], we temporally discretize equation (40) using a Cranck-
Nicholson scheme for the relaxation terms
(${1}/{De}((\gamma-1)/{\gamma}+\beta)$) and treating the advection
($\mathbf{u}\cdot\nabla\gamma$) and stretching
(${2}/{L^{2}}\boldsymbol{\Lambda}:\nabla\mathbf{u}$) terms explicitly. This
leads to the following quadratic equation for $\gamma^{n+1}$,

$(\gamma^{n+1})^{2}+\gamma^{n+1}\Delta t\Big{(}-\frac{\gamma^{n}}{\Delta
t}+\mathbf{u}^{n}\cdot\nabla\gamma^{n}+\frac{2\boldsymbol{\Lambda^{n}}:\nabla\mathbf{u}^{n}}{L^{2}}+\frac{\gamma^{n}-0.5}{De\gamma^{n}}+\frac{1}{DeL^{2}}\alpha_{1}^{n}\Big{)}-\frac{\Delta
t}{2De}\alpha_{2}=0,$ (41)

where $\alpha_{1}^{n}=3/(L^{2}-3)$ for FENE-P and $1.5/\gamma^{n}$ for FENE-
CR, and, $\alpha_{2}=1$ for FENE-P and $1-3/L^{2}$ for FENE-CR. This equation
leads to two real roots with opposite signs [20]. The negative root is
unphysical since it implies $\text{tr}(\boldsymbol{\Lambda}^{n+1})>L^{2}$.
Choosing the positive root ensures that the polymer stretch is upper-bounded
by the maximum extensibility, $L$. At each time step, $n+1$, we calculate
$f^{n+1}=1/\gamma^{n+1}$ from this treatment and use this value of $f^{n+1}$
in the relaxation term of the discretized polymer constitutive equation (30).

###  3.4 Evolution of particle orientation and velocity boundary conditions

As mentioned in section 2, we solve the governing equations in a particle
fixed reference frame. In our simulations, the inertial (or laboratory) frame
is defined either from the initial particle orientation (e.g., a particle
rotating about its axis in quiescent fluid) or through the geometry of the
imposed flow (e.g., a reference frame fixed with the imposed simple shear
flow, uniaxial extensional flow or a uniform flow field). The particle
orientation is defined using quaternions,
$\mathbf{q}=\begin{bmatrix}q_{1}&q_{2}&q_{3}&q_{4}\end{bmatrix}^{T}$, (see
chapter 8 of [56]) that are related to the Euler angles, $\theta,\phi,\psi$,
between the particle-fixed and inertial reference frames,

$\displaystyle\begin{split}q_{1}=&\sin(\theta/2)\cos((\phi-\psi)/2),\hskip
14.45377ptq_{2}=\sin(\theta/2)\sin((\phi-\psi)/2),\\\
q_{3}=&\cos(\theta/2)\sin((\phi+\psi)/2),\hskip
14.45377ptq_{4}=\cos(\theta/2)\cos((\phi+\psi)/2).\end{split}$ (42)

The sequence of three rotations that define the Euler angles are described in
[56]: $\theta=\sin^{-1}(-X_{3})$, $\phi=\sin^{-1}(Y_{3}/\sqrt{1-X_{3}^{2}})$,
$\psi=\sin^{-1}(X_{2}/\sqrt{1-X_{3}^{2}})$. Here, $X_{2}$ and $X_{3}$ are the
projections of the $X$ (or 1) axis of the particle-fixed reference frame on
the $2$ and $3$ axis, respectively, of the inertial reference frame. $Y_{2}$
and $Y_{3}$ are the same projections of the $Y$ (or 2) axis of the particle-
fixed reference frame. The quaternion formulation has previously been used to
study fluid flows around prolate spheroids in [26, 57]. The evolution equation
for quaternions is related to the particle’s angular velocity,
$\boldsymbol{\omega}_{p}=\begin{bmatrix}\omega_{1}&\omega_{2}&\omega_{3}\end{bmatrix}^{T}$,

$\frac{d}{dt}\begin{bmatrix}q_{1}\\\ q_{2}\\\ q_{3}\\\
q_{4}\end{bmatrix}=\frac{1}{2}\begin{bmatrix}q_{4}&-q_{3}&q_{2}&q_{1}\\\
q_{3}&q_{4}&-q_{1}&q_{2}\\\ -q_{2}&q_{1}&q_{4}&q_{3}\\\
-q_{1}&-q_{2}&-q_{3}&q_{4}\end{bmatrix}\begin{bmatrix}\omega_{1}\\\
\omega_{2}\\\ \omega_{3}\\\ 0\end{bmatrix}.$ (43)

We consider second and third-order accurate Adams-Bashforth schemes to
discretize equation (43) temporally. The transformation matrix from the
inertial frame to the particle-fixed frame is,

$\boldsymbol{A}(t)=2\begin{bmatrix}q_{1}^{2}+q_{4}^{2}-1/2&q_{1}q_{2}+q_{3}q_{4}&q_{1}q_{3}-q_{2}q_{4}\\\
q_{1}q_{2}-q_{3}q_{4}&q_{2}^{2}+q_{4}^{2}-1/2&q_{2}q_{3}+q_{1}q_{4}\\\
q_{1}q_{3}+q_{2}q_{4}&q_{2}q_{3}-q_{1}q_{4}&q_{3}^{2}+q_{4}^{2}-1/2\end{bmatrix}.$
(44)

$\boldsymbol{A}(t)$ is an orthonormal matrix. A vector in the particle-fixed
frame, $\mathbf{b}_{\text{particle}}$ is transformed into the inertial
reference frame with
$\mathbf{b}_{\text{inertial}}=\boldsymbol{A}(t)^{T}\cdot\mathbf{b}_{\text{particle}}$.
The major axis or the center-line of the spheroidal particle is along the
$x_{3}$ axis of the particle-fixed frame. At a particular time, in the
inertial reference frame the particle orientation vector i.e. the particle
center line is
$\mathbf{p}=\begin{bmatrix}q_{1}q_{3}+q_{2}q_{4}&q_{2}q_{3}-q_{1}q_{4}&q_{3}^{2}+q_{4}^{2}-1/2\end{bmatrix}$.
At the end of each time step, after updating the quaternions we update the
velocity boundary conditions. Consider a general case where the imposed fluid
motion consists of a uniform velocity, $\mathbf{u}_{far}^{\text{inertial}}$
and a velocity gradient, $\boldsymbol{\Gamma}_{\text{far}}^{\text{inertial}}$
in the laboratory frame. The particle’s translation and angular velocities
(numerical methods to evaluate these are discussed in section 3.5) in the
frame aligned with the particle are $\mathbf{u}_{p}$ and
$\boldsymbol{\omega}_{p}$ respectively. In the frame of reference aligned with
the particle coordinates, the appropriately rotated imposed fluid motion at
the outer boundary is,

$\mathbf{u}_{\infty}(\mathbf{r},t)=\boldsymbol{A}(t)\cdot\mathbf{u}_{\text{far}}^{\text{inertial}}+\mathbf{r}\cdot\boldsymbol{A}(t)\cdot\boldsymbol{\Gamma}_{\text{far}}^{\text{inertial}}\cdot\boldsymbol{A}(t)^{T}-\mathbf{u}_{p}-\boldsymbol{\omega}_{p}\times\mathbf{r},\hskip
14.45377pt\mathbf{r}=\mathbf{r}_{\infty}.$ (45)

###  3.5 Particle velocities

As mentioned in section 2 the particle’s angular and translation velocities,
$\boldsymbol{\omega}_{p}$ and $\mathbf{u}_{p}$ may be a prescribed quantity in
a study of interest, so that the procedure mentioned above has all the
required quantities to solve the equations. We are also interested in studying
the scenarios where the particle is free to move due to the imposed and fluid-
generated forces and torques. The case when particle inertia is zero must be
dealt with differently than the one with finite particle inertia, and we will
consider these two cases next.

####  3.5.1 Finite particle inertia

At each time step the particle motion must satisfy the Newton’s equations (6)
and (7). These equations are discretized using a second order implicit Crank-
Nicholson scheme,

$\displaystyle\mathbf{u_{p}}^{n+1}=$
$\displaystyle\mathbf{u_{p}}^{n}+\frac{\Delta
t}{2\rho_{p}V_{p}}(\mathbf{f}^{n+1}_{\text{net}}+\mathbf{f}^{n}_{\text{net}})$
(46) $\displaystyle\boldsymbol{\omega}_{p}^{n+1}=$
$\displaystyle\boldsymbol{\omega}_{p}^{n}+\frac{5\Delta
t}{2\rho_{p}V_{p}r_{\text{minor}}^{2}}\begin{bmatrix}\frac{1}{1+\kappa^{2}}&0&0\\\
0&\frac{1}{1+\kappa^{2}}&0\\\
0&0&0.5\end{bmatrix}\cdot(\mathbf{q}_{\text{net}}^{n+1}+\mathbf{q}_{\text{net}}^{n}),$
(47)

where $V_{p}=4\pi\kappa r_{minor}^{3}/3$,

$\displaystyle\begin{split}\mathbf{f}^{n+1}_{\text{net}}=&{\mathbf{f}_{\text{fluid}}^{n+1}+\mathbf{f}_{\text{ext.}}^{n+1},\hskip
14.45377pt\mathbf{f}_{\text{fluid}}^{n+1}=\mathbf{f}(\boldsymbol{\tau})^{n+1}+\mathbf{f}(\boldsymbol{\Pi})^{n+1}}\\\
\mathbf{q}_{\text{net}}^{n+1}=&\mathbf{q}_{\text{fluid}}^{n+1}+\mathbf{q}_{\text{ext.}}^{n+1},\hskip
14.45377pt\mathbf{q}_{\text{fluid}}^{n+1}=\mathbf{q}(\boldsymbol{\tau})^{n+1}+\mathbf{q}(\boldsymbol{\Pi})^{n+1}.\end{split}$
(48)

are the net force and torque acting on the particle at time step $n+1$.
$\mathbf{f}_{\text{fluid}}^{n+1}$ and $\mathbf{q}_{\text{fluid}}^{n+1}$ are
evaluated from equations (9) and (10) once the fluid stress on the particle
surface at time step $n+1$,
$\boldsymbol{\sigma}^{n+1}=\boldsymbol{\tau}^{n+1}+\boldsymbol{\Pi}^{n+1}$, is
available. $\mathbf{f}_{\text{ext.}}^{n+1}$ and
$\mathbf{q}_{\text{ext.}}^{n+1}$ are the externally imposed (possibly time-
varying) force and torques in the frame aligned with the particle at time step
$n+1$. We use the first order explicit Euler scheme in the first time step.

####  3.5.2 Zero particle inertia

As mentioned earlier in section 1 and 3.2, we are also interested in the
scenario of zero particle inertia or zero $\rho_{p}$. In our studies, such a
case arises when we want to completely remove the inertial effects and study
the influence of viscoelasticity, but it could also be used to investigate
massless particles in a fluid with finite inertia. In both of these cases,
Newton’s equations governing the particle motion at each time step are,

$\displaystyle\begin{split}\mathbf{f}^{n+1}_{\text{net}}=&{\mathbf{f}_{\text{fluid}}^{n+1}+\mathbf{f}_{\text{ext.}}^{n+1}}=0,\hskip
14.45377pt\mathbf{f}_{\text{fluid}}^{n+1}=\mathbf{f}(\boldsymbol{\tau})^{n+1}+\mathbf{f}(\boldsymbol{\Pi})^{n+1},\\\
\mathbf{q}_{\text{net}}^{n+1}=&\mathbf{q}_{\text{fluid}}^{n+1}+\mathbf{q}_{\text{ext.}}^{n+1}=0,\hskip
14.45377pt\mathbf{q}_{\text{fluid}}^{n+1}=\mathbf{q}(\boldsymbol{\tau})^{n+1}+\mathbf{q}(\boldsymbol{\Pi})^{n+1},\end{split}$
(49)

and the time marching used for finite $\rho_{p}$ in (48) cannot be employed.
The governing equations (49) are viewed as force and torque constraints that
velocity, pressure and polymer stress field must satisfy at each time step to
yield an appropriate $\mathbf{u_{p}}^{n+1}$ and
$\boldsymbol{\omega}_{p}^{n+1}$. Padhy et al. (2013) [58] used a secant method
to iteratively impose the torque-free constraint on a sphere rotating in a
cross shear flow. We first consider this method for imposing force- and
torque-free constraints for a finite fluid inertia case. Subsequently, we show
that a novel decomposition of inertia-less fluid’s momentum equation can be
used to impose these constraints in a non-iterative and hence computationally
efficient manner. These two techniques are discussed next.

Secant iteration method for finite fluid inertia and zero particle inertia: At
each time step or within each inner $k$ iteration (section 3.1) the polymer
stress, $\boldsymbol{\Pi}$ and hence the polymer torque and force,
$\mathbf{q}(\boldsymbol{\Pi})^{n+1}$ and $\mathbf{f}(\boldsymbol{\Pi})^{n+1}$,
from equations equations (9) and (10), are fixed. $\mathbf{u_{p}}^{n+1}$ and
$\boldsymbol{\omega}_{p}^{n+1}$ are iterated along with
$\widetilde{\mathbf{u}}^{n+1}$ and $\widetilde{p}^{n+1}$ to generate the
appropriate Newtonian solvent torque and force,
$\mathbf{q}(\boldsymbol{\tau})^{n+1}$ and
$\mathbf{f}(\boldsymbol{\tau})^{n+1}$, that ensures the torque and force
balance in equation (49). In practice, once the polymer constitutive equation
(30) is solved, we first obtain the polymeric torque and force,
$\mathbf{q}(\boldsymbol{\Pi})^{n+1}$ and $\mathbf{f}(\boldsymbol{\Pi})^{n+1}$,
from equation (9). The mass-momentum system described by equation (32) depends
on the particle’s angular velocity $\boldsymbol{\omega}_{p}$ (see equations
(15), (16), and (33)). In the secant iterations method, iterations proceed by
solving the mass-momentum system by the Schur complement method described in
section 3.2 and obtaining the Newtonian torque and force,
$\mathbf{q}(\boldsymbol{\tau})^{n+1}$ and
$\mathbf{f}(\boldsymbol{\tau})^{n+1}$, from equation (9). After each secant
iteration, $s$, the particle’s angular and translational velocities
$\boldsymbol{\omega}_{p}$ and $\mathbf{u}_{p}$ are updated component-wise,

$\omega_{p,i}^{s+1}=\omega_{p,i}^{s}-q_{net,i}^{s}\frac{\omega_{p,i}^{s}-\omega_{p,i}^{s-1}}{q_{net,i}^{s}-q_{net,i}^{s-1}},\hskip
14.45377ptu_{p,i}^{s+1}=u_{p,i}^{s}-f_{net,i}^{s}\frac{u_{p,i}^{s}-u_{p,i}^{s-1}}{f_{net,i}^{s}-f_{net,i}^{s-1}},\hskip
14.45377pti\in[1,3].$ (50)

$\boldsymbol{\omega}_{p}^{s+1}$ and $\mathbf{u}_{p}^{s+1}$ are used in the
velocity boundary condition (equation (45)) for the next secant iteration. The
secant iterations are stopped once the magnitudes of all components of the net
torque and force on the particle, $\mathbf{q}_{\text{net}}^{n+1}$ and
$\mathbf{f}_{\text{net}}^{n+1}$, are below prescribed tolerances. The
numerical solution at this point consists of the velocity, pressure, and
polymer stress fields and the particle’s required angular and translational
velocity that satisfy the force and torque constraints of equation (49).

Novel resistivity formulation for zero fluid and particle inertia in non-
Newtonian fluids: The system of mass and momentum equations and the torque-
and force- free constraint is quasi-steady when both fluid and particle
inertia are neglected. All the variables from here until the end of this
section are taken at time step $n+1$ and we will supress the superscript for
clarity. Upon neglecting fluid inertia, the momentum equation (2) becomes

$\nabla\cdot\boldsymbol{\sigma}=0,$ (51)

subject to the boundary conditions mentioned in equation (5) and the stress
tensor, $\boldsymbol{\sigma}$, in equation (4). In this case, the momentum and
mass (1) conservation equations are linear in the velocity and pressure.
Therefore, using

$\mathbf{u}=\mathbf{u}^{\text{M}}+\mathbf{u}^{\text{P}},\hskip
14.45377ptp=p^{\text{M}}+p^{\text{P}},$ (52)

the system of mass and momentum equations along with the associated boundary
conditions is linearly decomposed into two parts. $\mathbf{u}^{\text{M}}$ and
$p^{\text{M}}$ represent the motion-induced velocity and pressure fields and
$\mathbf{u}^{\text{P}}$ and $p^{\text{P}}$ the polymer-induced fields. The
equations governing the motion-induced part are,

$\nabla\cdot\mathbf{u}^{\text{M}}=0,\hskip
14.45377pt\nabla\cdot\boldsymbol{\sigma}^{\text{M}}=0,$ (53)

with,

$\boldsymbol{\sigma}^{\text{M}}=\boldsymbol{\tau}^{\text{M}}=-p^{\text{M}}\boldsymbol{\delta}+\mu(\nabla\mathbf{u}^{\text{M}}+(\nabla\mathbf{u}^{\text{M}})^{T})$
(54)

and the boundary conditions,

$\displaystyle\begin{split}\mathbf{u}^{\text{M}}=\mathbf{0},\text{ on particle
surface },\hskip
14.45377pt\mathbf{u}^{\text{M}}=\mathbf{u}_{\infty}(\mathbf{r})=\mathbf{r}\cdot\boldsymbol{\Gamma}+\mathbf{u}_{0}-\mathbf{u}_{p}-\boldsymbol{\omega}_{p}\times\mathbf{r},\text{
as }|\mathbf{r}|\rightarrow\infty\approx\mathbf{r}_{\infty}.\end{split}$ (55)

The polymer-induced part is governed by,

$\displaystyle\nabla\cdot\mathbf{u}^{\text{P}}=0,\hskip
14.45377pt\nabla\cdot\boldsymbol{\sigma}^{\text{P}}=0,$ (56)

with,

$\boldsymbol{\sigma}^{\text{P}}=\boldsymbol{\tau}^{\text{P}}+\boldsymbol{\Pi}=-p^{\text{P}}\boldsymbol{\delta}+\mu(\nabla\mathbf{u}^{\text{P}}+(\nabla\mathbf{u}^{\text{P}})^{T}+\boldsymbol{\Pi})$
(57)

and the boundary conditions,

$\mathbf{u}^{\text{P}}=\mathbf{0},\text{ on particle surface },\hskip
14.45377pt\mathbf{u}^{\text{P}}=\mathbf{0},\text{ as
}|\mathbf{r}|\rightarrow\infty\approx\mathbf{r}_{\infty}.$ (58)

The hydrodynamic force and torque acting on the particle are also decomposed
into a motion- and polymer-induced part,

$\mathbf{f}_{\text{fluid}}=\mathbf{f}_{\text{fluid}}^{\text{M}}+\mathbf{f}_{\text{fluid}}^{\text{P}},\hskip
7.22743pt\mathbf{q}_{\text{fluid}}=\mathbf{q}_{\text{fluid}}^{\text{M}}+\mathbf{q}_{\text{fluid}}^{\text{P}},$
(59)

with,

$\mathbf{f}_{\text{fluid}}^{\text{M}}=\mathbf{f}(\boldsymbol{\sigma}^{\text{M}}),\hskip
3.61371pt\mathbf{f}_{\text{fluid}}^{\text{P}}=\mathbf{f}(\boldsymbol{\sigma}^{\text{P}}),\hskip
3.61371pt\mathbf{q}_{\text{fluid}}^{\text{M}}=\mathbf{q}(\boldsymbol{\sigma}^{\text{M}}),\hskip
3.61371pt\text{and},\hskip
3.61371pt\mathbf{q}_{\text{fluid}}^{\text{P}}=\mathbf{q}(\boldsymbol{\sigma}^{\text{P}}),$
(60)

where, $\mathbf{f}(\boldsymbol{\sigma})$ and $\mathbf{q}(\boldsymbol{\sigma})$
are defined in equation (9). This way of decomposing the momentum equation
attributes the entire effect of particle motion (via the boundary conditions
in equation (55)) to the motion-induced part of the governing equations. This
part has no influence of the polymer stress. The polymer-induced pressure and
velocity fields are forced by the effect of the polymer stress, but are not
explicitly affected by the particle’s motion. A key observation allowing us to
circumvent an iterative procedure in calculating the particle’s motion is that
the motion-induced equations represent a Stokes flow (inertia-less Newtonian
flow). Hence, the motion induced force and torque are simply,

$\displaystyle\begin{split}\mathbf{f}_{\text{fluid}}^{M}=\boldsymbol{F}_{p}\cdot\mathbf{u}_{p}+\mathbf{f}_{\infty},\hskip
3.61371pt\mathbf{q}_{\text{fluid}}^{M}=\boldsymbol{Q}_{p}\cdot\boldsymbol{\omega}_{p}+\mathbf{q}_{\infty}.\end{split}$
(61)

The tensors $\boldsymbol{F}_{p}$ and $\boldsymbol{Q}_{p}$ depend upon only on
the particle shape, and the vectors $\mathbf{f}_{\infty}$ and
$\mathbf{q}_{\infty}$ depend upon both the particle shape and the imposed
flow. These are either evaluated analytically for simple particle shapes and
imposed flows or calculated by considering only the motion-induced mass and
momentum equations with appropriate boundary conditions. For example, the 21
component of $\boldsymbol{F}_{p}$ is simply the second component of
$\mathbf{f}(\boldsymbol{\sigma}^{\text{M}})$ for a boundary condition
$\mathbf{u}^{\text{M}}=\mathbf{0},\text{ on particle surface and
}\mathbf{u}^{\text{M}}=[1\hskip 3.61371pt0\hskip 3.61371pt0]^{T},\text{ as
}|\mathbf{r}|\rightarrow\infty\approx\mathbf{r}_{\infty}$.
$\mathbf{f}_{\infty}$ and $\mathbf{q}_{\infty}$ are simply the hydrodynamic
force and torque on a fixed particle in the imposed flow
($\mathbf{r}\cdot\boldsymbol{\Gamma}+\mathbf{u}_{0}$) of an inertia-less
Newtonian fluid. From the force- and torque-free constraints, we obtain

$\mathbf{u}_{p}=-\boldsymbol{F}_{p}^{-1}\cdot(\mathbf{f}_{\infty}+\mathbf{f}_{\text{fluid}}^{\text{P}}+\mathbf{f}_{\text{ext}}),\text{
and,
}\boldsymbol{\omega}_{p}=-\boldsymbol{Q}_{p}^{-1}\cdot(\mathbf{q}_{\infty}+\mathbf{q}_{\text{fluid}}^{\text{P}}+\mathbf{q}_{\text{ext}})$
(62)

in a non-iterative way. This method of obtaining a particle’s motion has long
been used in micro-hydrodynamics of Newtonian fluids and is called a
resistivity formulation [59]. However, with the novel decomposition of the
inertia-less momentum equation we have shown that a similar principal can be
applied in a computationally useful manner for non-Newtonian problems.

Unlike the polymer induced pressure and velocity fields, the determination of
the motion induced velocity and pressure fields does not require application
of the Schur complement approach described in section 3.2 at each time step of
the simulation. Once the particle’s translational and angular velocity are
available, these fields are either taken from the analytical solution of
Stokes flow around the particle or from a linear superposition of pre-
calculated numerical solutions of a few fundamental incompressible flows
around the particle. The boundary conditions at a given time instant can be
represented by a superposition of the boundary conditions in these fundamental
flows. In the most general case when the velocity at the outer boundary in
equation (55) has non-zero components for all the components of the effective
(incompressible) imposed velocity gradient
($\boldsymbol{\Gamma+\epsilon\cdot\omega}_{p}$) and the effective imposed
velocity ($\mathbf{u}_{0}-\mathbf{u}_{p}$), 11 fundamental Stokes flows have
to be pre-calculated. While we are primarily considering linear velocity
boundary conditions, this method can be extended to other types of boundary
conditions where the Stokes flows required to be pre-calculated will be
different. Using this resistivity formulation, our simulations become
significantly faster relative to the secant iteration method.

The two most computationally intensive components of our method are the
weighted Jacobi method to iteratively solve equation (30) for
$\boldsymbol{\Psi}^{n+1}$ and the algebraic multigrid method to invert
$\boldsymbol{\mathcal{L}}$. The domain decomposition method implemented in
Message Passing Interface (MPI) is used to parallelize the code in the
$\xi_{1}$ and $\xi_{2}$ directions (the $\xi_{3}$ direction could also be
parallelized in future). The weighted Jacobi method is implemented completely
in-house and shows good strong and weak scaling for 100s of processing units.
The parallel scaling of the algebraic multigrid method is dependent on the
external vendor’s capability. We use the academic version of AGMG [54] for
which good scaling was only obtained up to about 20 processing units. The
examples considered in section 5 are completed within a reasonable time of up
to four days using 20 processing units. However, for more stringent cases,
AGMG can be replaced with other open source algebraic multigrid solvers
available in libraries such as the BoomerAMG of HYPRE [60] that scales well up
to a large number of processing units, but requires more parameter tuning.

In this section 3, we have described the algorithmic structure of our
numerical method in a coordinate system free manner. In the next section, we
will delve into the specific choice of the prolate spheroidal coordinate
system and the finite difference schemes used to discretize the spatial
gradients in the governing equations. This will allow us to discretize the
computational space and various operators mentioned above.

##  4 Spatial discretization

One of the benefits of using body-fitted grids and expressing equations in the
particle reference frame is that the discrete space or mesh remains fixed as
the particle rotates, and the discretized spatial operators defined here need
to be calculated only once before the equations are evolved in time.

###  4.1 Prolate spheroidal coordinate system for spatial discretization of
computational domain

We use a prolate spheroidal coordinate system to spatially discretize the
particle surface and the fluid region around it. The transformation between
the prolate spheroidal ($\boldsymbol{\xi}$) and cartesian ($\mathbf{x}$)
coordinates, defined for the particle-fixed reference frame, are

$x_{1}=f\sinh(\xi_{1})\sin(\xi_{2})\sin(\xi_{3}),\hskip
14.45377ptx_{2}=f\sinh(\xi_{1})\sin(\xi_{2})\cos(\xi_{3}),\hskip
14.45377ptx_{3}=f\cosh(\xi_{1})\cos(\xi_{2}),$ (63)

where $\xi_{2}\in[0,\pi]$ and $\xi_{3}\in[0,2\pi]$. $f$ is the focal length of
the prolate spheroidal particle. Here, $\xi_{1}$, $\xi_{2}$, and $\xi_{3}$ are
similar to the radial, polar and azimuthal directions respectively in a
spherical coordinate system. The surface of the particle with an aspect ratio,
$\kappa$, is exactly modeled as one of the coordinate surface i.e.
$\xi_{1}=\xi_{1}^{\text{surface}}$. We consider a prolate spheroid spheroid
with minor radius, $r_{\text{minor}}$. The surface, $\xi_{1}^{\text{surface}}$
and focal length, $f$, are,

$\xi_{1}^{\text{surface}}=\frac{1}{2}\log\Big{(}\frac{\kappa+1}{\kappa-1}\Big{)},\hskip
14.45377ptf=\frac{r_{\text{minor}}}{\sinh(\xi_{1}^{\text{surface}})}.$ (64)

The outer surface of the computational domain, $\mathbf{r}_{\infty}$, is also
a prolate spheroidal surface that has a constant $\xi_{1}=\xi_{1}^{\infty}$,

$\xi_{1}^{\infty}=\sinh^{-1}(||\mathbf{r}^{\text{minor}}_{\infty}||_{2}/f),$
(65)

where $||\mathbf{r}^{\text{minor}}_{\infty}||_{2}$ (2-norm or the Euclidean
norm of $\mathbf{r}^{\text{minor}}$) is the minor radius of the outer surface.
It is a user prescribed parameter. The Euclidean distance of any point on a
prolate spheroid $||\mathbf{r}_{\infty}||_{2}$ is related to its minor axis
through
$||\mathbf{r}_{\infty}||_{2}^{2}=||\mathbf{r}^{\text{minor}}_{\infty}||_{2}^{2}+f^{2}\cos^{2}(\xi_{2})$.
If this surface is placed far from the particle, such as the outer surface,
$||\mathbf{r}_{\infty}||\approx||\mathbf{r}^{\text{minor}}_{\infty}||_{2}$. In
other words, the outer surface is a nearly spherical surface. This is shown
schematically in figure 1 and for an actual discretized example in the left
panel of figure 2. The computational domain is defined within the limits:
$\xi_{1}\in[\xi_{1}^{\text{surface}},\xi_{1}^{\infty}]$, $\xi_{2}\in[0,\pi]$,
$\xi_{3}\in[0,2\pi]$.

In order to prevent pressure aliasing that leads to the checkerboard effect
causing spurious pressure oscillations, we use a staggered grid arrangement
[61]. The three velocity components are stored at the same location that we
term velocity grid. Pressure and the components of the polymer configuration
tensor $\boldsymbol{\Lambda}$ and the latter’s matrix logarithm
$\boldsymbol{\Psi}$ are stored at a location that is staggered relative to the
velocity grid. These staggered locations form the pressure grid. The velocity
and pressure grids are

$\displaystyle\text{Velocity Grid:
}\xi_{1}^{\text{dist.,vel}}\cup\xi_{2}^{\text{dist.,vel}}\cup\xi_{3}^{\text{dist.,vel}},\hskip
14.45377pt\text{Pressure Grid:
}\xi_{1}^{\text{dist.,pres}}\cup\xi_{2}^{\text{dist.,pres}}\cup\xi_{3}^{\text{dist.,pres}},$
(66)

where

$\displaystyle\begin{split}\xi_{1,i}^{\text{dist.,vel}}=&\xi_{1}^{\text{surface}}+\frac{\xi_{1}^{\infty}-\xi_{1}^{\text{surface}}}{N_{1}-1}(i-1)[c_{1}\frac{(i-1)(i-0.5)}{(N_{1}-1)^{2}}+1],\hskip
7.22743pti=[1,2,\cdots,N_{1}],\\\
\xi_{1,1}^{\text{dist.,pres}}=&\xi_{1}^{\text{surface}},\hskip
7.22743pt\xi_{1,i}^{\text{dist.,pres}}=\frac{\xi_{1,i}^{\text{dist.,vel}}+\xi_{1,i-1}^{\text{dist.,vel}}}{2},\hskip
7.22743pti=[2,3,\cdots,N_{1}],\hskip
7.22743pt\xi_{1,N_{1}+1}^{\text{dist.,pres}}=\xi_{1}^{\infty},\\\
\xi_{2,j}^{\text{dist.,pres}}=&\pi\frac{j-1}{N_{2}}[c_{2}\frac{(j-1)(j-0.5)}{(N_{2}-1)^{2}}+1],\hskip
7.22743ptj=[2,3,\cdots,N_{2}+1],\\\
\xi_{2,j}^{\text{dist.,vel}}=&\frac{\xi_{2,j}^{\text{dist.,pres}}+\xi_{2,j+1}^{\text{dist.,pres}}}{2},\hskip
7.22743ptj=[1,2,\cdots,N_{2}],\\\
\xi_{3,k}^{\text{dist.,vel}}=&2\pi\frac{k-1}{N_{3}-1},\hskip
7.22743ptk=[1,2,\cdots,N_{3}],\hskip
7.22743pt\xi_{3,k}^{\text{dist.,pres}}=\frac{\xi_{3,k}^{\text{dist.,vel}}+\xi_{3,k+1}^{\text{dist.,vel}}}{2},\hskip
7.22743ptk=[1,2,\cdots,N_{3}-1],\\\ \end{split}$ (67)

when $c_{1}$ and $c_{2}$ are zero, we obtain a uniform grid in spheroidal
coordinates. A uniform spheroidal grid is naturally more clustered near the
particle surface in the Euclidean sense than a uniform Cartesian grid. When
$c_{1}<0$, further grid clustering towards the particle and the outer surface
is obtained. Another useful grid for this purpose to obtain grid clustering
near the particle surface only (and appropriate when $c_{1}\neq 0$) is
$\xi_{1,i}^{\text{dist.,vel}}=\xi_{1}^{\text{surface}}+({\xi_{1}^{\infty}-\xi_{1}^{\text{surface}}})\\{\exp[-c_{1}(i-1)/(N_{1}-1)]-1\\}/\\{\exp(-c_{1})-1\\},\hskip
7.22743pti=[1,2,\cdots,N_{1}]$. Setting $c_{2}<0$ allows us to cluster the
grids (in $\xi_{2}$ coordinate) towards the major axis of the particle (which
is useful in studies such as extensional flow around the particle with its
major axis aligned with the extensional axis). We only consider a uniform grid
in the azimuthal, $\xi_{3}$ direction as the particle is axisymmetric, and
there is no a priori preference of gradients at a particular $\xi_{3}$
location for a general linear flow around the particle. Any other function to
describe a non-uniform grid can be used with our method as the finite
difference and interpolation schemes described in section 4.3 do not assume a
grid type. Figure 2 shows the pressure grid for a spheroid with $\kappa=4$,
$||\mathbf{r}_{\infty}||_{2}\approx||\mathbf{r}^{\text{minor}}_{\infty}||_{2}=80$,
$N_{1}=90$, $N_{2}=N_{3}=71$ and $c_{1}=c_{2}=0$.

(a)

(b)

Fig. 2: Discretized computational domain ($\kappa=4$, $r_{\text{minor}}=1$,
$||\mathbf{r}_{\infty}||_{2}\approx||\mathbf{r}^{\text{minor}}_{\infty}||_{2}=80$,
$N_{1}=90$, $N_{2}=N_{3}=71$ and $c_{1}=c_{2}=0$): left panel shows the full
domain, and right panel shows a zoomed region near particle surface (red).
Pressure grid is shown here, velocity grid is staggered relative to this
according to definitions in equation (67).

###  4.2 Caretsian basis for vectors and tensors and discrete representation
of spatial operators

We use the Cartesian basis on this prolate spheroidal grid to represent the
velocity vector or the polymer configuration tensor at each location. This is
done to prevent the coupling of different velocity vector components in the
vector Laplacian operator (also see equation (31) and following discussion)
appearing in the momentum equations. The momentum equation written in
Cartesian basis on the prolate spheroidal grid is,

${\rho_{f}}\big{(}\frac{\partial\widetilde{u}_{i}}{\partial
t}+\widetilde{\text{ADV}}(\widetilde{\mathbf{u}},\mathbf{u}_{\infty},\boldsymbol{\omega}_{p};\mathbf{r})_{i}\big{)}=-\frac{\partial\xi_{k}}{\partial
x_{i}}\frac{\partial\widetilde{p}}{\partial\xi_{k}}+\nabla^{2}_{\text{SphinCart}}\widetilde{u}_{i}+\frac{\partial\xi_{k}}{\partial
x_{j}}\frac{\partial(\Pi_{ji}-\Pi_{ji,\infty})}{\partial\xi_{k}},$ (68)

where,

$\nabla^{2}_{\text{SphinCart}}=\frac{1}{h_{1}^{2}}\Big{(}\frac{1}{\tanh(\xi_{1})}\frac{\partial}{\partial\xi_{1}}+\frac{\partial^{2}}{\partial\xi_{1}^{2}}\Big{)}+\frac{1}{h_{2}^{2}}\Big{(}\frac{1}{\tan(\xi_{2})}\frac{\partial}{\partial\xi_{2}}+\frac{\partial^{2}}{\partial\xi_{2}^{2}}\Big{)}+\frac{1}{h_{3}^{2}}\frac{\partial^{2}}{\partial\xi_{3}^{2}},$
(69)

is a Laplacian operator that acts on individual Cartesian velocity component,
$u_{i}$ with derivatives defined in prolate spheroidal coordinates. In the
above equation

$h_{1}=h_{2}=f\sqrt{\sinh^{2}(\xi_{1})+\sin^{2}(\xi_{2})},h_{3}=f\sinh(\xi_{1})\sin(\xi_{2}).$
(70)

An advantage of this operator is that $\text{SphinCart}\widetilde{u}_{i}$ has
no cross-terms with other velocity components. This decoupling allows us to
solve the mass-momentum system of equations more efficiently. As mentioned in
section 3.2 the second step of the matrix-vector product in the Krylov
subspace procedure (equation (36)) within the Schur complement method involves
inversion of a discrete Laplacian operator acting on velocity vector. Unlike
in curvilinear basis, the vector Laplacian in Cartesian basis can be treated
as a set of three independent scalar Laplacians acting on three velocity
components. Hence, it allows us to use the algebraic multi-grid methods such
as AGMG [54], developed for scalar elliptic equations, to obtain the efficient
inversion of the vector Laplacian operator. In other words, this choice of
basis allows us to have a smaller bandwidth and number of off-diagonal
components of the discrete matrix operators in equation (31) as compared to
using a more natural choice of spheroidal basis.

As mentioned earlier in section 4.1, we use a staggered grid arrangement to
represent different variables: $\widetilde{\mathbf{u}}$ is stored on the
velocity and $\widetilde{p}$, $\boldsymbol{\Psi}$ and $\gamma$ on the pressure
grid. The momentum equation (29) is discretized on the velocity grid.
${\mathbf{u}}_{\infty}$ is an analytically known function that can be
evaluated at any location. The mass conservation equation (28), polymer
constitutive equation (30) and $\gamma$ equation (41) for FENE models are
discretized on the pressure grid. Similar to the finite difference
representation of the viscous terms by [24], we use a straightforward local
Lagrange polynomial representation of the spatial operators on a discretized
quantity. This allows flexibility in choosing the order of accuracy of finite
difference schemes in different directions and for different interpolation and
spatial derivative operators and is described in section 4.3. The spatially
and temporally discretized governing equations at each grid point and at time
step $n+1$ are,

$\displaystyle\widehat{\frac{\partial\xi_{k}}{\partial
x_{i}}}\widehat{\frac{\delta\widetilde{u}_{i}^{n+1}}{\delta\xi_{k}}}=0,$ (71)
$\displaystyle\Big{(}\frac{\rho_{f}}{\Delta
t}-\nabla^{2}_{\text{SphinCart}}\Big{|}^{\text{disc.}}\Big{)}\widetilde{u}^{n+1}_{i}+\overline{\frac{\partial\xi_{k}}{\partial
x_{i}}}\overline{\frac{\delta\widetilde{p}^{n+1}}{\delta\xi_{k}}}=g_{i}^{\text{disc.}},$
(72) $\displaystyle\frac{\Psi^{n+1}_{ij}}{\Delta
t}+\frac{1}{2}\hat{u}_{l}^{n+1}\widehat{\frac{\partial\xi_{k}}{\partial
x_{l}}}\frac{\delta{\Psi}_{ij}^{n+1}-\delta{\Psi}_{ij,\infty}^{n+1}}{\delta\xi_{k}}-\frac{\text{SR}(\boldsymbol{\Psi}^{n+1},\hat{\mathbf{u}}^{n+1})_{ij}^{\text{disc.}}}{2}=\frac{\text{SR}(\boldsymbol{\Psi}^{n},\hat{\mathbf{u}}^{n})_{ij}^{\text{disc.}}}{2}-\frac{1}{2}\hat{u}_{l}^{n}\widehat{\frac{\partial\xi_{k}}{\partial
x_{l}}}\frac{\delta{\Psi}_{ij}^{n}}{\delta\xi_{k}},$ (73)
$\displaystyle(\gamma^{n+1})^{2}+\gamma^{n+1}\Delta
t\Big{(}-\frac{\gamma^{n}}{\Delta
t}+\hat{u}_{j}^{n}\widehat{\frac{\partial\xi_{k}}{\partial
x_{j}}}\frac{\delta\gamma^{n}}{\delta\xi_{k}}+\frac{2{\Lambda^{n}_{ij}}}{L^{2}}\widehat{\frac{\partial\xi_{k}}{\partial
x_{i}}}\widehat{\frac{\partial
u_{j}}{\partial\xi_{k}}}+\frac{\gamma^{n}-0.5}{De\gamma^{n}}+\frac{1}{DeL^{2}}\alpha_{1}^{n}\Big{)}-\frac{\Delta
t}{2De}\alpha_{2}=0,$ (74)

$\mathbf{u}^{n+1}$ in the polymer equations (73) and (74) is the total
velocity i.e. the sum of velocity deviation $\widetilde{\mathbf{u}}^{n+1}$ and
undisturbed velocity $\mathbf{u}_{\infty}^{n+1}$ evaluated on the pressure
grid. In these equations the hat, $\hat{}$, represents the scenario when the
velocity or its spatial gradients are evaluated on the pressure grid.
Similarly, the overbar, $\bar{}$, represents the case when a pressure grid
variable or its gradient is evaluated on the velocity grid. Various terms
appearing in these equations are now defined. $\alpha_{1}^{n}=3/(L^{2}-3)$ for
FENE-P and $1.5/\gamma^{n}$ for FENE-CR model, and, $\alpha_{2}=1$ for FENE-P
and $1-3/L^{2}$ for FENE-CR model. These variables are irrelevant for
Oldroyd-B or Giesekus models as the $\gamma$ equation is unnecessary.

$\nabla_{\mathbf{x}}\boldsymbol{\xi}=\begin{bmatrix}\frac{\partial\xi_{1}}{\partial
x_{1}}&\frac{\partial\xi_{2}}{\partial x_{1}}&\frac{\partial\xi_{3}}{\partial
x_{1}}\\\ \frac{\partial\xi_{1}}{\partial
x_{2}}&\frac{\partial\xi_{2}}{\partial x_{2}}&\frac{\partial\xi_{3}}{\partial
x_{2}}\\\ \frac{\partial\xi_{1}}{\partial
x_{3}}&\frac{\partial\xi_{2}}{\partial x_{3}}&\frac{\partial\xi_{3}}{\partial
x_{3}}\end{bmatrix}=\begin{bmatrix}\frac{f\cosh(\xi_{1})\sin(\xi_{2})\cos(\xi_{3})}{h_{1}^{2}}&\frac{f\sinh(\xi_{1})\cos(\xi_{2})\cos(\xi_{3})}{h_{2}^{2}}&-\frac{\sin(\xi_{3})}{h_{3}}\\\
\frac{f\cosh(\xi_{1})\sin(\xi_{2})\sin(\xi_{3})}{h_{1}^{2}}&\frac{f\sinh(\xi_{1})\cos(\xi_{2})\sin(\xi_{3})}{h_{2}^{2}}&\frac{\cos(\xi_{3})}{h_{3}}\\\
\frac{f\sinh(\xi_{1})\cos(\xi_{2})}{h_{1}^{2}}&-\frac{f\cosh(\xi_{1})\sin(\xi_{2})}{h_{2}^{2}}&0\end{bmatrix},$
(75)

are spatial functions of the coordinate system that remain constant with time.
$\overline{\frac{\partial\xi_{i}}{\partial x_{j}}}$ and
$\widehat{\frac{\partial\xi_{k}}{\partial x_{i}}}$ are temporally constant
analytical functions evaluated at the velocity and pressure grid respectively.

$\nabla^{2}_{\text{SphinCart}}|^{\text{disc.}}\widetilde{u}^{n+1}_{i}=\frac{1}{h_{1}^{2}}\Big{(}\frac{1}{\tanh(\xi_{1})}\frac{\delta\widetilde{u}^{n+1}_{i}}{\delta\xi_{1}}+\frac{\delta^{2}\widetilde{u}^{n+1}_{i}}{\delta\xi_{1}^{2}}\Big{)}+\frac{1}{h_{2}^{2}}\Big{(}\frac{1}{\tan(\xi_{2})}\frac{\delta\widetilde{u}^{n+1}_{i}}{\delta\xi_{2}}+\frac{\delta^{2}\widetilde{u}^{n+1}_{i}}{\delta\xi_{2}^{2}}\Big{)}+\frac{1}{h_{3}^{2}}\frac{\delta^{2}\widetilde{u}^{n+1}_{i}}{\delta\xi_{3}^{2}},$
(76)

is discrete Laplacian operator acting at the velocity grid ($\xi_{1}$,
$\xi_{2}$, $h_{1}$, $h_{2}$ and $h_{3}$ are evaluated at each grid point on
the velocity grid).

${g}_{i}^{\text{disc.}}=\rho_{f}\big{(}\frac{{u}^{n}_{i}}{\Delta
t}-\frac{3}{2}\widetilde{\text{ADV}}(\widetilde{\mathbf{u}}^{n},\mathbf{u}^{n}_{\infty},\boldsymbol{\omega}_{p}^{n};\mathbf{r})_{i}+\frac{1}{2}\widetilde{\text{ADV}}(\widetilde{\mathbf{u}}^{n-1},\mathbf{u}^{n-1}_{\infty},\boldsymbol{\omega}_{p}^{n-1};\mathbf{r})\big{)}_{i}+\overline{\frac{\partial\xi_{k}}{\partial
x_{j}}}\overline{\frac{\delta(\Pi_{ji}^{n+1}-\Pi_{ji,\infty}^{n+1})}{\delta\xi_{k}}},$
(77)
$\widetilde{\text{ADV}}(\widetilde{\mathbf{u}}^{n},\mathbf{u}^{n}_{\infty},\boldsymbol{\omega}_{p}^{n};\mathbf{r})_{i}=(\widetilde{u}_{j}^{n}+{u}_{\infty,j}^{n})\overline{\frac{\partial\xi_{k}}{\partial
x_{j}}}\frac{\delta\widetilde{u}_{i}^{n}}{\delta\xi_{k}}+\widetilde{u}_{j}^{n}(\nabla\mathbf{u}_{\infty})_{ji}+2\epsilon_{ijk}{\omega_{j}^{n}}\widetilde{u}_{k}^{n},$
(78)

$\nabla\mathbf{u}_{\infty}$ is a spatially constant tensor, that is
analytically determined from the imposed flow and particle velocities along
with the particle orientation (equation (45)).
$\text{SR}(\boldsymbol{\Psi}^{n},\hat{\mathbf{u}}^{n})_{ij}|^{\text{disc.}},$
is defined in equation (26) with $\boldsymbol{\Psi}=\boldsymbol{\Psi}^{n}$ and
$\mathbf{u}=\hat{\mathbf{u}}^{n}$ where the velocity, $\hat{\mathbf{u}}^{n}$,
dependence appears through $\mathbf{L}=\nabla\mathbf{u}^{n}$ defined as,

${L}_{ij}={L}_{ij}^{n}=\widehat{\frac{\partial\xi_{k}}{\partial
x_{i}}}\widehat{\frac{\delta u_{j}^{n}}{\delta\xi_{k}}}.$ (79)

###  4.3 Finite Difference Schemes for spatial derivatives and interpolations

Values of various spatial derivatives and interpolated quantities appearing in
the discretized equations,
$\widehat{\frac{\delta\widetilde{u}_{i}^{n}}{\delta\xi_{j}}}$,
$\overline{\frac{\delta\widetilde{p}^{n}}{\delta\xi_{j}}}$,
$\frac{\delta({\Psi}_{ij}^{n+1}-{\Psi}_{ij,\infty}^{n+1})}{\delta\xi_{k}}$,
$\overline{\frac{\delta(\Pi_{ji}^{n+1}-\Pi_{ji,\infty}^{n+1})}{\delta\xi_{k}}}$,
$\frac{\delta\gamma^{n}}{\delta\xi_{j}}$, $\hat{u}_{i}^{n}$,
$\frac{\partial\widetilde{u}_{i}^{n}}{\partial\xi_{j}}$ and
$\frac{\delta^{2}\widetilde{u}^{n}_{i}}{\delta\xi_{j}^{2}}$, are approximated
by higher order finite difference or interpolation schemes. Lagrange
polynomials provide a convenient and flexible way to generate finite
difference weights, irrespective of the grid spacing and order of accuracy
[62]. Each variable, $\widetilde{u}_{i}$, $\widetilde{p}$, $\Psi_{ij}$,
$\Pi_{ij}$ and $\gamma$ is locally represented by an $(M-1)^{th}$ order
Lagrange polynomial by fitting the variable’s data on $M$ neighboring grid
points on which it is stored. Consider any of these variables to be
represented by $\phi(\xi_{1},\xi_{2},\xi_{3})$ which is stored at locations
$[\xi_{1,i}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}}]$,
$i,j,k,\in\mathbb{N}$ (which represent either the pressure or the velocity
grid defined in equation (67)). Interpolation along coordinate
$\xi_{i},i=[1,3]$, to a location $\widetilde{\xi}_{i}^{\text{dist.}}$, in
Lagrange basis polynomials is,

$\displaystyle\begin{split}\phi(\widetilde{\xi}_{1}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}})&=\Sigma_{q=n_{1}}^{M-n_{1}}l_{q}^{(0)}(\widetilde{\xi}_{1}^{\text{dist.}},{\xi}_{1,(.)}^{\text{dist.}})\phi(\xi_{1,q}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}}),\\\
\phi({\xi}_{1,i}^{\text{dist.}},\widetilde{\xi}_{2}^{\text{dist.}},\xi_{3,k}^{\text{dist.}})&=\Sigma_{q=n_{1}}^{M-n_{1}}l_{q}^{(0)}(\widetilde{\xi}_{2}^{\text{dist.}},{\xi}_{2,(.)}^{\text{dist.}})\phi(\xi_{1,i}^{\text{dist.}},\xi_{2,q}^{\text{dist.}},\xi_{3,k}^{\text{dist.}}),\\\
\phi({\xi}_{1,i}^{\text{dist.}},{\xi}_{2,j}^{\text{dist.}},\widetilde{\xi}_{3}^{\text{dist.}})&=\Sigma_{q=n_{1}}^{M-n_{1}}l_{q}^{(0)}(\widetilde{\xi}_{3}^{\text{dist.}},{\xi}_{3,(.)}^{\text{dist.}})\phi(\xi_{1,i}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,q}^{\text{dist.}}).\\\
&l_{q}^{(0)}(\widetilde{\xi}_{i}^{\text{dist.}},{\xi}_{i,(.)}^{\text{dist.}})=\prod\limits_{j=n_{1},j\neq
q}^{M-n_{1}}\frac{\widetilde{\xi}_{i}^{\text{dist.}}-\xi_{i,j}^{\text{dist.}}}{\xi_{i,q}^{\text{dist.}}-\xi_{i,j}^{\text{dist.}}},\end{split}$
(80)

where the set of discrete points
$\\{\xi_{i,q}^{\text{dist.}},q=n_{1},n_{1}+1,\cdots,M+n_{1},n_{1}\in\mathbb{N}\\}$
form an M point stencil around $\widetilde{\xi}_{i}^{\text{dist.}}$. The value
of the variable $\phi$ is thus interpolated from its representation at these
$M$ locations to a location represented by the point
$\widetilde{\xi}_{i}^{\text{dist.}}$. As expected, if
$\widetilde{\xi}_{i}^{\text{dist.}}$ is one of the locations in the $M$ point
stencil, $\xi_{i,q}^{\text{dist.}}$, the original value of $\phi$ is
recovered. Equation (80) is used to interpolate a variable from the pressure
or velocity grid to a point $\widetilde{\xi}_{i}^{\text{dist.}}$ which is not
one of $\xi_{i,q}^{\text{dist.}}$ locations but lies somewhere between these
in the continuous space. One can then choose to have
$\widetilde{\xi}_{i}^{\text{dist.}}$ represent a point on the other grid. The
discrete first and second derivatives in $\xi_{1}$ are,

$\displaystyle\begin{split}\frac{\delta\phi}{\delta\xi_{1}}(\widetilde{\xi}_{1}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}})&=\Sigma_{q=n_{1}}^{M+1-n_{1}}l_{q}^{(1)}(\widetilde{\xi}_{1}^{\text{dist.}},{\xi}_{1,(.)}^{\text{dist.}})\phi(\xi_{1,q}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}})\\\
\frac{\delta^{2}\phi}{\delta\xi_{1}^{2}}(\widetilde{\xi}_{1}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}})&=\Sigma_{q=n_{1}}^{M+1-n_{1}}l_{q}^{(2)}(\widetilde{\xi}_{1}^{\text{dist.}},{\xi}_{1,(.)}^{\text{dist.}})\phi(\xi_{1,q}^{\text{dist.}},\xi_{2,j}^{\text{dist.}},\xi_{3,k}^{\text{dist.}}),\end{split}$
(81)

where,

$\displaystyle\begin{split}&l_{q}^{(1)}(\widetilde{\xi}_{i}^{\text{dist.}},{\xi}_{i,(.)}^{\text{dist.}})=\Sigma_{p=n_{1},p\neq
q}^{M+1-n_{1}}\frac{1}{\xi_{i,q}^{\text{dist.}}-\xi_{i,p}^{\text{dist.}}}\prod\limits_{j=n_{1},j\neq(q,p)}^{M+1-n_{1}}\frac{\widetilde{\xi}_{i}^{\text{dist.}}-\xi_{i,j}^{\text{dist.}}}{\xi_{i,q}^{\text{dist.}}-\xi_{i,j}^{\text{dist.}}},\\\
&l_{q}^{(2)}(\widetilde{\xi}_{i}^{\text{dist.}},{\xi}_{i,(.)}^{\text{dist.}})=\Sigma_{r=n_{1},r\neq
q}^{M+1-n_{1}}\frac{1}{\xi_{i,q}^{\text{dist.}}-\xi_{i,r}^{\text{dist.}}}\Sigma_{p=n_{1},p\neq{(q,r)}}^{M+1-n_{1}}\frac{1}{\xi_{i,q}^{\text{dist.}}-\xi_{i,p}^{\text{dist.}}}\prod\limits_{j=n_{1},j\neq(q,p,r)}^{M+1-n_{1}}\frac{\widetilde{\xi}_{i}^{\text{dist.}}-\xi_{i,j}^{\text{dist.}}}{\xi_{i,q}^{\text{dist.}}-\xi_{i,j}^{\text{dist.}}},\\\
\end{split}$ (82)

The derivatives in $\xi_{2}$ and $\xi_{3}$ are similarly defined. In the first
and second order derivative operators of equation (82)
$\widetilde{\xi}_{i}^{\text{dist.}}$ may or may not be one of
$\xi_{i,q}^{\text{dist.}}$ locations. In the latter case it is somewhere
within these $\xi_{i,q}^{\text{dist.}}$ locations in the continuous space. The
interpolation and discrete differentiation in the above equations is expressed
as a discrete operator acting on the variable $\phi$.

Terms including variables on the other grid:
$\widehat{\frac{\delta\widetilde{u}_{i}^{n}}{\delta\xi_{j}}}$ is obtained by
the multiplying the first order discrete differential operator,
$l_{q}^{(1)}({\xi}_{j}^{\text{dist.,pres}},{\xi}_{j,(.)}^{\text{dist.,vel}})$
and two interpolation operators
$l_{q}^{(0)}({\xi}_{f}^{\text{dist.,pres}},{\xi}_{f,(.)}^{\text{dist.,vel}})$
for $\\{f=1,2,3,f\neq j\\}$ with $u_{i}^{n}$ in any order. Each operator in
this case acts on the velocity grid and produces an output on the pressure
grid in a particular dimension. Thus,
$\widehat{\frac{\delta\widetilde{u}_{i}^{n}}{\delta\xi_{j}}}$ is the
derivative of the velocity component $u_{i}$ at time step $n$ with respect to
the $\xi_{j}$ coordinate, evaluated on the pressure grid. Similarly,
$\overline{\frac{\delta\widetilde{p}^{n}}{\delta\xi_{j}}}$ is obtained by
multiplication of
$l_{q}^{(1)}({\xi}_{j}^{\text{dist.,vel}},{\xi}_{j,(.)}^{\text{dist.,pres}})$,
and
$l_{q}^{(0)}({\xi}_{f}^{\text{dist.,vel}},{\xi}_{f,(.)}^{\text{dist.,pres}})$
for $\\{f=1,2,3,f\neq j\\}$ with $\widetilde{p}^{n}$ in any order to give the
pressure derivative in $\xi_{j}$ direction, evaluated on the velocity grid.
$\overline{\frac{\delta(\Pi_{ji}^{n+1}-\Pi_{ji,\infty}^{n+1})}{\delta\xi_{k}}}$
is obtained in a way similar to pressure derivative after the polymer stress
$\Pi_{ji}^{n+1}$ and $\Pi_{ji,\infty}^{n+1}$ are obtained from the matrix
logarithm of the polymer configurations $\Psi_{ji}^{n+1}$ and
$\Psi_{ji,\infty}^{n+1}$ respectively using the method described just before
section 3.1. $\hat{u}_{i}^{n}$, the velocity component in $i^{th}$ coordinate
evaluated on the pressure grid, is obtained by three successive interpolations
of velocity deviation from the undisturbed value in each direction from
velocity to pressure grid and then adding the interpolated value to the
undisturbed velocity evaluated at the pressure grid. Therefore, the
interpolation errors are only incurred on the velocity deviation from the far
field. Since centered schemes introduce no numerical diffusion the stencils in
these cases consist of an even number of grid points centered (equal on each
side) around the output location, ${\xi}_{j}^{\text{dist.,vel}}$ or
${\xi}_{j}^{\text{dist.,pres}}$, i.e. $M$ is an even integer. Near the
boundaries of the computational domain, a centered stencil is not possible. In
this case we still use an $M$ point stencil, but with more points on one side
of the output location so as to keep all the stencil points within the
computation domain and avoid using any ghost nodes. This issue arises only in
the $\xi_{1}$ direction. $\xi_{3}$ coordinate is periodic about $\xi_{3}=0$ or
$2\pi$ and we utilize this periodicity. The $\xi_{2}$ coordinate also has no
external boundaries in the computational domain and we will discuss the
special treatment required near the coordinate system dependent internal
boundaries at $\xi_{2}=0$ and $\pi$ in section 4.3.1.

Non-convective terms including variables on the same grid: The terms
$\frac{\partial\widetilde{u}_{i}^{n}}{\partial\xi_{j}}$ and
$\frac{\delta^{2}\widetilde{u}^{n}_{i}}{\delta\xi_{j}^{2}}$ appearing in the
Laplacian operator,
$\nabla^{2}_{\text{SphinCart}}|^{\text{disc.}}{u}^{n+1}_{i}$, require a single
multiplication of the discrete derivative operator (first or second order)
with the respective variable. In this case, the input and output grid are the
same, i.e., the velocity grid. For
$\frac{\partial\widetilde{u}_{i}^{n}}{\partial\xi_{j}}$ and
$\frac{\delta^{2}\widetilde{u}^{n}_{i}}{\delta\xi_{j}^{2}}$ we use a stencil
of odd length as the output location is one of the stencil points and there
are equal number of points on either side of it. We deal with the $\xi_{1}$
boundaries in the similar way as described above for terms including variables
on the other grid.

Convective terms (HOUC schemes of Nourgaliev et al. (2007) [63]): The spatial
gradients appearing in the convective terms must be treated carefully to
prevent numerical instability. These include
$\frac{\delta({\Psi}_{ij}^{n+1}-{\Psi}_{ij,\infty}^{n+1})}{\delta\xi_{k}}$ and
$\frac{\delta\gamma^{n}}{\delta\xi_{j}}$ in equations (73) and (74), and,
$\frac{\delta\widetilde{u}_{i}^{n}}{\delta\xi_{k}}$ within
$\widetilde{\text{ADV}}(\widetilde{\mathbf{u}}^{n},\mathbf{u}^{n}_{\infty},\boldsymbol{\omega}_{p}^{n};\mathbf{r})_{i}$
of the discrete momentum equation (equations (72) and (77)). Using centered
difference schemes such as the one noted above for non-convective terms lead
to a numerical instability. A slight upwinding, where one more point upstream
is used relative to the downstream direction, allows one to maintain high
spatial accuracy while maintaining numerically smooth and stable solution [24,
63]. This is termed as a higher order upwind central scheme (HOUC). It is
straightforward to implement and introduces much less diffusion error compared
to WENO schemes used to discretize the convective term [63]. Near the
boundaries of the computational domain ($\mathbf{r}_{p}$ and
$\mathbf{r}_{\infty}$) we reduce the length of the stencil while maintaining
upwinding. The direction of the velocity term multiplying the convective
spatial gradient, $\tilde{u}_{j}^{n}\overline{\frac{\partial\xi_{k}}{\partial
x_{j}}}$, ${u}_{j,\infty}^{n}\overline{\frac{\partial\xi_{k}}{\partial
x_{j}}}$ or $\hat{u}_{j}^{n}\widehat{\frac{\partial\xi_{k}}{\partial x_{j}}}$,
determines the shape of the stencil used at the time step, $n$. This term is
equivalent to the velocity components in the spheroidal basis for $k=1$ and 2
i.e. the velocity along $\xi_{1}$ and $\xi_{2}$ directions respectively. For
$k=3$ this term differs from velocity in the azimuthal, $\xi_{3}$, direction
by a factor $h_{3}$ from equation (70) and instead we use the direction of
velocity in $\xi_{3}$ to determine the stencil shape. Therefore, similar to
the original HOUC proposition of [63] we choose the stencil shape based on the
direction of convection in each direction. HOUC was originally [63] proposed
for interface tracking in multiphase flows, but we have found it to improve
numerical stability of flow of viscoelastic fluid around a prolate spheroid as
well.

####  4.3.1 Internal boundaries in the $\xi_{2}$ direction

The only physical boundaries of the computational domain are the two $\xi_{1}$
surfaces at $\xi_{1}^{\text{surface}}$ and $\xi_{1}^{\infty}$. Due to the
choice of the coordinate system, $\xi_{2}$ is bounded by 0 and $\pi$ and
$\xi_{3}$ by 0 and $2\pi$, which requires one to consider appropriate boundary
conditions at these computational boundaries. As mentioned above, $\xi_{3}$ is
a periodic coordinate, and the points $\xi_{3}=0$ and $2\pi$ represent the
same physical location. Therefore periodic boundary condition in $\xi_{3}$
allows the interpolation mentioned above and discrete derivative operators in
the $\xi_{3}$ direction to be implemented near the boundaries in a way similar
to internal $\xi_{3}$ points. The boundary treatment of $\xi_{1}$ operators is
already described above. The $\xi_{2}$ coordinate has non-periodic internal
boundaries. Such a scenario appears at $r=0$ in the cylindrical coordinate
system and has previously been treated in the context of finite difference
schemes by several researchers [22, 23, 24, 64]. In this situation scenarios
the finite difference stencil near $r=0$ requires information at an unphysical
location $r<0$. In prolate spheroidal coordinates, the information is required
in the unphysical regions $\xi_{2}<0$ and $\xi_{2}>\pi$. Similar to these
previous works, a variable $\phi$ (representing $p$, $\gamma$ or a component
of $\mathbf{u}$ or $\mathbf{\Psi}$) at these locations is defined as,

$\displaystyle\begin{split}\phi(\xi_{1},\xi_{2},\xi_{3})&=\phi(\xi_{1},-\xi_{2},\xi_{3}+\pi),\hskip
14.45377pt\xi_{2}\leq 0\\\
\phi(\xi_{1},\xi_{2},\xi_{3})&=\phi(\xi_{1},2\pi-\xi_{2},\xi_{3}+\pi),\hskip
14.45377pt\xi_{2}\geq\pi.\end{split}$ (83)

In previous works, the coordinate system and choice of basis for the vectors
and tensors are the same, and care must be taken in the sign of certain
components of vectors or tensors when using these transformations around the
internal points/ axes [24]. However, in a Cartesian basis the transformations
of equation (83) are valid for any component of $\mathbf{u}$ or
$\mathbf{\Psi}$ (as well as $p$ or $\gamma$).

The issue of non-physical internal boundaries also appears in the spherical
coordinate system near the origin at $r=0$ and at the polar axis ($\theta=0$
and $\pi$). This has been treated in a recent work of Santelli et al. [21] by
writing momentum equations for a transformed vector $[u_{\phi}\hskip
3.61371ptu_{r}r^{2}\hskip 3.61371ptu_{\theta}\sin(\theta)]^{T}$ obtained from
the original velocity vector $[u_{\phi}\hskip 3.61371ptu_{r}\hskip
3.61371ptu_{\theta}]^{T}$. Zero transformed vector boundary conditions are
imposed on the internal boundaries for equations governing the transformed
vector by Santelli et al. [21].

####  4.3.2 Singular axis and boundary conditions

Another non-physical issue that arises due to the choice of the coordinate
system is the appearance of singular terms. $1/h_{3}$ (defined in equation
(70)), and $\frac{\partial\xi_{3}}{\partial x_{1}}$ and
$\frac{\partial\xi_{3}}{\partial x_{2}}$ (defined in equation (75)) have
$\sin(\xi_{2})$ in the denominator. These terms appear as a factor in front of
several terms in the governing equations. As expressed so far, these terms
lead to a coordinate system generated singularity at the axis on which
$\xi_{2}=0$ or $\pi$. However, these factors appear along with a spatial
gradient in the $\xi_{3}$ direction. The singular part of these terms can be
expressed as

$\text{Sing}=\frac{1}{\sin(\xi_{2})}\frac{\partial\phi}{\partial\xi_{3}}.$
(84)

On the singular axes $\phi$, representing either a scalar $p$ or $\gamma$ or a
Cartesian component of $\mathbf{u}$ or $\mathbf{\Psi}$, is unique. Therefore,

$\frac{\partial\phi}{\partial\xi_{3}}\Big{|}_{\xi_{2}=0}=\frac{\partial\phi}{\partial\xi_{3}}\Big{|}_{\xi_{2}=\pi}=0.$
(85)

It can be checked from the governing equations that wherever terms of the type
in equation (84) appear there is no other $\xi_{2}$ dependence in the
multiplicative factors. Therefore, using L’Hopital’s theorem

$\lim\limits_{\xi_{2}\rightarrow(0,\pi)}\text{Sing}=\frac{1}{\cos(\xi_{2})}\frac{\partial^{2}\phi}{\partial\xi_{2}\partial\xi_{3}},$
(86)

and the non-physical coordinate system generated singularity is removed at the
expense of introducing additional derivatives in the $\xi_{3}$ (azimuthal)
direction. This treatment of coordinate system generated axis singularities is
motivated by the treatment of the $1/r$ singularity at $r=0$ in cylindrical
coordinate system by [23] and [24]. The axis singularity appears in the
pressure grid where the points $\xi_{2,1}^{\text{dist.,pres}}=0$ and
$\xi_{2,N_{2}}^{\text{dist.,pres}}=\pi$ are defined. Therefore, the
singularity axis treatment is implemented in the discretization of the mass
conservation and polymer constitutive equations. Such singularity does not
appear on the velocity grid as it is staggered relative to the pressure grid.
This treatment is not required to discretize the momentum equation.

The previous methods of [23] and [24] that used L’Hopital’s theorem to account
for the coordinate system generated singularity used a curvilinear
(cylindrical) basis for the vectors. These vector components are multi-valued
at the singularity ($r=0$) and as a result the axis values of these components
were obtained using the relation between curvilinear and Cartesian components
on the axis. The Cartesian components were averaged/ interpolated using the
points off the singular axis. This leads to a difference between the
analytical equation obtained after L’Hopital’s theorem and that obtained using
the Taylor series expansion of the actual discretized equations (compare
equations (37) and (51) of [24]). As mentioned in those studies, strict energy
conservation cannot be obtained. Here, we use Cartesian components for the
vectors and tensors and do not encounter this problem.

This completes the description of the numerical method. In the next section,
we implement this method on several flows of inertia-less Newtonian and
viscoelastic fluids, along with examples of Newtonian fluids with inertia.

##  5 Numerical Tests

In this section, we use the numerical method described above to compute the
flow field around particles such as a sphere and a prolate spheroid in a
Newtonian fluid with and without (Stokes flow) inertia and an inertia-less
viscoelastic fluid. Several fluid-particle interaction problems are considered
for each class by changing the imposed velocity boundary conditions and torque
constraint on the particle. These include fixed and freely rotating spheres
and spheroids (including high aspect ratios) in uniform flow, simple shear
flow, extensional flow, and combinations of the first two. These studies show
forces, torques and stresslets, fluid streamlines, and the orientational
trajectories of freely suspended spheroids.

As mentioned in section 4.3 using Lagrange polynomials to represent finite
difference spatial discretization and interpolation allows us great
flexibility in choosing the order of spatial discretization schemes. In all
the examples presented, we have used a four-point stencil in $\xi_{1}$, a six-
point stencil in $\xi_{2}$ and an eight-point stencil in the $\xi_{3}$
direction. The higher-order discretization in $\xi_{2}$ and $\xi_{3}$ is
useful because for most cases (excluding section 5.3.2) we use less mesh
points in these directions as compared to $\xi_{1}$. After initial testing of
different orders of accuracy (not shown), we have found these orders of
accuracy to be adequate for validating the cases we present while obtaining
numerically stable simulations. In all the examples presented
$r_{\text{minor}}=1$.

###  5.1 Stokes flow: Motion of particles in inertia-less Newtonian fluid

One of our primary objectives in developing the numerical method described in
the previous sections is to study flows with zero to moderate inertia. The
interest in the former is in the presence of viscoelasticity. As discussed in
section 1, 3.2 and 3.5.2 zero inertia studies are important in isolating the
effect of viscoelastic fluids and our numerical method is suitable for
studying flow with zero particle and fluid inertia. In this section, we
compare our numerical method against the analytically available solution of
Stokes flow ($\rho_{f}=\rho_{p}=0$) of a Newtonian fluid ($c=0$) around
prolate spheroids.

####  5.1.1 Jeffery Orbits: A spheroid rotating in an inertia-less Newtonian
fluid

The orientation trajectory of a freely rotating prolate spheroid in a simple
shear flow of an inertia-less Newtonian fluid in the absence of particle
inertia is available from the work of Jeffery (1922) [65]. At a given time,
$t$, the polar angle, $\theta_{vort}$, with the vorticity axis of the imposed
shear flow, and the azimuthal angle, $\phi_{grad}$, with the gradient
direction in the flow-gradient plane, subtended by a prolate spheroid of
aspect ratio $\kappa$ are given by

$\tan(\phi_{grad})=\kappa\tan\Big{(}\frac{\dot{\gamma}t}{\kappa+\kappa^{-1}}\Big{)},\hskip
14.45377pt\tan(\theta_{vort})=\frac{C\kappa}{\sqrt{\kappa^{2}\cos^{2}(\phi_{grad})+\sin^{2}(\phi_{grad})}},$
(87)

where $\dot{\gamma}$ is the shear rate of the imposed flow and $C$ is an
orbital constant that depends on initial orientation. Hence, a prolate
spheroid undergoes three-dimensional periodic motion in a shear flow of a
Newtonian fluid, and the orbit’s shape is determined by initial orientation.
In figures 3a and 3b, we show these Jeffery orbits and the evolution of the
polar angle, $\phi_{grad}$, respectively, computed from our code along with
the analytical expressions from equation (87) for $\kappa=20$ at four starting
orientations. In the numerical solution, the domain size used is 10$\kappa$
and a uniform mesh is used with $N_{1}=120$, $N_{2}=N_{3}=71$ (these are
defined in equation (67)). Terms multiplying fluid and particle density are
ignored. At each time step, the resistivity formulation method described in
section 3.5.2 is used to determine the particle’s angular velocity, which
leads to zero net torque on the particle. The necessary Stokes flows (motion
induced part of the momentum equation) required in this formulation (as
mentioned in section 3.5.2) are pre-calculated numerically using the Schur
complement approach of section 3.2 before time evolving the particle’s
orientation. The numerical and the analytical curves shown in figure 3 are
indistinguishable.

(a)

(b)

Fig. 3: Jeffery Orbits: Prolate spheroid of aspect ratio $\kappa=20$ rotating
in simple shear flow of Newtonian fluid. Four different starting orientations
are plotted in (a), and the orientation trajectory follows one of the
degenerate periodic orbits depending upon the starting orientation. (b) shares
the same legend as (a) and shows the time evolution of the polar angle
$\phi_{grad}$ for the orbits of (a) (curves are indistinguishable as
$\phi_{grad}$ is independent of the orbit constant). We obtain good agreement
with the analytical prediction of Jeffery [65] for all orientations.

####  5.1.2 Forces on a fixed spheroid in a uniform flow

The drag, $C_{D}$, and the lift coefficient, $C_{L}$, of a prolate spheroid
with aspect ratio, $\kappa$, fixed in a uniform flow are given by [66, 67,
68],

$C_{L}=\frac{\text{Lift}}{\mu
U_{0}r_{\text{minor}}}=16\pi\kappa((K_{zz}-K_{xx})\cos^{2}(\theta)+K_{xx}),\hskip
14.45377ptC_{D}=\frac{\text{Drag}}{\mu
U_{0}r_{\text{minor}}}=16\pi\kappa(K_{xx}-K_{zz})\cos(\theta)\sin(\theta),$
(88)

where drag and lift are the hydrodynamic forces acting parallel and
perpendicular to the imposed flow in the plane of the imposed flow and the
particle center line (major axis). $\mu$ and $U_{0}$ are the fluid viscosity
and imposed speed, and $\theta$ is the angle of imposed flow relative to the
major axis of the particle.

$\begin{split}&K_{xx}=\frac{1}{\xi+\alpha},\hskip
14.45377ptK_{zz}=\frac{1}{\xi+\kappa^{2}\gamma},\\\
\alpha=\frac{\kappa^{2}}{\kappa^{2}-1}+\frac{\kappa}{2(\kappa^{2}-1)^{1.5}}\eta,&\hskip
14.45377pt\gamma=-\frac{2}{\kappa^{2}-1}-\frac{\kappa}{(\kappa^{2}-1)^{1.5}}\eta,\hskip
14.45377pt\eta=\log\Big{(}\frac{\kappa-\sqrt{\kappa^{2}-1}}{\kappa+\sqrt{\kappa^{2}-1}}\Big{)}.\end{split}$
(89)

(a)  (b)

Fig. 4: Lift and drag coefficients (89) on a prolate spheroid fixed in a
uniform flow at various angle of attacks ($\theta$) and particle aspect ratios
($\kappa$) in the inertia-less limit ($\rho_{f}=0$). We obtain good agreement
with the analytical prediction available in [66, 67, 68] for each $\kappa$ and
$\theta$.

We show the lift and drag coefficient for the flows past spheroids of three
different aspect ratios, $\kappa=6$, 10 and 20 in figure 4 for seven different
angles of attack between the imposed uniform flow and the particle’s major
axis. Good quantitative agreement between our numerical results and the
analytical expressions from equations (88) and (89) is obtained. The domain
size used is 100$\kappa$. A large computational domain is required for low
inertia flows past a solid particle because the velocity disturbance created
by the particle decays as $1/r$ at large distances ($r$) from the particle
[59, 45]. This issue has marred previous computations of Andersson et al. [69]
who mentioned the requirement of a large computational domain as the Reynolds
number in their study of $\kappa=6$ particle was reduced to about $0.1$. They
used a rectangular computational domain of size $64\kappa\times 64\kappa\times
42.7\kappa$. Depending upon the angle of attack, they reported a deviation of
0.6 to 15% from the analytical $C_{D}$ and $10.3-24.9$ % from the analytical
$C_{L}$ values for zero inertia. In contrast, using a very large computational
domain of size 100$\kappa$ with $N_{1}=250$, $N_{2}=80$ and $N_{3}=51$ we
obtain highly accurate results for $\kappa=6$, 10 and 20 shown in figure 4.
The largest absolute error is at $\theta=45^{\circ}$ for $\kappa=20$. It is
only 1.09% and 0.23% in $C_{L}$ and $C_{D}$ respectively. An accurate
simulation of these cases with a reasonable number of mesh points is possible
because while the required size of the computational domain is large, the grid
can be relatively sparse in the far-field and the prolate spheroidal grid is
naturally more clustered near the particle surface. Therefore, prolate
spheroidal coordinates employed in our simulations are an appropriate choice
for such studies, particularly when we test theories designed to be a
perturbation from Stokes flow. This example will be shown in section 5.2.1.

###  5.2 Particles in Newtonian fluid with finite inertia

####  5.2.1 Uniform flow past a fixed prolate spheroid

A prolate spheroid sediments in an inertial-less Newtonian fluid without any
change in orientation. Equivalently, a uniform flow of inertial-less Newtonian
fluid past a fixed spheroid, i.e., the cases of section 5.1.2, exerts no
hydrodynamic torque on the particle. Inertia, however, leads to a finite
torque on a fixed spheroid. In the low inertia limit, at the steady-state,
Dabade et al. (2015) [70] calculated the coefficient of inertial torque,
$C_{T}$, on a particle with aspect ratio, $\kappa$, fixed in a Newtonian
fluid. When the angle of attack of the oncoming flow relative to the particle
axis is $\theta$, they find,

$C_{T}=\frac{Re}{2}F(\kappa)\sin(2\theta),$ (90)

where, $C_{T}$ and the Reynolds number, $Re$, are,

$C_{T}=\frac{\text{Torque}}{\mu U_{0}r_{\text{minor}}^{2}\kappa^{2}},\hskip
14.45377ptRe=\frac{\rho_{f}U_{0}r_{\text{minor}}\kappa}{\mu},$ (91)

and $U_{0}$ is the imposed velocity. $F(\kappa)$ is a non-dimensional
parameter that depends on the particle aspect ratio and $F(6)=0.5458$ [70]. In
figure 5, we show the steady-state values of $C_{T}/Re$ at various $\theta$
for $\kappa=6$ and three different $Re=0.3$, 3.0 and 30.0 from our numerical
results along with those of Andersson and Jiang (2019) [69] and Jiang et al.
(2021) [71]. We also show the analytical expression $F(6)\sin(2\theta)/2$ of
Dabade et al. (2015) [70].

Due to the inertial screening (see [72, 73, 74] for numerical and experimental
evidence of this mechanism in a dilute suspension of sedimenting spheres), the
velocity disturbance due to the particle in the presence of fluid inertia
decays at a rate faster than the $1/r$ ($r$ is the distance from the particle)
decay characteristic of the inertia-less or Stokes limit (see section 5.1.2
and [59, 45]). However, the inertial screening length may still be large if
Reynolds number is small and a large computational domain may be required to
quantitatively assess the validity of theories developed for small inertial
corrections. In figure 5 we observe that our numerical results, performed with
a large computational domain size of
$||\mathbf{r}_{\infty}||_{2}\approx||\mathbf{r}^{\text{minor}}_{\infty}||_{2}=100\kappa$
are closer to the analytical prediction of Dabade et al. (2015) [70] for
$Re=0.3$ than the numerical results of Andersson, Jiang and co-workers
(2019,2021) [69, 71] performed with a rectangular domain of size
$64\kappa\times 64\kappa\times 42.7\kappa$ [69] or $34\kappa\times
34\kappa\times 34\kappa$ [71]. For 45∘, simulations of [69, 71] have a
deviation of $17\%$ from the analytical prediction of [70] at $Re=0.3$. We
have a deviation of $2.7\%$. In addition to larger domain size, we use a
straightforward boundary condition of a constant uniform velocity on the outer
boundary. However, the discussion in [69] points to a more involved boundary
condition treatment. At higher $Re=3.0$ and 30.0 our computations agree with
that of [69, 71] as shown in figure 5. For $Re=3.0$ we show another simulation
result at
$||\mathbf{r}_{\infty}||_{2}\approx||\mathbf{r}^{\text{minor}}_{\infty}||_{2}=50\kappa$.
The resolution of our computational grid is $N_{1}=200$, $N_{2}=131$ and
$N_{3}=65$ for both computational domains. The agreement with the results of
[69, 71] is better for the smaller of our computational domains at $Re=3.0$, a
further indication of the importance of domain size at small to moderate
inertia. For $\theta=45^{\circ}$ at $Re=3.0$ and 30.0 we simulated additional
cases with increased resolution $N_{1}=300$, $N_{2}=201$ and $N_{3}=91$ (not
shown) and found similar results as with our lower resolution. Due to the
nature of spheroidal coordinates (see equation (63) and more detailed
discussion in section 1 and 4.1 about Euclidean spacing of the spheroidal
grid), the domain size can be greatly increased without decreasing the
resolution near the particle surface significantly. Therefore, our method is
well equipped to study flow around a sedimenting particle or a particle fixed
in uniform flow in the limit of small to moderate inertia.

We also show the steady-state streamlines of the flow for $Re=0$ and $Re=30$
in figure 6 for $\theta=90^{\circ}$. The presence of inertia leads to the
formation of trailing edge vortices. The streamline pictures for $Re=30$ are
qualitatively similar to those shown by Andersson et al. [69].

(a)  Fig. 5: Uniform flow past a fixed spheroid: Steady-state torque
coefficient, $C_{T}$ normalized with Reynolds number, $Re$, from our numerical
simulations and those of Andersson and Jiang (2019) [69] and Jiang et al.
(2021) [71] at different angles of attack, $\theta$ and Reynolds numbers,
$Re=0.3$, 3.0 and 30 along with the analytical calculations of Dabade et al.
[70] is also shown. Our results agree with that of [69, 71] at larger $Re$ and
with [70] at $Re=0.3$.

(a)

(b)

(c)

(d)

Fig. 6: Streamlines of the Newtonian fluid flow around a fixed prolate
spheroid of aspect ratio, $\kappa=6$ at Reynolds number $Re=0$ (Stokes flow)
and $Re=30$ with major axis at 90∘ to the oncoming flow (left to right). For
each $Re$, two different views are shown. The $Re=30$ streamlines are the same
as those observed i in figures 5(d) and 6(d) of [69].

####  5.2.2 Effect of inertia on a freely rotating sphere in simple shear flow

In the previous section, we demonstrated the validity of our method to capture
the inertial effect on a fixed particle. In this section, we show an example
where particle and fluid inertia are moderate, and the particle is allowed to
rotate due to the hydrodynamic stresses acting on the particle surface as
treated by the method described in section 3.5.1. For this purpose, we
validate our code with the simulations of Bagchi and Balachandar (2002) [12]
where the particle is a sphere. We use $\kappa=1.001$ to represent a sphere in
the prolate spheroidal coordinate system. The sphere is not allowed to
translate but can rotate freely under the influence of the hydrodynamic
torque. The imposed fluid velocity is,

$\mathbf{u}_{\text{imposed}}^{\text{inertial}}=[U_{0}+\dot{\gamma}y\hskip
7.22743pt0\hskip 7.22743pt0]^{T},$ (92)

in the laboratory frame. Two different Reynolds numbers based on the uniform
flow speed, $U_{0}$, or the shear rate, $\dot{\gamma}$ can be defined in this
case. One of them i.e. the translational Reynolds number ($Re$) is the same as
defined earlier in equation (91) and the shear Reynolds number is,

$Re_{\dot{\gamma}}=\frac{\rho_{f}\dot{\gamma}r^{2}}{\mu},$ (93)

where $r$ ($=r_{major}$) is the radius of the sphere. $Re$ and
$Re_{\dot{\gamma}}$ are related by a non-dimensional factor,
$\dot{\gamma}r/U_{0}$. In Stokes flow ($Re=Re_{\dot{\gamma}}=0$) the sphere
rotates at the angular velocity of the fluid, $\omega_{f}=\dot{\gamma}/2$. The
presence of inertia lowers the rotation rate of the sphere. In figure 7a we
show the steady state rotation rate of the sphere, $\omega_{st}$, normalized
with $\omega_{f}$ for a range of $Re$ at two different
$\dot{\gamma}r/U_{0}=0.05$ and 0.1 from our simulations and those of Bagchi
and Balachandar (2002) [12]. We find a good, albeit not exact, agreement of
$\omega_{st}/\omega_{f}$ between the two simulations for all the cases shown
up to $Re=100$. Bagchi and Balachandar (2002) [12] found the inertial slow-
down, $\omega_{st}/\omega_{f}$, of a sphere’s rotation to be independent of
$\dot{\gamma}r/U_{0}$ for a constant $Re$. This is also captured by our
results in figure 7a as the results for the two $\dot{\gamma}r/U_{0}$ shown
nearly collapse. In figure 7b we show the streamlines of the fluid flow around
the sphere for $Re=100$ and $\dot{\gamma}r/U+0=0.1$. These are qualitatively
similar to the streamlines shown in figure 8(b) of Bagchi and Balachandar
(2002) [12].

(a)

(b)

Fig. 7: Comparison of the inertial slow down of a freely rotating sphere in a
flow $\mathbf{u}=[U_{0}+\dot{\gamma}y\hskip 7.22743pt0\hskip 7.22743pt0]^{T}$
from our simulations and those of Bagchi and Balachandar (2002) [12] ($U_{0}$
and $\dot{\gamma}$ are constants). The sphere is allowed to rotate freely but
not translate. The left figure shows the ratio of steady state rotation rate
of a sphere, $\omega_{st}$ to the fluid rotation rate,
$\omega_{f}=\dot{\gamma}/2$ at various translations Reynolds numbers and at
two different $\dot{\gamma}r/U_{0}=0.05$ and 0.1. Quantitative agreement in
$\omega_{st}/\omega_{f}$ between our simulations and those of [12] is found at
the $Re$ and $\dot{\gamma}r/U_{0}$ shown. The right figure shows the
streamlines for Re=100 and $\dot{\gamma}r/U_{0}=0.1$ that are qualitatively
similar to the respective streamlines in figure 8(b) of [12].

This example shows the capability of our numerical methodology to handle
moving particles in the presence of inertia.

###  5.3 Motion of particles in inertia-less viscoelastic fluid

In this section, we compare our numerical method with other numerical and
semi-analytical results of the flow of viscoelastic fluids around spheres and
prolate spheroids in various linear flows.

####  5.3.1 Torque-free rotation in simple shear flow

In section 5.2.2, we discussed the influence of inertia in slowing down the
rotation of a sphere in simple shear flow from its value in the Stokes limit.
Avino et al. (2008) [75] reported a similar effect due to viscoelasticity
where increasing $De$ in an inertia-less viscoelastic fluid lowers the
rotation rate of a spherical particle. A subset of these authors reported the
steady-state angular velocity of the sphere in Avino et al. (2014) [26]. Here
$De$ is the non-dimensional parameter known as the Deborah number,
representing the product of the imposed shear rate and the particle relaxation
time. The fluid and particle inertia are ignored, and the net torque due to
the fluid stresses acting on the surface of a freely moving particle is zero.
Hence the particle motion in the cases considered in this section is obtained
through the method described in section 3.5.2. The Giesekus constitutive
relation (equations (18) and (19)) with $c=10.0$ and $\alpha=0.2$ is used to
model the polymer stress in the viscoelastic fluid. In figure 8 we show that
the results from our simulations for a sphere rotating in a simple shear flow
of this fluid are almost identical to that of Avino et al. (2014) [26]. In our
simulations, a prolate spheroid with aspect ratio $\kappa=1.001$ represents
the sphere. Avino et al. (2014) [26] also shown the average angular velocity
for the case when the major axis of a prolate spheroid rotates in the flow-
gradient plane of a simple shear flow of the same viscoelastic fluid. We
compare the results for $\kappa=2.0$ from our simulations with that of Avino
et al. (2014) [26] in figure 8 and find an excellent agreement at all $De$
shown. For a spheroid in a Newtonian fluid, i.e. at $De=0.0$, according to
Jeffery (1922) [65] $\bar{\omega}/\omega_{f}=2.0\kappa/(1+\kappa^{2})=$ 1.0
and 0.8 for $\kappa=1.0$ (spheres) and 2.0 respectively. These analytical
estimates are the numerical values shown in figure 8 at $De=0.0$. We show the
streamlines around a freely rotating sphere in a Newtonian fluid and
viscoelastic fluid with $De=1.0$ in figures 9a and 9b. These are qualitatively
similar to those shown in figures 9(a) and 9(d) of Avino et al. (2008) [75].
The effect of viscoelasticity is to distort the region of closed streamlines
around the sphere that extends to infinity for Newtonian fluid. At $De=1.0$,
both the simulations find stagnation points on either side of the sphere in
the flow direction, which marks the end of the closed streamline region. By
showing streamlines for $De=$ 0.1, 0.3 and 1.0 Avino et al. (2008) [75]
established the movement of hte stagnation points closer to the sphere upon
increasing $De$. We also show the continuation of this trend to $De=3.0$. At
this higher $De$, an asymmetry of the region of closed streamlines and a
reverse wake about the flow (horizontal) direction is also observed.

Fig. 8: Slower particle rotation due to viscoelasticity. Period average
angular velocity, $\bar{\omega}$, (normalized with the fluid rotation rate,
$\omega_{f}=\dot{\gamma}/2$) of a torque-free sphere and a prolate spheroid
with aspect ratio, $\kappa=2.0$ in the flow-gradient plane of a simple shear
flow of Giesekus viscoelastic fluid ($c=10$ and $\alpha=0.2$) at various $De$.
Our results are quantitatively similar to that of Avino et al. (2014) [26] at
all $De$.

(a)

(b)

(c)

Fig. 9: Streamlines around a torque-free sphere rotating in a simple shear
flow of a Newtonian fluid (left) and of a Giesekus viscoelastic fluid (right)
with $c=10.0$, $\alpha=0.2$ and $De=$1.0. The changes in streamlines due to
viscoelasticity are consistent with the observations of Avino et al. (2008)
[75] (figures 9(a) and 9(d) of [75]). Streamlines for a Giesekus fluid with
$De=3.0$ are also shown that extend the conclusions of Avino et al. (2008)
[75] to higher $De$ than explored in their study.

Now we consider three dimensional rotation of a $\kappa=4.0$ prolate spheroid
in shear flow of the same viscoelastic fluid. As shown by the numerical
simulations of Avino et al. (2014) [26], in the presence of viscoelasticity
the orientation of a prolate spheroid no longer follows the Jeffery orbits as
in Newtonian Stokes flow [65] (section 5.1.1). In figure 10 we compare the
orientation trajectory of a $\kappa=4.0$ prolate spheroid rotating in an
inertia-less viscoelastic fluid from our simulations with that of Avino et al.
(2014) [26]. We show two different starting orientations for $De=1.0$ and 3.0
in figure 10. The resolution used in these simulations is $N_{1}=150$,
$N_{2}=71$ and $N_{3}=57$ (equation (67)). The size of the computational
domain used is
$||\mathbf{r}_{\infty}||_{2}\approx||\mathbf{r}^{\text{minor}}_{\infty}||_{2}=20\kappa$.
Our simulation results qualitatively agree with those of Avino et al. (2014)
[26] with small, subtle differences. Similar to Avino et al. (2014) [26], we
find the final orientation behavior of the particle at $De=1$ to be spiraling
towards the vorticity direction of the imposed simple shear flow, irrespective
of the initial orientation. At $De=3.0$, the particle in both simulations
settles to a location very close to the flow direction, irrespective of the
starting orientation.

(a)

(b)

Fig. 10: Orientation trajectory of a torque-free prolate spheroid with aspect
ratio, $\kappa=4$ released in a simple shear flow of an inertia-less Giesekus
viscoelastic fluid ($c=10$ and $\alpha=0.2$) at two different initial
orientations in the gradient-vorticity plane at $De=1.0$ and 3.0. Solid orange
and purple lines are from our simulations and the dashed grey and black are
from Avino et al. (2014) [26]. Grey arrows indicate the imposed flow in the
shearing plane. Good qualitative agreement is obtained between the two
results. The orientation trajectories of Avino et al. (2014) [26] were
obtained through personal communication with the authors.

####  5.3.2 Rheology of dilute suspensions of spheres

As mentioned in section 1, fluid flow around a particle in an unbounded fluid
is useful in studying the rheology of a dilute suspension of particles where
the inter-particle interaction is rare. In Newtonian fluids, the presence of
particles leads to an additional stress $n\boldsymbol{S}$ in the suspension,
where $n$ is the number of particles per unit volume and $\boldsymbol{S}$ is
termed stresslet. In rheology studies of incompressible suspensions, the
deviatoric or traceless part of the suspension stress is most interesting as
the trace of the stress can be absorbed in the modified pressure. The
deviatoric part of the stresslet, $\hat{\boldsymbol{S}}$, is an area integral
over the particle surface, $\mathbf{r}_{p}$, of a tensor product of the fluid
stress [76],

$\displaystyle\hat{\text{{S}}}(\boldsymbol{\sigma})=\int_{\mathbf{r}=\mathbf{r}_{p}}\text{d}A\big{\\{}\frac{1}{2}[{\text{{nn}}}\cdot\boldsymbol{\sigma}+\text{{n}}\cdot\boldsymbol{\sigma}\text{{n}}]-\frac{1}{3}\boldsymbol{\delta}\text{{n}}\cdot\boldsymbol{\sigma}\cdot\text{{n}}\big{\\}},$
(94)

where n is the unit surface normal pointing into the fluid. The stresslet also
appears in the particle suspension of viscoelastic fluids where