# On the Lattice Ground State of Densely Packed Hard Ellipses

S. Wagner1, G. Kahl1, R. Melnyk2 and A. Baumketner2<EMAIL_ADDRESS>1Institute for Theoretical Physics, TU Wien, Wiedner Hauptstraße 8-10, A-1040
Wien, Austria,
2Institute for Condensed Matter Physics, National Academy of Sciences of
Ukraine, 1 Svientsistsky Str., Lviv, UA-79011, Ukraine

(August 28, 2024)

###### Abstract

Among lattice configurations of densely packed hard ellipses, Monte Carlo
simulations are used to identify the so-called parallel and diagonal lattices
as the two favourable states. The free energies of these two states are
computed for several system sizes employing the Einstein Crystal method. An
accurate calculation of the free energy difference between the two states
reveals the parallel lattice as the state with lowest free energy. The origin
of the entropic difference between the two states is further elucidated by
assessing the roles of the translational and rotational degrees of freedom.

Since the early 1950s, when computer simulations were first used to tackle
problems of statistical mechanics numerically [1, 2], these methods have
provided an extremely valuable insight into the pivotal role that entropy
plays in self-assembling systems. Notably, the long-standing notion of entropy
being a thermodynamic potential driving any system into disorder was instantly
challenged when simulations of hard spheres first indicated [3, 4] (and later
confirmed [5, 6, 7, 8]) their crystallization, i.e., an ordering transition.
Self-assembly was seen to occur into two alternative states, fcc and hcp
lattices, with the question which of these two states possesses a lower free
energy evolving into a decade-long scientific debate - discussed vigorously in
many theoretical studies using a variety of numerical and analytical methods
[9, 10, 11]\- before being settled in favor of fcc lattice [12, 13, 14, 15,
16].

In systems of hard ellipses the situation becomes even more involved, as in
addition to the translational degrees of freedom (DOFs) also a rotational DOF
comes into play. Onsager was the first to point out the interplay between
these two types of DOFs in the system of hard rods [17], for which rotations
were seen to play an increasingly important role as the spatial elongation of
the particles increased, leading to a transition into the nematic phase, i.e.
a partially ordered phase in which the orientations of the particles are all
aligned while the respective positions remain disordered. Nematic phases are
also the focal point of many theoretical studies of hard-core ellipses [18,
19, 20, 21]. However, considerably less attention has been devoted to the
crystallization transition in this system that occurs upon increasing the
density and results in the formation of a densely packed ordered state. In
contrast to their shape-isotropic counterparts of hard disks, ellipses are
able to arrange into a multitude of different lattice states, some of which
were considered earlier [22, 20]. However, the nature of these different
lattice conformations is not yet understood. In particular, the questions of
how many stable structures can be expected and which of these is the most
stable one have not been addressed so far.

The present contribution is dedicated to answering exactly these questions.
Among all possible ellipse lattice states related to hexagonal order [18], we
identify - by means of Monte Carlo (MC) simulations - two stable candidates:
the so-called diagonal and the parallel lattice states, each of which is
characterized by the particles’ relative positions to one another as well as
by their orientations with respect to the lattice direction. We then proceed
to evaluate the free energy difference between the identified states by
explicit free energy methods, similarly to the approaches employed for the
hard sphere crystal problem [12]. Finally, one lattice is identified with the
lowest free energy - or equivalently with the highest entropy - which is
analyzed in terms of the particles’ translational and orientational DOFs
contributing to the lattice’s thermodynamic stability.

We consider crystalline, high-density states of hard elliptic particles of
aspect ratio $\kappa=2$. We name the described lattice states hexagonal-like
lattices as these structures can be derived via deforming a hexagonal lattice
of disks. In these hexagonal-like lattices all ellipses are oriented in the
same direction having six nearest neighbours each. The transformation of a
lattice of disks into a lattice of ellipses was first described by Vieillard-
Baron [18] at close-packed density $\rho_{\mathrm{cp}}$. Here we generalize
this approach to arbitrary densities $\rho^{*}$ (defining
$\rho^{\ast}=\frac{\rho}{\rho_{\mathrm{cp}}}$) by performing two subsequent
transformations, each of them being characterized by a parameter, which allows
us to parameterize the resulting lattice states as seen in Figure 1.

Figure 1:  Transformation of a close packed lattice of disks (left column)
into hexagonal-like lattices of ellipses. Each panel shows the unit cell of
the lattice, as defined by lattice parameter $\omega$ (rows) and $\tau$
(columns), and spanned by the unit cell vectors; see text. The colored
particles mark a pair of ellipses that is moved closer to one another via the
$\tau$-transformation.

First, we apply an affine transformation, $\mathbf{H_{1}}$ that deforms disks
into ellipses of aspect ratio $\kappa$, see Equation 1. The angle $\omega$ at
which this transformation is applied to the hexagonal lattice of disks can be
chosen arbitrarily with $\omega\in[0,\pi/6)$ as the lattice is characterized
by a dihedral symmetry $\mathcal{D}_{6}$. We introduce the angle $\omega$
related to the rotation $\mathbf{R}(\omega)$ as the first lattice parameter,
defining the transformation:

$\mathbf{H_{1}}(\kappa)\cdot\mathbf{R}(\omega)=\Big{(}\begin{smallmatrix}\sqrt{\kappa}&0\\\
0&\frac{1}{\sqrt{\kappa}}\end{smallmatrix}\Big{)}\cdot\Big{(}\begin{smallmatrix}\cos{\omega}&-\sin{\omega}\\\
\sin{\omega}&\cos{\omega}\end{smallmatrix}\Big{)}\text{.}$ (1)

Upon expanding the close-packed lattice state to lower densities
$\rho^{\ast}<1$, a second lattice parameter, $\tau$, is introduced as the
ratio between the expansion rates along the two spatial directions, visualized
in Figure 1. The transformation capturing the non-uniform expansion along the
two lattice directions is given by:

$\mathbf{H_{2}}(\rho^{\ast},\tau)=\frac{1}{\sqrt{\rho^{\ast}}}\Big{(}\begin{smallmatrix}\sqrt{\tau}&0\\\
0&\frac{1}{\sqrt{\tau}}\end{smallmatrix}\Big{)}.$ (2)

Figure 2:  PDFs as obtained from lattice MC simulations for determining the
optimal values $\omega_{\mathrm{opt}}$ and $\tau_{\mathrm{opt}}$ for various
densities $\rho^{\ast}$ (as labeled). Top panel: PDFs as functions of $\omega$
from simulations with $\omega$-moves at fixed $\tau=1$. We find two distinct
states, the diagonal ($\omega^{\mathrm{opt}}=0$ in blue) and the parallel
state ($\omega_{\mathrm{opt}}=\frac{\pi}{6}$ in red). Central and bottom
panels: PDFs as functions of $\tau$ from simulations with $\tau$-moves at
fixed $\omega_{\mathrm{opt}}=0$ (central) and
$\omega_{\mathrm{opt}}=\frac{\pi}{6}$ (bottom). Optimal values
$\tau_{\mathrm{opt}}$ are identified as maxima of the respective PDF.

We now are interested in the following question: which state, defined by
$(\omega_{\mathrm{opt}},\tau_{\mathrm{opt}})$, is entropically the most
favorable one? To answer this question, we apply a specifically designed
lattice MC method [13], that probes different lattice states by proposing
lattice-parameter MC-moves, $\Delta\omega$ and $\Delta\tau$, with a rate of 5%
during he course of a single particle move MC-simulation that involves
particles’ displacements and rotations. In a first step, we implement lattice
moves for $\omega$ at fixed $\tau=1$. In an ensemble of $N=100$ particles, we
perform a total of $60$ parallel runs, with initial lattice parameters
$\omega_{\mathrm{init}}$ evenly distributed within the range $[0,\pi/6)$. The
resulting probability density functions (PDFs) are shown in the top panel of
Figure 2. We see that stability occurs at two highest symmetry states, i.e.
for $\omega_{\mathrm{opt}}=0$ and $\omega_{\mathrm{opt}}=\frac{\pi}{6}$, which
we term the diagonal and parallel lattice, respectively. In contrast,
probabilities for most of the other $\omega$-values vanish. With this clear
indication of two very distinct optimal lattice states, the second
optimization step with respect to the parameter $\tau$ is applied, focusing on
the two optimal $\omega$-parameters only. This second optimization step leads
to PDFs shown in the central and the bottom panel of Figure 2, considering the
diagonal and parallel lattice states both in rectangular box geometries and at
various densities.

It is seen that as the density is increased, the distributions become more
sharply peaked around certain optimal values. This trend is consistent between
simulations with constant $\tau$ or $\omega$. Specifically, lattices were
constructed by replicating the unit cells of ellipses (as shown in Figure 1)
along the directions of the unit cell vectors, thereby creating skewed
simulation boxes defined by box side lengths $L_{\mathrm{x}}$ and
$L_{\mathrm{y}}$, as well as an enclosing angle $\alpha$. Instead of
optimizing the lattice state by varying the lattice parameters $\omega$ and
$\tau$, we now let the simulation box geometry adjust according to the
ellipses’ most favourable arrangement. We introduce MC-moves acting solely on
the box parameters - $L_{\mathrm{x}}$, $L_{\mathrm{y}}$ and $\alpha$ \- while
keeping the volume constant [23]. Our simulations uncover two dominant lattice
configurations at high densities: a) those defined by $\alpha_{\textrm{min}}$
that correspond to the diagonal state and b) those defined by
$\alpha_{\textrm{max}}$ corresponding to the parallel state. It can be shown
[23] that $\alpha_{\textrm{min}}=\cos^{-1}\sqrt{\frac{3}{3+\kappa^{2}}}$ and
$\alpha_{\textrm{max}}=\frac{\pi}{2}-\cos^{-1}\sqrt{\frac{3\kappa^{2}}{3+\kappa^{2}}}$,
where $\kappa$ is the aspect ratio.

In a subsequent step, we used the umbrella sampling method [24] to determine
the free energy difference between the diagonal and the parallel state. We
applied a biasing potential on $\alpha$ in order to force the lattice to find
its equilibrium configuration at the pinned value of the box angle. Multiple
bins were considered so as to construct a path connecting the two states.
Simulations were performed for each bin and then combined in the re-weighting
procedure [25] to produce the free energy profile as a function of the box
angle $\alpha$, shown in Figure 3 as $\Delta
A(\alpha)=A(\alpha)-A(\alpha_{\textrm{max}})$. Its appearance, featuring two
minima and one maximum confirms that there are two stable states. Furthermore,
the parallel state appears as the global free energy minimum while the
diagonal state corresponds to a local minimum. We find the free energy
difference (FD) between the two states, $\Delta
A_{\textrm{d}-\textrm{p}}=\Delta A(\alpha_{\textrm{min}})$ to be 2.0$\pm$0.9
$k_{\mathrm{B}}T\leavevmode\nobreak\ $.

Figure 3:  Free energy profile $\Delta A(\alpha)=A(\alpha)-A(\alpha_{\rm
max})$ shifted to zero at $\alpha=\alpha_{\rm max}$ as obtained for the system
composed of $N$=36 particles at $\rho^{*}=0.98$ in umbrella sampling
simulations. Error bars were estimated as standard deviations over 10
independent trajectories.

This value for the FD is small but statistically significant. Still, the FD
may result from a systematic error stemming from the finite size of the
simulated system ($N$=36) rather than representing the true free energy
difference between the two states. In an effort to determine which of the two
scenarios takes place - i.e., either finite-size effects or a real FD
sustained in the thermodynamic limit - we performed a finite-size scaling
analysis. To that end, $\Delta A_{\textrm{d}-\textrm{p}}$ was computed for
several system sizes $N$, with the aim of performing an extrapolation towards
the limit of $N\to\infty$. As the umbrella sampling simulations become cost-
prohibitive for large $N$, we instead used the Einstein crystal method [26],
more specifically an implementation proposed by Vega et al [27, 28]. In this
method, the total free energy $A$ of the target state, defined by Hamiltonian
$\mathcal{H}$, is computed by connecting it via a thermodynamic path to a
reference state of known free energy. This reference state is an Einstein
crystal (EC), i.e. a lattice of independent, classical harmonic oscillators
defined by an Hamiltonian $\mathcal{H}_{\mathrm{EC}}$ and a free energy
$A_{0}$. The EC is linked to the target lattice state via an intermediate
mixed state defined by the Hamiltonian
$\mathcal{H}_{\lambda}=\mathcal{H}+(1-\lambda)\mathcal{H}_{\mathrm{EC}}$, with
$\lambda$ as a tuning parameter. The free energy change $\Delta A_{1}$
associated with passing from the EC to the mixed state $\mathcal{H}_{\lambda}$
with $\lambda=0$ is computed using the umbrella sampling method [29]. In a
second step, the harmonic potential is switched off (via $\lambda\to 1$), thus
taking the system to the original Hamiltonian and completing the integration
path. The associated free energy cost of the last step is expressed via the
following integral:

$\displaystyle\Delta
A_{2}=-\int_{\lambda=0}^{1}\Big{\langle}\sum_{i=1}^{N}(\Lambda_{\rm E}^{\rm
t}\Delta\mathbf{r}_{i}^{2}+\Lambda_{\rm E}^{\rm
r}\Delta\phi_{i}^{2})\Big{\rangle}_{\lambda}d\lambda\text{ ,}$ (3)

where $\Delta\mathbf{r}_{i}=\mathbf{r}_{i}-\mathbf{r}_{i,0}$ and
$\Delta\phi_{i}=\phi_{i}-\phi_{i,0}$ describe the translational and
orientational displacements of particle $i$ with respect to the defined
equilibrium lattice states $\mathbf{r}_{i,0}$ and $\phi_{i,0}$. The constants
$\Lambda_{\textrm{E}}^{\textrm{t}}$ and $\Lambda_{\textrm{E}}^{\textrm{r}}$
are adjustable parameters that set the strength of the harmonic force in the
EC. Finally, the total free energy is expressed as the sum of three terms $A$
= $A_{0}$ \+ $\Delta A_{1}$ \+ $\Delta A_{2}$.

Since the EC is the same reference state for both the parallel and diagonal
lattice, the term $A_{0}$ drops from the difference $\Delta A_{\mathrm{d-p}}$.
Further, the contribution $\Delta A_{1}$ can be made numerically very small by
an appropriate choice of the constants $\Lambda_{\textrm{E}}^{\textrm{t}}$ and
$\Lambda_{\textrm{E}}^{\textrm{r}}$. In our computations, this leaves the
integral $\Delta A_{2}$ as the only non-zero contribution to the free energy
difference $\Delta A_{\mathrm{d-p}}$. We compute this difference for boxes of
two different geometries, a rectangular and a skewed simulation box, at a
density $\rho^{\ast}=0.98$. For the rectangular box we construct lattices
using the optimal lattice parameters $\omega_{\mathrm{opt}}$ and
$\tau_{\mathrm{opt}}$. The considered systems varied in size from $N=32$ to
$N=398$ particles. The number of particles in each lattice was chosen such
that a box aspect ratio of $\sim 1$ is achieved, making the boxes almost
square-shaped. For the skewed boxes the optimal geometry was extracted from
simulations with fixed box angles $\alpha$ at $\alpha_{\textrm{min}}$ and
$\alpha_{\textrm{max}}$. The number of particles was chosen as $N=n^{2}$,
where $n$ varies from 6 to 18.

We evaluate $\Delta A_{2}$ in Equation 3 numerically with the results shown in
Figure 4. For each box type we perform a finite-size analysis by fitting
$\beta\Delta A_{\mathrm{d-p}}(N)$ to $f(N)=kN^{b-1}+c$, where constants $k$,
$b$, and $c$ are treated as adjustable parameters. For both boxes the constant
$b$ was found to be close to 2, suggesting that the linear function
$f(N)=kN+c$ is able to represent the simulation data best. We therefore
proceed with a linear regression model. In the thermodynamic limit
$N\rightarrow\infty$, the FD per particle $\Delta
A^{\ast}_{\mathrm{d-p}}=N^{-1}\Delta A_{\mathrm{d}-\mathrm{p}}$ equates to the
constant $k$. The linear regression applied to our simulated data as shown in
Figure 4 yields $\Delta A^{\ast}_{\mathrm{d}-\mathrm{p}}=(0.0136\pm
0.0023)$$k_{\mathrm{B}}T\leavevmode\nobreak\ $for the square box and $\Delta
A^{\ast}_{\mathrm{d-p}}=(0.0176\pm
0.0022)$$k_{\mathrm{B}}T\leavevmode\nobreak\ $for the skewed box with a
confidence of $90\%$. Confirming the high credibility of the reported results,
the estimates are in agreement within error bars, strongly suggesting that the
FD between the diagonal and parallel state is real and not an artefact of the
employed simulation protocol.

Figure 4:  FD between the diagonal (d) and parallel (p) lattices at
$\rho^{\ast}=0.98$ for different system sizes $N$ and in two box geometries.
At finite sizes the box geometry has a significant effect on the FD. When
applying linear regression models to obtain the thermodynamic limit, given by
constant $k$ results in both box geometries agree within error bars.

Next, we focus on the microscopic origin of this FD. Equation 3 provides an
excellent tool for this task as it enables a clear separation of the total
free energy into contributions originating from the translational and the
orientational DOFs, i.e. $\Delta
A^{\mathrm{t}}_{2}=-\int\langle\sum_{i}^{N}\Lambda^{\rm t}_{\rm
E}\Delta\mathbf{r}_{i}^{2}\rangle_{{}_{\lambda}}d\lambda$ and $\Delta
A^{\mathrm{r}}_{2}=-\int\langle\sum_{i}^{N}\Lambda^{\rm r}_{\rm
E}\Delta\phi_{i}^{2}\rangle_{{}_{\lambda}}d\lambda$, respectively. This
separation will facilitate a mechanistic understanding of the origin of the
FD. It is important to point out that there is no strict distinction between
the translational and rotational contributions to the FD as the DOFs are
coupled. Nevertheless, $\Delta A^{\mathrm{t}}_{2}$ and $\Delta
A^{\mathrm{r}}_{2}$ can be seen as a useful measure of the phase-space volume
roughly accessible to the respective DOFs. We also note that these integrals
depend on the actual values of $\Lambda^{\textrm{t}}_{\textrm{E}}$ and
$\Lambda^{\textrm{r}}_{\textrm{E}}$; still, we find that this dependence is
weak and does not affect the conclusions at the qualitative level. As an
example, we report that $\Delta A_{2}^{\mathrm{t,d-p}}=\Delta
A_{2}^{\mathrm{t},\textrm{d}}-\Delta A_{2}^{\textrm{t},\textrm{p}}=(0.260\pm
0.041$) $Nk_{\mathrm{B}}T\leavevmode\nobreak\ $for the translational DOFs and
$\Delta A_{2}^{\mathrm{r,d-p}}=\Delta A_{2}^{\textrm{r},\textrm{d}}-\Delta
A_{2}^{\textrm{r},\textrm{p}}=(-0.235\pm 0.022$)
$Nk_{\mathrm{B}}T\leavevmode\nobreak\ $for the rotational DOF for the system
with $N=392$ particles and a particular set of $\Lambda_{\rm E}$’s.
Remarkably, the contributions come with a differing sign, suggesting that the
translational DOFs enhance the stability of the parallel state while the
rotational DOF seem to destabilize it.

We performed further tests to elucidate the role of different DOFs. The FD for
a system of $N=36$ particles in a skewed simulation box is $\Delta
A_{\textrm{d}-\textrm{p}}=(2.0\pm 0.1)$$k_{\mathrm{B}}T\leavevmode\nobreak\ $.
To probe the effect of particle rotations on this value, we repeated the
computations for the same system but allowing only particle translations, i.e.
freezing rotational DOFs. Here we obtained a FD $\Delta
A_{\textrm{d}-\textrm{p}}=(5.6\pm 0.1)$$k_{\mathrm{B}}T\leavevmode\nobreak\ $,
indicating that the parallel state does offer a much greater translational
freedom to the ellipses. In a similar manner, when we freeze the translational
DOFs, the resulting FD becomes a non-negative value $\Delta
A_{\textrm{d}-\textrm{p}}=(0.6\pm 0.1)$$k_{\mathrm{B}}T\leavevmode\nobreak\ $,
opposed to the analysis based on Equation 3. Interestingly, in the system with
frozen translational DOFs the parallel lattice turns out to offer a larger
free volume accessible through rotations, when in the system with all DOFs
active it is the diagonal lattice. This result highlights the impact of the
coupling of the rotational to the translational DOFs as it drastically raises
the rotational entropy within the diagonal lattice. Overall, however, the
higher rotational entropy in the diagonal lattice cannot outweigh the
advantage achieved by the translational DOFs within the parallel lattice.

The importance of the rotational DOF is further highlighted when we again
compute FD in a system of ellipses whose orientational DOFs are frozen, but
now in a simulation box that has relaxed to the equilibrium geometry specific
to the system with frozen orientational DOFs. Here we found that the FD among
all states, including the parallel and diagonal one, vanishes completely.

Taken together, these results show that ellipses require both translations and
rotations in order to develop a preferred thermodynamic state at high
densities. Further, they elevate the role of coupling of rotational and
translational DOFs in letting the parallel lattice state develop this
thermodynamic preference.

## References

  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Journal of Chemical Physics 21, 1087 (1953).
  * Rosenbluth and Rosenbluth [1954] M. N. Rosenbluth and A. W. Rosenbluth, Journal of Chemical Physics 22, 881 (1954).
  * Alder and Wainwright [1957] B. J. Alder and T. E. Wainwright, Journal of Chemical Physics 27, 1208 (1957).
  * Wood and Jacobson [1957] W. W. Wood and J. D. Jacobson, Journal of Chemical Physics 27, 1207 (1957).
  * Pusey and van Megen [1986] P. N. Pusey and W. van Megen, Nature 320, 340 (1986).
  * Harland and van Megen [1997] J. Harland and W. van Megen, Phys. Rev. E 55, 3054 (1997).
  * Simeonova and Kegel [2004] N. B. Simeonova and W. K. Kegel, Physical Review Letters 93, 035701 (2004).
  * Zaccarelli  _et al._ [2009] E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, M. E. Cates, and P. N. Pusey, Physical Review Letters 103, 135704 (2009).
  * Radin and Sadun [2005] C. Radin and L. Sadun, Physical Review Letters 94, 015502 (2005).
  * Woodcock [1997a] L. V. Woodcock, Letters to Nature 343, 141 (1997a).
  * Woodcock [1997b] L. V. Woodcock, Nature 388, 236 (1997b).
  * Noya and Almarza [2015] E. G. Noya and N. G. Almarza, Molecular Physics 113, 1061 (2015).
  * Bruce  _et al._ [1997] A. D. Bruce, N. B. Wilding, and G. J. Ackland, Phys. Rev. Lett. 79, 3002 (1997).
  * Bolhuis  _et al._ [1997] P. G. Bolhuis, D. Frenkel, S. C. Mau, and D. Huse, Nature 388, 235 (1997).
  * Mau and Huse [1999] S. C. Mau and D. A. Huse, Physical Review E 59, 4396 (1999).
  * de Miguel  _et al._ [2007] E. de Miguel, R. G. Marguta, and E. M. del Rio, Journal of Chemical Physics 127, 154512 (2007).
  * Onsager [1949] L. Onsager, Annals of the New York Academy of Sciences 51, 627 (1949).
  * Vieillard-Baron [1972] J. Vieillard-Baron, Journal of Chemical Physics 56, 24501 (1972).
  * Frenkel and Eppenga [1985] D. Frenkel and R. Eppenga, Physical Review A 31, 1776 (1985).
  * Cuesta and Frenkel [1990] J. A. Cuesta and D. Frenkel, Physical Review A 42, 2126 (1990).
  * Maeda and Maeda [2003] H. Maeda and Y. Maeda, Physical Review Letters 90, 4 (2003).
  * Bautista-Carbajal and Odriozola [2014] G. Bautista-Carbajal and G. Odriozola, Journal of Chemical Physics 140, 204502 (2014).
  * Baumketner [2022] A. Baumketner, Condensed Matter Physics 25, 33201 (2022).
  * Torrie and Valleau [1977] G. M. Torrie and J. P. Valleau, Journal of Computational Physics 23, 187 (1977).
  * Kumar  _et al._ [1992] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and J. M. Rosenberg, Journal of Computational Chemistry 13, 1011 (1992).
  * Frenkel and Ladd [1984] D. Frenkel and A. J. C. Ladd, Journal of Chemical Physics 81, 3188 (1984).
  * Vega and Noya [2007] C. Vega and E. G. Noya, Journal of Chemical Physics 127, 154113 (2007).
  * Noya  _et al._ [2008] E. G. Noya, C. Vega, and E. D. Miguel, Journal of Chemical Physics 128, 154507 (2008).
  * Frenkel and Smit [2023] D. Frenkel and B. Smit,  _Understanding Molecular Simulation_ (Elsevier Science, 2023).


# Implicit-explicit schemes for incompressible flow problems with variable
viscosity

Gabriel R. Barrenechea 1, Ernesto Castillo 2
and Douglas R. Q. Pacheco 3, 000Email address<EMAIL_ADDRESS>corresponding author
1 Department of Mathematics and Statistics, University of Strathclyde,
Glasgow, Scotland
2 Department of Mechanical Engineering, University of Santiago de Chile,
Santiago, Chile
3 Department of Mathematical Sciences, NTNU, Trondheim, Norway

###### Abstract

In this work we study different Implicit-Explicit (IMEX) schemes for
incompressible flow problems with variable viscosity. Unlike most previous
work on IMEX schemes, which focuses on the convective part, we here focus on
treating parts of the diffusive term explicitly to reduce the coupling between
the velocity components. We present different, both monolithic and fractional-
step, IMEX alternatives for the variable-viscosity Navier–Stokes system,
analysing their theoretical and algorithmic properties. Stability results are
proven for all the methods presented, with all these results being
unconditional, except for one of the discretisations using a fractional-step
scheme, where a CFL condition (in terms of the problem data) is required for
showing stability. Our analysis is supported by a series of numerical
experiments.

###### keywords:

Incompressible flow , IMEX methods , Variable viscosity , Generalised
Newtonian fluids , Finite element method , Temporal stability

††journal: ArXiv

##  1 Introduction

Variable-viscosity flow problems are relevant in many physical and
technological processes. Viscosity variations can be produced, for example, by
temperature or pressure gradients, by non-Newtonian rheological behaviour, or
by the interaction of multiple fluid phases. Numerically, non-constant
viscosity often leads to extra non-linearities and ill-conditioning in the
equations, which can affect the convergence and the performance of linear and
nonlinear solvers [1, 2, 3]. These issues get amplified in the time-dependent
situation, where those solvers are used multiple (even thousands) of times as
time advances. As a consequence, the last few years have seen an expansion in
the literature on numerical methods for variable-viscosity incompressible flow
problems [2, 3, 4, 5, 6, 7].

Over the last few decades, several alternatives have been proposed to reduce
the computational complexity of fluid simulations. Those date back to the
seminal work [8], as well as operator-splitting schemes [9] and more general
projection methods, such as fractional-step schemes [4, 10, 11, 12, 13, 14]. A
more recent trend, which can also be combined with fractional stepping, are
implicit-explicit (IMEX) methods. These have become increasingly popular in
applications ranging from waves to fluids [15, 16, 17, 18, 19].

In the context of incompressible flows, IMEX usually refers to temporal
discretisations that treat convection explicitly or semi-implicitly, while
keeping viscous terms fully implicit [20], which is the case for most of the
IMEX literature for flow problems. To the best of our knowledge, only few
articles so far have considered IMEX methods in the presence of variable
viscosity. It was recently proposed in [21, 18] to augment the viscous term
with a grad-div stabilisation [22], treated explicitly. Another alternative
consists on adding and subtracting a large enough constant-coefficient viscous
term on one side of the equation, and then applying the IMEX time-marching
method to the equivalent equation. This is the approach used in [23, 24, 25],
showing good numerical performance, although the analysis carried out so far
is limited to scalar equations with constant coefficients.

For flow problems with variable viscosity, treating the diffusive term semi-
implicitly can present significant computational advantages. As a matter of
fact, for those flow problems the viscous term cannot be written as the
standard Laplacian. This induces a coupling of the different components of the
velocity (due to the presence of the symmetric gradient), which makes the
linear solvers more involved and memory- and time-consuming. Thus, the
possibility of devising IMEX discretisations that make the transposed velocity
gradient explicit, leading to linear problems with the same sparsity pattern
as a scalar transport equation, has practical (and theoretical) appeal.

In this work we further develop an idea explored in [3], where an alternative
formulation using the viscosity gradient was used to propose an IMEX scheme
decoupling the velocity components; the results indicate an improved temporal
stability of the modified formulation in comparison to “naive” IMEX
discretisations. More precisely, we study a variety of IMEX discretisations,
both monolitic and fractional-step, arising from making one part of the
viscous term explicit, while keeping the “Laplacian” part implicit. We prove
stability for all the methods studied, but the norm with respect to which the
stability is proven varies according to the formulation. In fact, we show that
the most “natural” way of splitting the viscous term leads to a weak
stability, while a more careful rewriting of the term to be made explicit
enhances it, which is confirmed by our numerical studies. To keep the
technical details to a minimum, and focus on the effects of the time
discretisation, we only analyse the semi-discretised (in time) problem, and
consider only temporally first-order schemes, both for the monolithic and the
fractional-step alternatives.

The rest of the manuscript is organised as follows. In Section 2 we introduce
the model problem, notation and useful analytical tools. In Sections 3 and 4
we derive stability estimates for monolithic and fractional-step methods,
respectively. Section 5 summarises the theoretical results of the paper, and
in Section 6 we briefly discuss spatial discretisation matters. We test and
compare numerically the different alternatives presented in Section 7, and we
draw some final remarks in Section 8.

##  2 Preliminaries

###  2.1 Model problem

Let us consider a finite time interval $(0,T]$ and a fluid domain
$\Omega\subset\mathbb{R}^{d}$, $d=2$ or $3$, with Lipschitz boundary
$\Gamma=\partial\Omega$. As a model problem, we consider the variable-
viscosity incompressible Navier–Stokes equations:

$\displaystyle\partial_{t}\mbox{\boldmath$u$}+(\nabla\mbox{\boldmath$u$})\mbox{\boldmath$u$}-\nabla\cdot[2\nu(\mbox{\boldmath$x$},t)\nabla^{\mathrm{s}}\mbox{\boldmath$u$}]+\nabla
p$ $\displaystyle=\mbox{\boldmath$f$}$ $\displaystyle\text{in}\ \
\Omega\times(0,T]=:Q\,,$ (1) $\displaystyle\nabla\cdot\mbox{\boldmath$u$}$
$\displaystyle=0$ $\displaystyle\text{in}\ \ \Omega\times(0,T]\,,$ (2) $u$
$\displaystyle=\mbox{\boldmath$0$}$ $\displaystyle\text{on}\ \
\Gamma\times(0,T]\,,$ (3) $u$ $\displaystyle=\mbox{\boldmath$u$}_{0}$
$\displaystyle\text{at}\ \ t=0\,,$ (4)

where the unknowns are the velocity $u$ and the pressure $p$, while the
remaining quantities are problem data. The viscosity field
$\nu(\mbox{\boldmath$x$},t)$ is assumed to satisfy

$\nu(\mbox{\boldmath$x$},t)\geq\nu_{\mathrm{min}}>0\ \ \text{in}\ \bar{Q}\,,$
(5)

where $\nu_{\mathrm{min}}$ is a known constant that represents, e.g., the
minimum viscosity of a shear-thinning fluid. To simplify notation, we shall
use $\nu$ instead of $\nu(\mbox{\boldmath$x$},t)$ throughout this work (but
will keep in mind its variable character).

###  2.2 Alternative formulations of the viscous term

In incompressible flows, the divergence-free constraint (2) leads to multiple
ways of describing the viscous term. Although popular in the constant
viscosity scenario, the Laplacian form $\nu\Delta\mbox{\boldmath$u$}$ is not
appropriate for problems with variable viscosity. For that reason, the full
stress-divergence (SD) form

$\nabla\cdot\left(2\nu\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\right)=\nabla\cdot(\nu\nabla\mbox{\boldmath$u$})+\nabla\cdot(\nu\nabla^{\top}\mbox{\boldmath$u$})$
(6)

is normally used when $\nabla\nu\not=\mbox{\boldmath$0$}$. Yet, other
possibilities exist, such as [18]

$\displaystyle\nabla\cdot\left(2\nu\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\right)=\nabla\cdot\left(2\nu\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\right)-\nabla(\nu\nabla\cdot\mbox{\boldmath$u$})\,,$
(7)

or also [6]

$\displaystyle\nabla\cdot\left(2\nu\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\right)$
$\displaystyle=2\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\nabla\nu-\nu\nabla\times(\nabla\times\mbox{\boldmath$u$})\,.$
(8)

In this work we consider the following rewriting of the viscous term

$\nabla\cdot\left(2\nu\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\right)\equiv\nabla\cdot\left(\nu\nabla\mbox{\boldmath$u$}\right)+\nabla^{\top}\mbox{\boldmath$u$}\nabla\nu+\nu\nabla(\nabla\cdot\mbox{\boldmath$u$})\\\
=\nabla\cdot\left(\nu\nabla\mbox{\boldmath$u$}\right)+\nabla^{\top}\mbox{\boldmath$u$}\nabla\nu\,,$
(9)

which we denote as generalised Laplacian (GL). Although equivalent at the
continuous level, the writings (6)-(9) can (and do) lead to dramatically
different numerical performances. As mentioned in the introduction, the
symmetric gradient leads to a coupling of the $d$ velocity components
(especially if $\nabla\nu\not=\mbox{\boldmath$0$}$). In this work our aim is
to make the parts of the viscous term responsible for this coupling explicit,
and to analyse the implications of such IMEX approaches.

###  2.3 Useful notation, identities and inequalities

We consider usual notation for Hilbert and Lebesgue spaces aligned, e.g., with
[26]. More precisely, we denote by $L^{2}_{0}(\Omega)$ the space of
$L^{2}(\Omega)$ functions with zero mean on $\Omega$. The inner product in
$L^{2}(\Omega)$ is denoted by $(\cdot,\cdot)$, and we make no distinction
between the product of scalar or vector/tensor-valued functions. Moreover, we
denote by $\|\cdot\|$ and $\|\cdot\|_{\infty}$ the $L^{2}(\Omega)$ and
$L^{\infty}(\bar{Q})$ norms, respectively. In our analyses, we will assume
that $\mbox{\boldmath$f$}\in[L^{2}(Q)]^{d}$ and
$\mbox{\boldmath$u$}_{0}\in[H^{1}_{0}(\Omega)]^{d}$.

We will carry out a discrete-in-time analysis. Approximate or discrete values
of the variables at the different instants will be defined with sub-indices:
$\mbox{\boldmath$u$}_{n}$, for instance, denotes the velocity approximation at
the $n$-th time step. We will often use the identity

$\displaystyle
2(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$u$}_{n+1})=\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}\,,$
(10)

where
$\delta\mbox{\boldmath$u$}_{n+1}:=\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n}$.
A classical and useful fact is the skew-symmetry of the convective term. More
precisely, if $\mbox{\boldmath$v$}\in[H^{1}_{0}(\Omega)]^{d}$ is solenoidal,
there holds

$\displaystyle((\nabla\mbox{\boldmath$w$})\mbox{\boldmath$v$},\mbox{\boldmath$w$})=0\
\ \text{for all $\mbox{\boldmath$w$}\in[H^{1}_{0}(\Omega)]^{d}$.}$ (11)

If $v$ is not divergence-free, we will consider the following skew-symmetrised
form $\mathbf{c}$ of the convective derivative:

$\displaystyle\mathbf{c}(\mbox{\boldmath$v$},\mbox{\boldmath$w$}):=(\nabla\mbox{\boldmath$w$})\mbox{\boldmath$v$}+\frac{1}{2}(\nabla\cdot\mbox{\boldmath$v$})\mbox{\boldmath$w$}\,,$
(12)

which enjoys the property

$\displaystyle(\mathbf{c}(\mbox{\boldmath$v$},\mbox{\boldmath$w$}),\mbox{\boldmath$w$})=0\
\ \text{for all}\
\mbox{\boldmath$w$},\mbox{\boldmath$v$}\in[H^{1}_{0}(\Omega)]^{d}\,.$ (13)

For details on different forms of the convective term and their properties, we
refer the reader to, e.g., [20].

To construct IMEX schemes, we can combine backward differentiation formulas in
time with extrapolation rules of matching order of consistency. Our analysis
will be restricted to the first-order case, for which extrapolation means
simply replacing a certain quantity at $t_{n+1}$ by its previous value at
$t_{n}$, and

$\displaystyle\partial_{t}\mbox{\boldmath$u$}|_{t=t_{n+1}}\approx\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n})\,,$

where $\tau>0$ denotes the time-step size, which will be considered constant,
for simplicity – although all the results in this paper extend to variable
$\tau$ in a straightforward manner.

We will analyse both monolithic and fractional-step schemes. The stability
analysis will rely on different versions of the discrete Gronwall inequality.
The proof of the two results below can be found in [27, Lemma 5.1] (see also
[20, Lemma A.56]).

######  Lemma 2.1 (Discrete Gronwall inequality)

Let $N\in\mathbb{N}$, and $\alpha,B,a_{n},b_{n},c_{n}$ be non-negative numbers
for $n=1,\ldots,N$. Let us suppose that these numbers satisfy

$\displaystyle a_{N}+\sum_{n=1}^{N}b_{n}\leq
B+\sum_{n=1}^{N}c_{n}+\alpha\sum_{n=1}^{N-1}a_{n}\,.$ (14)

Then, the following inequality holds:

$\displaystyle a_{N}+\sum_{n=1}^{N}b_{n}\leq\mathrm{e}^{\alpha
N}\left(B+\sum_{n=1}^{N}c_{n}\right)\ \ \text{for}\ \ N\geq 1\,.$ (15)

######  Lemma 2.2 (Discrete Gronwall lemma, conditional version)

Let $N\in\mathbb{N}$, and
$\alpha,B,a_{n},b_{n},c_{n}$ be non-negative numbers for $n=1,\ldots,N$. Let
us suppose that these numbers satisfy

$\displaystyle a_{N}+\sum_{n=1}^{N}b_{n}\leq
B+\sum_{n=1}^{N}c_{n}+\alpha\sum_{n=1}^{N}a_{n}\,,$ (16)

with $\alpha<1$. Then, the following inequality holds:

$\displaystyle a_{N}+\sum_{n=1}^{N}b_{n}\leq\mathrm{e}^{\frac{\alpha
N}{1-\alpha}}\left(B+\sum_{n=1}^{N}c_{n}\right)\,.$ (17)

Notice that the only difference between (14) and (16) is where the last sum on
the right-hand side stops ($N-1$ or $N$).

##  3 Monolithic schemes

In this section we will study three different IMEX schemes. All of them share
the property that the velocity and the pressure are computed simultaneously at
each time step (hence the term “monolithic” to describe them). The proof of
their stability will be based on the discrete Gronwall inequality given in
Lemma 2.1. Hence, every time we refer to the discrete Gronwall inequality in
this section, we will be referring to Lemma 2.1.

###  3.1 Implicit stress-divergence formulation

We will start by considering the simplest possible IMEX scheme, where (other
than the convective velocity) only the viscosity is made explicit. In such a
case, the time-discrete momentum equation reads

$\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n})-\nabla\cdot\left(2\nu_{n}\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n+1}\right)+(\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n}+\nabla
p_{n+1}=\mbox{\boldmath$f$}_{n+1}\,.$ (18)

The standard variational problem for the $(n+1)$-th time step is given by:
find $(\mbox{\boldmath$u$}_{n+1},p_{n+1})\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$ such that

$\begin{split}\frac{1}{\tau}\left(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$}\right)+\left(2\nu_{n}\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n+1},\nabla^{\mathrm{s}}\mbox{\boldmath$w$}\right)+\big{(}(\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$}\big{)}-\left(p_{n+1},\nabla\cdot\mbox{\boldmath$w$}\right)&=\left(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$w$}\right),\\\
\left(q,\nabla\cdot\mbox{\boldmath$u$}_{n+1}\right)&=0\,,\end{split}$ (19)

for all $(\mbox{\boldmath$w$},q)\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$.

###### Remark 3.1

The decision to make the viscosity explicit in (19) stems from the fact that
in realistic situations (such as generalised Newtonian flows), the viscosity
is seldom known beforehand, being instead velocity-dependent. This is why we
consider the first-order extrapolation $\nu_{n+1}\approx\nu_{n}$.

The scheme (19) is rather classical, and well-understood. Nevertheless, we now
reproduce its stability result (and its proof) for completeness, and to
highlight the main differences that appear in more involved schemes presented
later.

######  Lemma 3.1 (Stability of the implicit stress-divergence formulation)

Let us suppose that the time-step size is given by $\tau=T/N$, $N\geq 2$.
Then, the scheme (19) satisfies the following stability inequality:

$\|\mbox{\boldmath$u$}_{N}\|^{2}+4\sum_{n=1}^{N}\tau\big{\|}\sqrt{\nu_{n-1}}\,\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n}\big{\|}^{2}\leq\mathrm{e}^{\varepsilon}\Bigg{[}\left(1+\varepsilon\frac{\tau}{T}\right)\|\mbox{\boldmath$u$}_{0}\|^{2}+\left(\tau+\frac{T}{\varepsilon}\right)\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\Bigg{]}\,,$
(20)

where $\varepsilon>0$ is an arbitrary constant.

Proof. Setting $\mbox{\boldmath$w$}=2\tau\mbox{\boldmath$u$}_{n+1}$ and
$q=2\tau p_{n+1}$ in (19) and adding the resulting equations yields

$\displaystyle
2(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$u$}_{n+1})+4\tau\big{\|}\sqrt{\nu_{n}}\,\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}+2\tau((\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$u$}_{n+1})=2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})\,.$

The convective term vanishes due to its skew-symmetric property (11). Applying
(10) and the Cauchy-Schwarz and Young inequalities, we obtain

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+4\tau\big{\|}\sqrt{\nu_{n}}\,\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}$
$\displaystyle\leq 2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})$
$\displaystyle=2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n}+\delta\mbox{\boldmath$u$}_{n+1})$
$\displaystyle\leq
2\tau\|\delta\mbox{\boldmath$u$}\|\,\|\mbox{\boldmath$f$}_{n+1}\|+2\tau\|\mbox{\boldmath$f$}_{n+1}\|\,\|\mbox{\boldmath$u$}_{n}\|$
$\displaystyle\leq\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}+\Big{(}\tau+\frac{T}{\varepsilon}\Big{)}\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}+\tau\frac{\varepsilon}{T}\|\mbox{\boldmath$u$}_{n}\|^{2},$

where $\varepsilon>0$ is arbitrary. Rearranging terms we get

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+4\tau\big{\|}\sqrt{\nu_{n}}\,\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}\leq\varepsilon\frac{\tau}{T}\|\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\tau^{2}+\frac{\tau
T}{\varepsilon}\Big{)}\|\mbox{\boldmath$f$}_{n+1}\|^{2}\,.$

Adding this equality from $n=0$ to $n=N-1$ yields

$\displaystyle\|\mbox{\boldmath$u$}_{N}\|^{2}+4\tau\sum_{n=1}^{N}\big{\|}\sqrt{\nu_{n-1}}\,\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n}\big{\|}^{2}\leq\left(1+\varepsilon\frac{\tau}{T}\right)\|\mbox{\boldmath$u$}_{0}\|^{2}+\varepsilon\frac{\tau}{T}\sum_{n=1}^{N-1}\|\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\tau+\frac{T}{\varepsilon}\Big{)}\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2},$

and the result follows as a direct application of the discrete Gronwall lemma.

###### Remark 3.2

The above result applied to $\mbox{\boldmath$f$}=\mbox{\boldmath$0$}$ leads to
the fact that (19) is strongly stability preserving, i.e.,

$\|\mbox{\boldmath$u$}_{N}\|^{2}+4\sum_{n=1}^{N}\tau\big{\|}\sqrt{\nu_{n-1}}\,\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n}\big{\|}^{2}\leq\|\mbox{\boldmath$u$}_{0}\|^{2}$

The main question, to be addressed next, is whether and how we can preserve
stability even if a part of the viscous term is treated explicitly.

###  3.2 Implicit-explicit stress-divergence formulations

The simplest, naive way to decouple the velocity components is to split the
viscous term into two parts, as in (6), and make the second term explicit.
This choice leads to the following semi-discrete scheme:

$\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n})-\nabla\cdot\left(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1}\right)+(\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n}+\nabla
p_{n+1}=\nabla\cdot\left(\nu_{n}\nabla^{\top}\mbox{\boldmath$u$}_{n}\right)+\mbox{\boldmath$f$}_{n+1}\,,$
(21)

with the corresponding weak form: find
$(\mbox{\boldmath$u$}_{n+1},p_{n+1})\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$ such that

$\begin{split}\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+\left(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1},\nabla\mbox{\boldmath$w$}\right)+\big{(}(\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$}\big{)}-\left(p_{n+1},\nabla\cdot\mbox{\boldmath$w$}\right)&=\\\
-\left(\nu_{n}\nabla^{\top}\mbox{\boldmath$u$}_{n},\nabla\mbox{\boldmath$w$}\right)+(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$w$})\,,&\\\
\left(q,\nabla\cdot\mbox{\boldmath$u$}_{n+1}\right)&=0\,,\end{split}$

for all $(\mbox{\boldmath$w$},q)\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$.

Testing with the same functions as in the proof of the last result, we get to

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}+2\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}+\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}=\|\mbox{\boldmath$u$}_{n}\|^{2}-2\tau(\nu_{n}\nabla^{\top}\mbox{\boldmath$u$}_{n},\nabla\mbox{\boldmath$u$}_{n+1})+2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})\,.$

In this case, an extra term (that requires bounding) appears on the right-hand
side:

$-2(\nu_{n}\nabla^{\top}\mbox{\boldmath$u$}_{n},\nabla\mbox{\boldmath$u$}_{n+1})\leq
2\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}\,\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}\leq\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}\,,$

where we have used the Cauchy-Schwarz and Young inequalities. Then, estimating
the forcing term as done for the implicit case, we get the following bound at
step $(n+1)$:

$\|\mbox{\boldmath$u$}_{n+1}\|^{2}+\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}\leq\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\frac{\tau}{T}\|\mbox{\boldmath$u$}_{n}\|^{2}+(\tau^{2}+\tau
T)\|\mbox{\boldmath$f$}_{n+1}\|^{2}\,.$ (22)

At this point, there is no guarantee that the viscous term on the left-hand
side will be able to control the one on the right-hand side, since the
viscosity on the left-hand side term is $\nu_{n}$, not $\nu_{n+1}$. Some
stability might still be attained if $\nu_{n+1}\leq\nu_{n}$, or if $\tau$ is
small enough so that $\nu_{n+1}\approx\nu_{n}$; otherwise, the scheme could
become unstable, especially for rapidly varying viscosity fields. This may
explain the poor stability observed for a related implicit-explicit SD
formulation in [3], see also Section 7.2 below.

According to the analysis above, guaranteeing stability of the IMEX-SD
formulation would require knowing $\nu_{n+1}$ already at the $(n+1)$-th time-
step. Although this is seldom realistic, a scheme that requires $\nu_{n+1}$
can still make practical sense if combined, for instance, with a Picard-like
method linearising the current viscosity. Therefore, let us now present a
scheme that assumes $\nu_{n+1}$ as already known at the point of computing
$(\mbox{\boldmath$u$}_{n+1},p_{n+1})$. Consider the first-order extrapolation

$\displaystyle\nu_{n+1}\nabla^{\top}\mbox{\boldmath$u$}_{n+1}=\sqrt{\nu_{n+1}}\sqrt{\nu_{n+1}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n+1}\approx\sqrt{\nu_{n+1}}\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n}\,,$

from which we can propose the following alternative to (21):

$\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n})-\nabla\cdot\left(\nu_{n+1}\nabla\mbox{\boldmath$u$}_{n+1}\right)+(\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n}+\nabla
p_{n+1}=\nabla\cdot\left(\sqrt{\nu_{n+1}}\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n}\right)+\mbox{\boldmath$f$}_{n+1}$

and the corresponding weak form: find
$(\mbox{\boldmath$u$}_{n+1},p_{n+1})\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$ such that

$\begin{split}\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+\left(\nu_{n+1}\nabla\mbox{\boldmath$u$}_{n+1},\nabla\mbox{\boldmath$w$}\right)+((\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})-\left(p_{n+1},\nabla\cdot\mbox{\boldmath$w$}\right)&=\\\
-\left(\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n},\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$w$}\right)+(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$w$})\,,&\\\
(q,\nabla\cdot\mbox{\boldmath$u$}_{n+1})&=0\,,\end{split}$ (23)

for all $(\mbox{\boldmath$w$},q)\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$. Based on that, we will prove the following stability
estimate.

######  Theorem 3.2 (Stability of an IMEX stress-divergence formulation)

Let us suppose that the time-step size is given by $\tau=T/N$, $N\geq 2$.
Then, the scheme (23) satisfies the following stability inequality

$\|\mbox{\boldmath$u$}_{N}\|^{2}+\tau\big{\|}\sqrt{\nu_{N}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}+\varepsilon\sum_{n=1}^{N-1}\frac{\tau^{2}}{T}\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}\leq\mathrm{e}^{\varepsilon}\Bigg{[}C_{0}+\Big{(}\tau+\frac{T}{\varepsilon}\Big{)}\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\Bigg{]},$
(24)

where $\varepsilon>0$ is an arbitrary constant, and

$\displaystyle
C_{0}=\left(1+\varepsilon\frac{\tau}{T}\right)\|\mbox{\boldmath$u$}_{0}\|^{2}+\tau\|\sqrt{\nu_{0}}\,\nabla\mbox{\boldmath$u$}_{0}\|^{2}\,.$

Proof. Using the same techniques as in the proof of (22) we can show that

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau\big{\|}\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}\leq\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\varepsilon\frac{\tau}{T}\|\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\tau+\frac{T}{\varepsilon}\Big{)}\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}.$

Adding up from $n=0$ to $n=N-1$ gives

$\displaystyle\|\mbox{\boldmath$u$}_{N}\|^{2}+\tau\big{\|}\sqrt{\nu_{N}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}\leq$
$\displaystyle\left(1+\varepsilon\frac{\tau}{T}\right)\|\mbox{\boldmath$u$}_{0}\|^{2}+\tau\big{\|}\sqrt{\nu_{0}}\,\nabla\mbox{\boldmath$u$}_{0}\big{\|}^{2}+\varepsilon\frac{\tau}{T}\sum_{n=1}^{N-1}\|\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\tau+\frac{T}{\varepsilon}\Big{)}\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\,.$
(25)

The proof is completed by adding

$\displaystyle\varepsilon\frac{\tau}{T}\sum_{n=1}^{N-1}\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}$

to both sides of (25) and using the discrete Gronwall inequality with
$a_{n}=\|\mbox{\boldmath$u$}_{N}\|^{2}+\tau\big{\|}\sqrt{\nu_{N}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}$,
$b_{n}=\varepsilon\frac{\tau}{T}\tau\big{\|}\sqrt{\nu_{N}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}$,
and $c_{n}$ gathering the terms involving $\mbox{\boldmath$u$}_{0}$ and $f$ on
the right-hand side of (25).

###  3.3 Implicit-explicit generalised Laplacian formulation

Due to the squared time-step size accompanying the viscous sum in (24), the
stability estimate proven in Theorem 3.2 is somewhat weaker than the one
proved for (19). It is, still, an unconditional stability attained despite the
explicit treatment of the coupling term. Of course, this comes at the cost of
requiring a-priori knowledge of the viscosity field. To circumvent this
requirement, we now present a different IMEX scheme based on replacing the
viscous term by the generalised Laplacian version (9). In that case, treating
the coupling term explicitly will lead to the semi-discrete momentum equation

$\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n})-\nabla\cdot\left(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1}\right)+(\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n}+\nabla
p_{n+1}=\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n}+\mbox{\boldmath$f$}_{n+1}\,.$

The corresponding variational problem is: find
$(\mbox{\boldmath$u$}_{n+1},p_{n+1})\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$ such that

$\begin{split}\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1},\nabla\mbox{\boldmath$w$})+((\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})-(p_{n+1},\nabla\cdot\mbox{\boldmath$w$})&=\\\
(\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n},\mbox{\boldmath$w$})+(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$w$})\,,&\\\
(q,\nabla\cdot\mbox{\boldmath$u$}_{n+1})&=0\,,\end{split}$ (26)

for all $(\mbox{\boldmath$w$},q)\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega)$. The stability of (26) is stated next.

######  Theorem 3.3 (Stability of the IMEX generalised Laplacian formulation)

Let us suppose that $\nabla\nu\in[L^{\infty}(\bar{Q})]^{d}$, and that the
time-step size is given by $\tau=T/N$ with $N\geq 2$. Then, the IMEX scheme
(26) satisfies the following stability inequality

$\|\mbox{\boldmath$u$}_{N}\|^{2}+\sum_{n=1}^{N}\tau\big{\|}\sqrt{\nu_{n-1}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}\leq
M\Bigg{[}C_{0}+\left(2\tau+\frac{T}{\varepsilon}\right)\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\Bigg{]},$
(27)

where

$\displaystyle
M=\mathrm{exp}\Big{(}\varepsilon+\frac{2T\|\nabla\nu\|_{\infty}^{2}}{\nu_{\mathrm{min}}}\Big{)}\,,\quad
C_{0}=\Big{(}1+\frac{\varepsilon\tau}{T}+\frac{2\|\nabla\nu\|^{2}_{\infty}\tau}{\nu_{\mathrm{min}}}\Big{)}(\|\mbox{\boldmath$u$}_{0}\|^{2}+\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{0}\|^{2})\,,$

and $\varepsilon>0$ is an arbitrary constant.

Proof. Set $\mbox{\boldmath$w$}=2\tau\mbox{\boldmath$u$}_{n+1}$ and $q=2\tau
p_{n+1}$ in (26) to get

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}+2\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}$
$\displaystyle=2\tau(\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n},\mbox{\boldmath$u$}_{n+1})+2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})\,.$
(28)

We now bound the terms on the right-hand side of the inequality above. Using
the Cauchy-Schwarz, Hölder, and Young inequalities, we get

$\displaystyle
2\tau(\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n},\mbox{\boldmath$u$}_{n+1})$
$\displaystyle=2\tau(\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n},\delta\mbox{\boldmath$u$}_{n+1}+\mbox{\boldmath$u$}_{n})$
$\displaystyle\leq
2\tau\|\delta\mbox{\boldmath$u$}_{n+1}\|\,\|\nabla\nu\|_{\infty}\|\nabla\mbox{\boldmath$u$}_{n}\|+2\tau\|\nabla\nu\|_{\infty}\|\nabla\mbox{\boldmath$u$}_{n}\|\,\|\mbox{\boldmath$u$}_{n}\|$
$\displaystyle\leq\frac{1}{2}\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}+2\tau^{2}\|\nabla\nu\|_{\infty}^{2}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\frac{\tau\|\nabla\nu\|_{\infty}^{2}}{\nu_{\mathrm{min}}}\|\mbox{\boldmath$u$}_{n}\|^{2}\,,$
(29)

and

$\displaystyle 2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})$
$\displaystyle=2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n}+\delta\mbox{\boldmath$u$}_{n+1})$
$\displaystyle\leq
2\tau\|\delta\mbox{\boldmath$u$}_{n+1}\|\,\|\mbox{\boldmath$f$}_{n+1}\|+2\tau\|\mbox{\boldmath$f$}_{n+1}\|\,\|\mbox{\boldmath$u$}_{n}\|$
$\displaystyle\leq\frac{1}{2}\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}+\left(2\tau+\frac{T}{\varepsilon}\right)\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}+\frac{\varepsilon\tau}{T}\|\mbox{\boldmath$u$}_{n}\|^{2}\,.$
(30)

Inserting (29) and (30) in (28) gives

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+2\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}$
$\displaystyle\leq\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\frac{\varepsilon\tau}{T}+\frac{\tau\|\nabla\nu\|_{\infty}^{2}}{\nu_{\mathrm{min}}}\Big{)}\|\mbox{\boldmath$u$}_{n}\|^{2}+\frac{2\tau\|\nabla\nu\|_{\infty}^{2}}{\nu_{\mathrm{min}}}\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\left(2\tau+\frac{T}{\varepsilon}\right)\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}$
$\displaystyle\leq\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\frac{\varepsilon\tau}{T}+\frac{2\tau\|\nabla\nu\|_{\infty}^{2}}{\tau\nu_{\mathrm{min}}}\Big{)}(\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2})+\left(2\tau+\frac{T}{\varepsilon}\right)\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}\,.$

Rearranging terms leads to

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau\Big{(}\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}-\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}\Big{)}+\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}$
$\displaystyle\leq$
$\displaystyle\left(\frac{\varepsilon\tau}{T}+\frac{2\tau\|\nabla\nu\|_{\infty}^{2}}{\nu_{\mathrm{min}}}\right)(\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2})+\left(2\tau+\frac{T}{\varepsilon}\right)\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}\,.$

Adding up from $n=0$ to $n=N-1$ gives

$\displaystyle\|\mbox{\boldmath$u$}_{N}\|^{2}+\tau\sum_{n=0}^{N-1}\Big{(}\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}-\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}\Big{)}+\sum_{n=1}^{N}\tau\big{\|}\sqrt{\nu_{n-1}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}$
$\displaystyle\leq\ $
$\displaystyle\|\mbox{\boldmath$u$}_{0}\|^{2}+\Big{(}\frac{\varepsilon\tau}{T}+\frac{2\tau\|\nabla\nu\|_{\infty}^{2}}{\nu_{\mathrm{min}}}\Big{)}\sum_{n=0}^{N-1}\left(\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}\right)+\left(2\tau+\frac{T}{\varepsilon}\right)\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\,.$

Finally, using

$\displaystyle\sum_{n=0}^{N-1}\Big{(}\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}-\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}\Big{)}\geq\big{\|}\sqrt{\nu_{N-1}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}-\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{0}\|^{2}\,,$

estimate (27) follows directly from the discrete Gronwall lemma.

With the stability result just proven, we eliminate the need to know the
viscosity at time step $n+1$ and still have an unconditionally stable scheme.
The price to pay is a higher regularity requirement on the viscosity field,
which must now be Lipschitz continuous in space. This higher requirement may
not be satisfied in the fully discrete case for some practical scenarios, but
in Section 6.2 we will discuss possible remedies.

##  4 Fractional-step schemes

This section is devoted to extending the IMEX methods from Section 3 to a
class of fractional-step schemes: the incremental pressure correction method.
We shall consider this prototypical scheme in its standard variant, which
allows for variable viscosity [28].

###  4.1 IMEX stress-divergence formulation

For the Navier–Stokes equations in stress-divergence form, consider the
following first-order IMEX fractional-step scheme: find the solution
$(\mbox{\boldmath$u$}_{n+1},\hat{\mbox{\boldmath$u$}}_{n+1},p_{n+1})$ of the
system

$\displaystyle\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\hat{\mbox{\boldmath$u$}}_{n})-\nabla\cdot(\nu_{n+1}\nabla\mbox{\boldmath$u$}_{n+1})+\mathbf{c}(\mbox{\boldmath$u$}_{n},\mbox{\boldmath$u$}_{n+1})$
$\displaystyle=\nabla\cdot(\sqrt{\nu_{n+1}}\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n})-\nabla
p_{n}+\mbox{\boldmath$f$}_{n+1}\,,$
$\displaystyle\frac{1}{\tau}(\hat{\mbox{\boldmath$u$}}_{n+1}-\mbox{\boldmath$u$}_{n+1})+\nabla(p_{n+1}-p_{n})$
$\displaystyle=\mbox{\boldmath$0$}\,,$
$\displaystyle\nabla\cdot\hat{\mbox{\boldmath$u$}}_{n+1}$ $\displaystyle=0\,,$

with boundary conditions $\mbox{\boldmath$u$}_{n+1}=\mbox{\boldmath$0$}$ and
$\mbox{\boldmath$n$}\cdot\hat{\mbox{\boldmath$u$}}_{n+1}=0$ on $\Gamma$. Since
$u$ is no longer divergence-free, the skew-symmetrised form $\mathbf{c}$ of
the convective term (12) is used. In addition, since the initial condition
$\mbox{\boldmath$u$}_{0}$ is divergence-free, we assume $p_{0}=0$ for the rest
of this section.

We can use simple algebraic manipulation to eliminate the end-of-step velocity
$\hat{\mbox{\boldmath$u$}}$ from the system, reformulating the
$(\hat{\mbox{\boldmath$u$}},p)$-subproblem as a pressure Poisson equation. The
resulting weak formulation is: find
$(\mbox{\boldmath$u$}_{n+1},p_{n+1})\in[H^{1}_{0}(\Omega)]^{d}\times[L^{2}_{0}(\Omega)\cap
H^{1}(\Omega)]$ such that

$\begin{split}\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+(\nu_{n+1}\nabla\mbox{\boldmath$u$}_{n+1},\nabla\mbox{\boldmath$w$})+((\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+\frac{1}{2}((\nabla\cdot\mbox{\boldmath$u$}_{n})\mbox{\boldmath$u$}_{n+1},\mbox{\boldmath$w$})&=\\\
(2p_{n}-p_{n-1},\nabla\cdot\mbox{\boldmath$w$})-\left(\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n},\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$w$}\right)+(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$w$})\,,&\\\
\tau(\nabla(p_{n+1}-p_{n}),\nabla
q)+\left(q,\nabla\cdot\mbox{\boldmath$u$}_{n+1}\right)&=0\,,\end{split}$ (31)

for all
$(\mbox{\boldmath$w$},q)\in[H^{1}_{0}(\Omega)]^{d}\times[L^{2}_{0}(\Omega)\cap
H^{1}(\Omega)]$. This split-step scheme is also unconditionally stable, as the
following result shows.

######  Theorem 4.1 (Stability of an IMEX fractional-step scheme in SD form)

Let the time step be given by $\tau=T/N$ with $N\geq 2$. Then, the fractional-
step method (31) satisfies the following stability estimate

$\begin{split}\|\mbox{\boldmath$u$}_{N}\|^{2}+\tau\big{\|}\sqrt{\nu_{N}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}+\tau^{2}\|\nabla
p_{N}\|^{2}+\frac{\tau}{T}\sum_{n=1}^{N}\Big{(}\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\tau^{2}\|\nabla
p_{n}\|^{2}\Big{)}\\\ \leq
M\Bigg{[}\|\mbox{\boldmath$u$}_{0}\|^{2}+\tau\|\sqrt{\nu_{0}}\,\nabla\mbox{\boldmath$u$}_{0}\|^{2}+T\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\Bigg{]}\,,\end{split}$
(32)

where

$\displaystyle M$
$\displaystyle=\mathrm{exp}\Bigg{(}\frac{1}{1-\tau/T}\Bigg{)}=\mathrm{exp}\Bigg{(}\frac{N}{N-1}\Bigg{)}\leq\mathrm{e}^{2}\,.$

Proof. Setting $\mbox{\boldmath$w$}=2\tau\mbox{\boldmath$u$}_{n+1}$ in (31)
yields

$\begin{split}&\|\mbox{\boldmath$u$}_{n+1}\|^{2}+\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}-\|\mbox{\boldmath$u$}_{n}\|^{2}+2\tau\big{\|}\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}\\\
&=2\tau\left[\left(2p_{n}-p_{n-1},\nabla\cdot\mbox{\boldmath$u$}_{n+1}\right)-\left(\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n},\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$u$}_{n+1}\right)+(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})\right]\\\
&\leq
2\tau\left[\left(2p_{n}-p_{n-1},\nabla\cdot\mbox{\boldmath$u$}_{n+1}\right)+\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\big{\|}\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}+(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})\right].\end{split}$
(33)

Compared to the monolithic cases presented in the last section, we have an
additional pressure term to estimate. Denoting $\delta
p_{n+1}:=p_{n+1}-p_{n}$, the second equation in (31) implies

$\tau(\nabla\delta p_{n+1}-\nabla\delta p_{n},\nabla
q)=-(\nabla\cdot\delta\mbox{\boldmath$u$}_{n+1},q)=\left(\delta\mbox{\boldmath$u$}_{n+1},\nabla
q\right)\,.$ (34)

Now, setting $q=2\tau(2p_{n}-p_{n-1})$ in (31) gives

$\displaystyle 2\tau(2p_{n}-p_{n-1},\nabla\cdot\mbox{\boldmath$u$}_{n+1})$
$\displaystyle=-2\tau^{2}(\nabla(2p_{n}-p_{n-1}),\nabla(p_{n+1}-p_{n}))$
$\displaystyle\equiv-2\tau^{2}\big{(}\nabla[p_{n+1}-(\delta p_{n+1}-\delta
p_{n})],\nabla(p_{n+1}-p_{n})\big{)}$ $\displaystyle=-2\tau^{2}(\nabla
p_{n+1},\nabla p_{n+1}-\nabla p_{n})+2\tau^{2}(\nabla\delta
p_{n+1}-\nabla\delta p_{n},\nabla\delta p_{n+1}).$

For the first term on the right-hand side above we have

$\displaystyle-2(\nabla p_{n+1},\nabla p_{n+1}-\nabla p_{n})=\|\nabla
p_{n}\|^{2}-\|\nabla p_{n+1}\|^{2}-\|\nabla\delta p_{n+1}\|^{2}\,,$

and for the second one we set $q=\tau\delta p_{n+1}$ in (34), so that

$\displaystyle 2\tau^{2}(\nabla\delta p_{n+1}-\nabla\delta p_{n},\nabla\delta
p_{n+1})=2\tau(\delta\mbox{\boldmath$u$}_{n+1},\nabla\delta
p_{n+1})\leq\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}+\tau^{2}\|\nabla\delta
p_{n+1}\|^{2}\,.$

Collecting these last results, we estimate the pressure term in (33) as

$\displaystyle
2\tau(2p_{n}-p_{n-1},\nabla\cdot\mbox{\boldmath$u$}_{n+1})\leq\tau^{2}\|\nabla
p_{n}\|^{2}-\tau^{2}\|\nabla
p_{n+1}\|^{2}+\|\delta\mbox{\boldmath$u$}_{n+1}\|^{2}\,.$ (35)

Moreover,

$\displaystyle
2\tau(\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})\leq\frac{\tau}{T}\|\mbox{\boldmath$u$}_{n+1}\|^{2}+T\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}\,.$
(36)

Replacing (35) and (36) in (33) we get

$\displaystyle\left(\|\mbox{\boldmath$u$}_{n+1}\|^{2}+\big{\|}\sqrt{\nu_{n+1}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}+\tau^{2}\|\nabla
p_{n+1}\|^{2}\right)-\left(\|\mbox{\boldmath$u$}_{n}\|^{2}+\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\tau^{2}\|\nabla
p_{n}\|^{2}\right)$
$\displaystyle\leq\frac{\tau}{T}\|\mbox{\boldmath$u$}_{n+1}\|^{2}+T\tau\|\mbox{\boldmath$f$}_{n+1}\|^{2}\,,$

and adding from $n=1$ to $n=N-1$ leads to

$\displaystyle\|\mbox{\boldmath$u$}_{N}\|^{2}+\tau\big{\|}\sqrt{\nu_{N}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}+\tau^{2}\|\nabla
p_{N}\|^{2}$
$\displaystyle\leq\|\mbox{\boldmath$u$}_{0}\|^{2}+\tau\big{\|}\sqrt{\nu_{0}}\,\nabla\mbox{\boldmath$u$}_{0}\big{\|}^{2}+\tau^{2}\|\nabla
p_{0}\|^{2}+\frac{\tau}{T}\sum_{n=1}^{N}\|\mbox{\boldmath$u$}_{n}\|^{2}+T\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\,.$

Finally, adding

$\displaystyle\frac{\tau}{T}\sum_{n=1}^{N}\Big{(}\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\tau^{2}\|\nabla
p_{n}\|^{2}\Big{)}$

to both sides of the last inequality, and noticing that $\tau/T<1$, we
conclude by using the Gronwall inequality from Lemma 2.2.

As it was the case with Theorem 3.2, the stability result just shown is in a
weak norm, and also requires the knowledge of $\nu_{n+1}$. To remedy this, we
will present below an IMEX fractional-step scheme based on the generalised
Laplacian rewriting of the diffusive term.

###  4.2 IMEX generalised Laplacian formulation

After elimination of the end-of-step velocity, the weak problem for the GL
version of the IMEX incremental projection scheme is: find
$(\mbox{\boldmath$u$}_{n+1},p_{n+1})\in[H^{1}_{0}(\Omega)]^{d}\times[L^{2}_{0}(\Omega)\cap
H^{1}(\Omega)]$ such that

$\begin{split}\frac{1}{\tau}(\mbox{\boldmath$u$}_{n+1}-\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+\left(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1},\nabla\mbox{\boldmath$w$}\right)+((\nabla\mbox{\boldmath$u$}_{n+1})\mbox{\boldmath$u$}_{n},\mbox{\boldmath$w$})+\frac{1}{2}((\nabla\cdot\mbox{\boldmath$u$}_{n})\mbox{\boldmath$u$}_{n+1},\mbox{\boldmath$w$})&=\\\
\left(2p_{n}-p_{n-1},\nabla\cdot\mbox{\boldmath$w$}\right)+\left(\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n},\mbox{\boldmath$w$}\right)+\left(\mbox{\boldmath$w$},\mbox{\boldmath$f$}_{n+1}\right)\,,&\\\
\tau(\nabla(p_{n+1}-p_{n}),\nabla
q)+\left(\nabla\cdot\mbox{\boldmath$u$}_{n+1},q\right)&=0\,,\end{split}$ (37)

for all
$(\mbox{\boldmath$w$},q)\in[H^{1}_{0}(\Omega)]^{d}\times[L^{2}_{0}(\Omega)\cap
H^{1}(\Omega)]$. The following result states the stability for this scheme.

######  Theorem 4.2 (Stability of the IMEX fractional-step scheme in GL form)

Let
us assume that $\nabla\nu\in[L^{\infty}(\bar{Q}]^{d}$. Then, the fractional-
step IMEX-GL scheme (37) yields the stability estimate

$\begin{split}\|\mbox{\boldmath$u$}_{N}\|^{2}+\gamma\tau\big{\|}\sqrt{\nu_{N-1}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}+\tau^{2}\|\nabla
p_{N}\|^{2}+\sum_{n=1}^{N}\tau\Big{[}(2-\gamma)\big{\|}\sqrt{\nu_{n-1}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}+\alpha\tau^{2}\|\nabla
p_{n}\|^{2}\Big{]}\\\ \leq\mathrm{exp}\left(\frac{\alpha
T}{1-\alpha\tau}\right)\left[\|\mbox{\boldmath$u$}_{0}\|^{2}+\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{0}\|^{2}+\varepsilon
T\sum_{n=1}^{N}\tau\|\mbox{\boldmath$f$}_{n}\|^{2}\right]\,,\end{split}$ (38)

for any positive time-step size $\tau$ satisfying

$\tau<\left(\frac{\|\nabla\nu\|_{\infty}^{2}}{\gamma\nu_{\mathrm{min}}}+\frac{1}{\varepsilon
T}\right)^{-1}=:\alpha^{-1}\,,$ (39)

with $\gamma\in(0,2)$ and $\varepsilon>0$, both chosen arbitrarily.

Proof. Compared to the stress-divergence version, the only term estimated
differently is

$\displaystyle
2\tau(\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n}+\mbox{\boldmath$f$}_{n+1},\mbox{\boldmath$u$}_{n+1})$
$\displaystyle\leq
2\tau\|\nabla\nu\|_{\infty}\|\nabla\mbox{\boldmath$u$}_{n}\|\,\|\mbox{\boldmath$u$}_{n+1}\|+2\tau\|\mbox{\boldmath$f$}_{n+1}\|\,\|\mbox{\boldmath$u$}_{n+1}\|$
$\displaystyle\leq\gamma\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\frac{\tau\|\nabla\nu\|_{\infty}^{2}}{\gamma\nu_{\mathrm{min}}}+\frac{\tau}{\varepsilon
T}\Big{)}\|\mbox{\boldmath$u$}_{n+1}\|^{2}+\tau\varepsilon
T\|\mbox{\boldmath$f$}_{n+1}\|^{2},$

in which $0<\gamma<2$ and $\varepsilon>0$. Hence,

$\displaystyle\|\mbox{\boldmath$u$}_{n+1}\|^{2}+\tau^{2}\|\nabla
p_{n+1}\|^{2}+2\tau\big{\|}\sqrt{\nu_{n}}\,\nabla\mbox{\boldmath$u$}_{n+1}\big{\|}^{2}$
$\displaystyle\leq\|\mbox{\boldmath$u$}_{n}\|^{2}+\tau^{2}\|\nabla
p_{n}\|^{2}+\gamma\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{n}\|^{2}+\Big{(}\frac{\tau\|\nabla\nu\|_{\infty}^{2}}{\gamma\nu_{\mathrm{min}}}+\frac{\tau}{\varepsilon
T}\Big{)}\|\mbox{\boldmath$u$}_{n+1}\|^{2}+\tau\varepsilon
T\|\mbox{\boldmath$f$}_{n+1}\|^{2}.$

Adding up from $n=0$ to $n=N-1$ and using that $p_{0}=0$ yields

$\begin{split}&\|\mbox{\boldmath$u$}_{N}\|^{2}+\gamma\tau\big{\|}\sqrt{\nu_{N-1}}\,\nabla\mbox{\boldmath$u$}_{N}\big{\|}^{2}+\tau^{2}\|\nabla
p_{N}\|^{2}+(2-\gamma)\tau\sum_{n=1}^{N}\big{\|}\sqrt{\nu_{n-1}}\,\nabla\mbox{\boldmath$u$}_{n}\big{\|}^{2}\\\
\leq\
&\|\mbox{\boldmath$u$}_{0}\|^{2}+\gamma\tau\nu_{\mathrm{min}}\|\nabla\mbox{\boldmath$u$}_{0}\|^{2}+\tau\varepsilon
T\sum_{n=1}^{N}\|\mbox{\boldmath$f$}_{n}\|^{2}+\Big{(}\frac{\tau\|\nabla\nu\|_{\infty}^{2}}{\gamma\nu_{\mathrm{min}}}+\frac{\tau}{\varepsilon
T}\Big{)}\sum_{n=1}^{N}\|\mbox{\boldmath$u$}_{n}\|^{2}\,.\end{split}$ (40)

Finally, by adding

$\displaystyle\Big{(}\frac{\tau\|\nabla\nu\|_{\infty}^{2}}{\gamma\nu_{\mathrm{min}}}+\frac{\tau}{\varepsilon
T}\Big{)}\sum_{n=1}^{N}\tau^{2}\|\nabla p_{n}\|^{2}$

to both sides of (40) and choosing $\tau$ according to (39), we get the
stability estimate (38) by applying Lemma 2.2.

###### Remark 4.1

The last stability result requires a CFL condition on the time-step size
$\tau$. It is worth noticing that this condition includes only the viscosity
field (i.e., the problem data). In addition, since the constants $\gamma$ and
$1/\varepsilon$ in (39) can be chosen arbitrarily close to $2$ and $0$,
respectively, this restriction can be simplified to

$\tau<\frac{2\nu_{\mathrm{min}}}{\|\nabla\nu\|_{\infty}^{2}}\,,$ (41)

which can also be obtained by letting $T\rightarrow\infty$, that is, letting
the problem tend to a steady-state. Interestingly, (41) is related to the
condition imposed on the reaction term in [6] to analyse a mixed finite
element method for the linear, stationary Oseen equation.

##  5 Summary of theoretical results

Let us briefly collect the results of our analysis. We have proved the
stability of several IMEX discretisations, and the corresponding results are
summarised on Table 1. For four discretisations we have unconditional
stability (albeit in different norms, some stronger than others), and for
another one we have proved stability under a time-step restriction (41). The
only method where stability has not been proven is the SD formulation with
fully explicit viscosity, that is, (21). Therefore, the main purpose of our
numerical examples will be to verify whether our analysis is sharp with
respect to condition (41) and to the (lack of) stability of the “naive” IMEX-
SD discretisation.

Table 1: Theoretical results on the stability of different IMEX treatments of the viscous term. Type | Left-hand side term | Right-hand side term | Stability
---|---|---|---
Monolithic | $-\nabla\cdot(2\nu_{n}\nabla^{\mathrm{s}}\mbox{\boldmath$u$}_{n+1})$ |  | unconditional
Monolithic | $-\nabla\cdot\big{(}\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1})$ | $\nabla\cdot(\nu_{n}\nabla^{\top}\mbox{\boldmath$u$}_{n}\big{)}$ | not proved
Monolithic | $-\nabla\cdot(\nu_{n+1}\nabla\mbox{\boldmath$u$}_{n+1})$ | $\nabla\cdot\big{(}\sqrt{\nu_{n+1}}\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n}\big{)}$ | unconditional
Monolithic | $-\nabla\cdot(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1})$ | $\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n}$ | unconditional
Fractional | $-\nabla\cdot(\nu_{n+1}\nabla\mbox{\boldmath$u$}_{n+1})$ | $\nabla\cdot\big{(}\sqrt{\nu_{n+1}}\sqrt{\nu_{n}}\,\nabla^{\top}\mbox{\boldmath$u$}_{n}\big{)}$ | unconditional
Fractional | $-\nabla\cdot(\nu_{n}\nabla\mbox{\boldmath$u$}_{n+1})$ | $\nabla^{\top}\mbox{\boldmath$u$}_{n}\nabla\nu_{n}$ | $\tau<\tau_{\mathrm{max}}(\nu)$

##  6 On the fully discrete problem

Although we restrict our stability analysis to the semi-discrete setting, some
brief comments should be drawn on the fully discrete problem. We will restrict
our discussion to conforming finite element spaces with continuous pressures,
since that choice does not modify the matrix structure of the semi-discrete
problems.

###  6.1 Matrix structure of the IMEX methods

One of the most appealing computational features of the IMEX schemes presented
herein is the resulting matrix structure. Consider a generic spatial
discretisation applied to $\mbox{\boldmath$u$}_{n}$ and $p_{n}$, denoting by
$\mbox{\boldmath$U$}_{n}$ and $\mbox{\boldmath$P$}_{n}$ the respective vectors
of degrees of freedom. Since we are interested in problems where the viscosity
may not be known, and also interested in IMEX schemes leading to solving only
linear systems at each time step, we shall only discuss the schemes that do
not require $\nu_{n+1}$.

In the monolithic case (e.g., the scheme (18)), the fully-discrete linear
system to be solved at each time step has the form

$\displaystyle\begin{bmatrix}\mathbf{A}_{n}&\mathbf{B}^{\top}\,\\\
\mathbf{B}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mbox{\boldmath$U$}_{n+1}\\\
\mbox{\boldmath$P$}_{n+1}\end{bmatrix}=\begin{bmatrix}\mbox{\boldmath$F$}_{n+1}\\\
\mbox{\boldmath$0$}\end{bmatrix},$

where $\mathbf{B}$ is (minus) the divergence matrix and
$\mbox{\boldmath$F$}_{n+1}$ is the right-hand side vector depending on
$\mbox{\boldmath$f$}_{n+1}$, $\mbox{\boldmath$U$}_{n}$, and $\nu_{n}$; the
velocity-velocity matrix $\mathbf{A}_{n}$ will consist of $d$ identical
blocks:

$\displaystyle\mathbf{A}_{n}=\begin{bmatrix}\mathbf{K}_{n}&\mbox{\boldmath$0$}\\\
\mbox{\boldmath$0$}&\mathbf{K}_{n}\end{bmatrix}\ \text{for}\ d=2\,,\
\,\text{or}\
\mathbf{A}_{n}=\begin{bmatrix}\mathbf{K}_{n}&\mbox{\boldmath$0$}&\mbox{\boldmath$0$}\\\
\mbox{\boldmath$0$}&\mathbf{K}_{n}&\mbox{\boldmath$0$}\\\
\mbox{\boldmath$0$}&\mbox{\boldmath$0$}&\mathbf{K}_{n}\end{bmatrix}\
\text{for}\ d=3\,,$ (42)

where
$\mathbf{K}_{n}:=\mathbf{M}+\mathbf{D}(\nu_{n})+\mathbf{C}(\mbox{\boldmath$U$}_{n})$,
with $\mathbf{M}$, $\mathbf{D}$ and $\mathbf{C}$ denoting, respectively, the
standard mass, diffusion, and convection matrices from a scalar transport
equation. The same structure will be observed for higher-order temporal
discretisations, with $\nu_{n}$ and $\mbox{\boldmath$U$}_{n}$ replaced by
their corresponding extrapolations. This block structure reduces the costs of
assembling and solving the algebraic system [20], being especially
advantageous for fractional-step schemes, as it allows us to split the
velocity update into $d$ scalar problems:

$\displaystyle\mathbf{K}_{n}\,\mbox{\boldmath$U$}^{i}_{n+1}=\mbox{\boldmath$F$}^{i}(\mbox{\boldmath$P$}_{n},\mbox{\boldmath$P$}_{n-1},\mbox{\boldmath$U$}_{n},\nu_{n})\,,\
\ i=1,...\,,d\,,$

where $\mbox{\boldmath$U$}^{i}_{n+1}$ contains the degrees of freedom of the
$i$-th velocity component. After computing all the components, the pressure is
updated by solving another scalar problem.

###### Remark 6.1

One may also wish to treat convection explicitly, to eliminate from
$\mathbf{K}_{n}$ the non-symmetric matrix $\mathbf{C}$, so that more efficient
solvers can be used. The price to pay for this choice is, in most cases, a CFL
condition, as it was the case in recent works [17, 19].

###  6.2 Generalised Newtonian fluids

So far in our analysis we have considered that the viscosity $\nu$ is fixed
and known. Nevertheless, in several interesting applications (such as
hemodynamic flows and nanofluids, see e.g. [29, 30, 31, 32]) the fluid is best
described using a generalised Newtonian constitutive law. In such a case, the
viscosity is given by

$\displaystyle\nu(\mbox{\boldmath$x$},t)=\eta(|\nabla^{\mathrm{s}}\mbox{\boldmath$u$}(\mbox{\boldmath$x$},t)|)\,,$

where $|\cdot|$ denotes the Euclidean norm, and $\eta$ is a continuous
function. A popular example is the Carreau model [33]

$\displaystyle\eta(|\nabla^{\mathrm{s}}\mbox{\boldmath$u$}|)=\nu_{\infty}+(\nu_{0}-\nu_{\infty})(1+2\lambda^{2}|\nabla^{\mathrm{s}}\mbox{\boldmath$u$}|^{2})^{-m}\,,$
(43)

where $\nu_{0}$, $\nu_{\infty}$ and $\lambda$ are positive, experimental
rheological parameters. For the common case of shear-thinning fluids
($\nu_{\infty}>\nu_{0}$), we get $\nu_{0}\leq\nu\leq\nu_{\infty}$.

The viscosity just described, although covered by our analysis, may cause
issues in the fully discrete case, where discontinuities in the velocity
gradient lead to discontinuous viscosity. Several approaches can be used to
overcome this. For example, local averaging can be used to produce a
continuous viscosity field. Alternatively, computing the
$L^{2}(\Omega)$-orthogonal projection onto a space of continuous finite
element functions would also fit our hypotheses. This last alternative has
been used, with satisfactory numerical results, for schemes related to the
ones presented herein [2, 3], and it is also the approach we follow in our
numerical experiments. The mathematical analysis of the resulting schemes has
not been carried out so far, and it is out of the scope of the present paper.

##  7 Numerical examples

In this section we present three series of numerical experiments testing the
difference in stability properties for the IMEX schemes. In accordance to the
discussion in the last section, we have considered conforming Taylor-Hood
elements in quadrilateral meshes (biquadratic velocities and continuous
bilinear pressures). For problems with a velocity-dependent viscosity we have
projected the viscosity onto the continuous bilinear finite element space.
Throughout this section, whenever an IMEX-SD formulation is mentioned, we will
be referring to one for which we have proved unconditional stability – unless
where otherwise stated.

###  7.1 Temporal convergence

Our first test aims at assessing the temporal accuracy/convergence of the IMEX
schemes presented herein. For this purpose, we use a problem with known exact
solution. The viscosity field is given by

$\nu(x,y,t)=xyf(t)+g(t)\,,$ (44)

where $f$ and $g$ are any two functions that guarantee $\nu>0$ in $\Omega$ for
all $t\in[0,T]$. With this choice, the exact solution of the Navier–Stokes
equation is given by

$\displaystyle\mbox{\boldmath$u$}=2f(t)\begin{pmatrix}y\\\
x\end{pmatrix}\quad\text{and}\quad p=(C-2xy)f^{\prime}(t)\,,$

with $\mbox{\boldmath$f$}\equiv\mbox{\boldmath$0$}$ and the boundary
conditions chosen accordingly; the constant $C$ is set so that $p\in
L^{2}_{0}(\Omega)$ for all $t$. We have chosen this solution since, in space,
it belongs to the finite element space, so the entire error stems from time
discretisation. We set $\Omega=(0,1)^{2}$, $T=10$, $f(t)=\sin^{2}t$ and
$g(t)\equiv 0.001$. The square domain is discretised uniformly with $4\times
4$ square elements, which completely resolves (at each time instant) the exact
solution in space. The time-step size is halved at each refinement level,
starting from $\tau=1$. In Figure 1 we depict the errors in velocity and
pressure at $t=10$. For all schemes considered, the convergence is at least
linear in both fields, with the monolithic GL being the most accurate for
velocity (although the remaining schemes show a higher order of convergence),
and a comparable performance of all methods for pressure.

Figure 1: Temporal convergence study for the IMEX schemes. The quadratic
velocity convergence of the fractional-step methods is probably an initial
superconvergence that would break down with further temporal refinement.

###### Remark 7.1

It is interesting to notice that all methods present a very stable behaviour
even for the larger time step sizes. This, despite the fact that the viscosity
field considered in this example is such that
$\nu_{\mathrm{min}}/\|\nabla\nu\|^{2}_{\infty}=5\times 10^{-4}$. According to
our analysis, this would require $\tau<0.001$ to guarantee the stability of
the fractional-step GL method. That, however, was not observed in this study.

###  7.2 Stability test

To put the stability of the fractional-step IMEX schemes to a more challenging
test, we tackle a problem with less smooth solution. The lid-driven cavity
flow is one of the most popular benchmarks for generalised Newtonian flow
solvers, and a common setup involves a power-law fluid:

$\displaystyle\eta(|\nabla^{\mathrm{s}}\mbox{\boldmath$u$}|)=\kappa\big{|}2^{\frac{1}{2}}\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\big{|}^{n-1}\,,$
(45)

where $\kappa$, $n$ and $\nu_{\mathrm{min}}$ are positive constants. For the
square cavity $\Omega=(0,1)^{2}$, the boundary conditions are a horizontal
velocity $U(t)$ on the lid ($y=1$), and no slip on the bottom and lateral
walls. We consider a reference shear-thickening fluid with $\kappa=0.01$ and
$n=1.5$ [34]. To avoid a zero viscosity, we replace (45) by

$\displaystyle\eta(|\nabla^{\mathrm{s}}\mbox{\boldmath$u$}|)=\mathrm{max}\left\\{10^{-5},\,\kappa\big{|}2^{\frac{1}{2}}\nabla^{\mathrm{s}}\mbox{\boldmath$u$}\big{|}^{n-1}\right\\}.$

The lid velocity is prescribed according to

$\displaystyle
U(t)=\frac{1}{2}\left\\{1+\text{sign}(t-t^{\star})+\left[1-\text{sign}(|2t-t^{\star}|-t^{\star})\right]\sin^{2}\Big{(}\frac{\pi
t}{2t^{\star}}\Big{)}\right\\},$ (46)

which ramps up smoothly from $U=0$ (at $t=0$) to $U=1.0$ (for $t\geq
t^{\star}=2$). We run the simulation up to $T=30$, aiming for a steady-state
flow that can be compared to reference stationary results [34, 35]. Two
uniform meshes are considered: one with $200\times 200$ elements, and one
twice as fine. To test the stability of the different schemes, we take the
time-step size rather large, namely $\tau=0.25$. The aim of this test is
twofold: firstly, we want to verify whether the time-step restriction (41)
found for the fractional-step GL scheme is sharp – and we do that by tackling
a problem where $\|\nabla\nu\|_{\infty}/\nu_{\mathrm{min}}\rightarrow\infty$.
Second, we aim to illustrate the poor stability of IMEX-SD schemes when the
coupling term $\nabla\cdot(\nu\nabla^{\top}\mbox{\boldmath$u$})$ is made fully
explicit (and, since the viscosity $\nu_{n+1}$ is an unknown of the problem,
we can only consider the scheme (21) as IMEX-SD). To do that, three schemes
are considered:

  1. 1.

Fractional-step IMEX, generalised Laplacian form (37).

  2. 2.

Fractional-step IMEX, stress-divergence form using only $\nu_{n}$ (not
$\nu_{n+1}$).

  3. 3.

Monolithic, fully coupled stress-divergence formulation (19).

The third variant is meant as a control case, since it is the method for which
we have the strongest stability estimate.

To verify the accuracy of our solution, we use the location of the primary
vortex as benchmark. Table 2 shows the comparison between our results (with
the coarser mesh) and reference ones [34, 35]. Among the three methods we
used, the two with proven stability are in good agreement with the references;
on the other hand, for the IMEX-SD (21) (with fully explicit viscosity) no
primary vortex position could be clearly identified. Figure 2 shows how the
kinetic energy $E(t):=\|\mbox{\boldmath$u$}\|^{2}/2$ evolves with time: while
the two stable schemes reach a steady-state level, IMEX-SD (21) exhibits a
clearly unstable energy growth (notice that the spatial refinement does not
remedy the instability). The two negative results for the IMEX-SD (21) are
important as they hint strongly that, in general, one cannot guarantee
unconditional stability of schemes such as (21). Furthermore, Figure 3 depicts
the viscosity field at the final time, where we see the spikes near the
corners (which lead to very high gradients). The time-step size used for all
methods in this experiment ($\tau=0.25$) clearly violates (41). Yet, the
fractional-step GL method again showed a stable behaviour, which hints at the
possibility that (41) is pessimistic.

Table 2: Lid-driven cavity flow with power-law fluid: steady-state position $(x,y)$ of the primary vortex for three schemes: fully coupled SD (19), fractional-step IMEX-GL (37) and fractional-step IMEX-SD using only $\nu_{n}$. For the latter, the solution did not converge to a steady-state, so it was not possible to determine clear vortex locations. Ref. [34] | Ref. [35] | Fully coupled | IMEX-GL | IMEX-SD, $\nu_{n}$
---|---|---|---|---
$(0.565,0.724)$ | $(0.5628,0.7282)$ | $(0.562,0.728)$ | $(0.560,0.728)$ | undetermined
Figure 2: Lid-driven cavity flow with power-law fluid: temporal evolution of
the kinetic energy, confirming that a stress-divergence scheme with fully
explicit viscosity can be unstable. Figure 3: Lid-driven cavity flow with
power-law fluid: steady-state viscosity field obtained using the IMEX-GL
method.

###  7.3 Idealised hemodynamic flow

Our last numerical test tackles a truly dynamic flow problem. A popular
application for variable viscosity is hemodynamics, with blood usually
modelled as a shear-thinning fluid. We have applied the Carreau model (43)
with the following parameters: $\nu_{\infty}=3.286\times 10^{-6}$ m2/s,
$\nu_{0}=53.33\times 10^{-6}$ m2/s , $\lambda=3.313$ s and $m=0.3216$ [36].
Combined with those, we consider an idealised two-dimensional aneurysm, whose
geometry is illustrated in Figure 4. The upper wall is described by the smooth
curve

$\displaystyle
y(x)=H+\frac{H}{4}[1-\mathrm{sign}(|2x-6H|-3H)]\cos^{2}\Big{(}\pi\frac{x-3H}{3H}\Big{)}\,,$

with $H=2.5\times 10^{-3}$ m.

Figure 4: Ccomputational domain representing an idealised ICA aneurysm.

To simulate a physiological regime [37], we prescribe a (ramped up) pulsating
horizontal inflow profile with period $T_{p}=0.917$ s:

$\displaystyle
u(y,t)|_{x=0}=\frac{1}{2}\Bigg{[}1-\Big{(}\frac{2y}{H}-1\Big{)}^{2}\Bigg{]}U(t)f(t)\,,$

with $U(t)$ as in (46), $t^{\star}=T_{p}$. The function $f(t)$ is a
trigonometric polynomial fitted from internal carotid artery (ICA)
measurements [38]:

$\displaystyle
f(t)=10^{-4}\Bigg{[}\frac{a_{0}}{2}-\sum_{k=1}^{12}a_{k}\cos\left(\frac{2k\pi
t}{T_{p}}\right)+b_{k}\sin\left(\frac{2k\pi t}{T_{p}}\right)\Bigg{]}\,,$

where $a_{0}=13349$, and the remaining coefficients are listed on Table 3.
That leads to a peak inflow speed of $U_{\mathrm{peak}}=0.5$ m/s. At the
outlet ($x=6H$), we prescribe the outflow condition

$\displaystyle(\nu\nabla\mbox{\boldmath$u$})\mbox{\boldmath$n$}-p\mbox{\boldmath$n$}=\mbox{\boldmath$0$}\,,$
(47)

while the bottom and top walls are no-slip boundaries
($\mbox{\boldmath$u$}=\mbox{\boldmath$0$}$).

Table 3: Trigonometric coefficients for the ICA inflow waveform ($a_{0}=13349$). $k$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---
$a_{k}$ | 695.98 | 905.24 | 452.42 | 754.46 | 164.6 | 133.57 | 79.4 | 74.40 | 36.72 | 38.80 | 42.34 | 3.67
$b_{k}$ | $-2067.65$ | $-496.92$ | $-340.74$ | 202.96 | 98.82 | 249.55 | 81.7 | 67.84 | 9.95 | 15.34 | $-3.49$ | $-3.22$

For this test we consider the GL formulation, as it offers a natural way to
enforce the outflow condition (47). We use the monolithic IMEX-GL scheme (26)
with different time-step sizes $\tau$ and a mesh with 38,400 elements. Figure
5 shows, over one cardiac cycle, the mean pressure drop across the aneurysm:

$\displaystyle\Delta
p=\frac{1}{H}\left(\int_{\Gamma_{\mathrm{in}}}p\,\mathrm{d}\Gamma-\int_{\Gamma_{\mathrm{out}}}p\,\mathrm{d}\Gamma\right),$

which is an important clinical marker [39]. The results confirm very good
temporal stability for the numerical solutions, even for time-step sizes as
large as $\tau=0.05$. For $t=6T$, Figure 6 shows the horizontal velocity
profiles at different sections of the vessel, showcasing also good spatial
stability.

Figure 5: Carreau fluid through Idealised ICA aneurysm: normalised pressure
difference (inlet minus outlet) over one cardiac cycle, for different time-
step sizes $\tau$. All three computations are based on the monolithic IMEX
scheme in GL form, which allows enforcing the outflow condition (47)
naturally. Figure 6: Carreau fluid through Idealised ICA aneurysm: horizontal
velocity profiles at $t=6T$, for a time-step size $\tau=0.002$.

##  8 Concluding remarks

In this work, we have presented, analysed and tested different implicit-
explicit schemes for the variable-viscosity Navier–Stokes equations. While
IMEX methods for constant viscosity focus on making the convective velocity
(or different stabilisation terms associated to convection) explicit, we split
also the viscous term into implicit and explicit parts. This renders the
velocity subsystem block-diagonal, thereby enabling simpler, more efficient
linear solvers. The temporal stability of the resulting schemes varies
dramatically according to the term made explicit, and in which form. To
investigate that, we have carried out a semi-discrete analysis of prototypical
first-order schemes. Our results show that, when considering a stress-
divergence (SD) description of the viscous term, naive extrapolations can
cause instabilities. While unconditionally stable IMEX-SD schemes are
possible, we have shown that they require a pre-computed viscosity field,
which is not realistic in many applications. On the other hand, a generalised
Laplacian (GL) formulation of the viscous term leads to stable IMEX methods
requiring only known viscosity values. For a monolithic and a fractional-step
IMEX-GL scheme, respectively, we have proved unconditional and conditional
stability. Our numerical tests, however, suggest that the time-step
restriction derived for the fractional-step variant may be an analytical
artefact.

Our theoretical and numerical results contribute towards breaking the paradigm
that considers inevitable the use of fully coupled, elasticity-like
formulations for variable-viscosity flows. The GL variants presented in this
work have provided accurate and stable results, performing well even for
problems with non-smooth solutions. Several questions remain open though.
Proposing and analysing second- and higher-order schemes is a challenging and
interesting open question. In addition, analysing the stability and
convergence of fully discrete versions of the schemes proposed in this work is
also of interest, especially for the case of nonlinear viscosity models.

## Acknowledgments

The work of GRB has been partially supported by the Leverhulme Trust Research
Project Grant No. RPG-2021-238. EC acknowledges the support given by the
Agencia Nacional de Investigación y Desarrollo (ANID) through the project
FONDECYT 1210156.

## References

  * Carey et al. [1989] G.F. Carey, K.C. Wang, and W.D. Joubert.  Performance of iterative methods for newtonian and generalized newtonian flows.  _Internat. J. Numer. Methods Fluids_ , 9(2):127–150, 1989.
  * Schussnig et al. [2021] R. Schussnig, D.R.Q. Pacheco, and T.-P. Fries.  Robust stabilised finite element solvers for generalised newtonian fluid flows.  _J. Comput. Phys._ , 442:110436, 10 2021.
  * Pacheco et al. [2021] D.R.Q. Pacheco, R. Schussnig, and T.-P. Fries.  An efficient split-step framework for non-Newtonian incompressible flow problems with consistent pressure boundary conditions.  _Comput. Methods Appl. Mech. Engrg._ , 382:113888, 2021\.
  * Deteix and Yakoubi [2018] J. Deteix and D. Yakoubi.  Improving the pressure accuracy in a projection scheme for incompressible fluids with variable viscosity.  _Appl. Math. Lett._ , 79:111–117, 2018.
  * Plasman et al. [2020] L. Plasman, J. Deteix, and D. Yakoubi.  A projection scheme for Navier-Stokes with variable viscosity and natural boundary condition.  _International Journal for Numerical Methods in Fluid_ , 92(12):1845–1865, 2020.
  * Anaya et al. [2021] V. Anaya, R. Caraballo, B. Gómez-Vargas, D. Mora, and R. Ruiz-Baier.  Velocity-vorticity-pressure formulation for the oseen problem with variable viscosity.  _Calcolo_ , 58:44, 12 2021.
  * Anaya et al. [2023] V. Anaya, R. Caraballo, R. Ruiz-Baier, and H. Torres.  On augmented finite element formulations for the navier–stokes equations with vorticity and variable viscosity.  _Comput. Math. Appl._ , 143:397–416, 8 2023.
  * Chorin [1967] A.J. Chorin.  The numerical solution of the navier-stokes equations for an incompressible fluid.  _Bull. Amer. Math. Soc._ , 73:928–931, 1967.
  * Glowinski [2003] Roland Glowinski.  Finite element methods for incompressible viscous flow.  In _Numerical Methods for Fluids (Part 3)_ , volume 9 of _Handbook of Numerical Analysis_ , pages 3–1176. Elsevier, 2003.
  * Chorin [1968] A.J. Chorin.  Numerical solution of the Navier-Stokes equations.  _Math. Comp._ , 22:745–762, 1968.
  * Temam [1969] R. Temam.  Sur l’approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II.  _Arch. Ration. Mech. Anal._ , 33:377–385, 1969.
  * Guermond et al. [2006] J.L. Guermond, P. Minev, and Jie Shen.  An overview of projection methods for incompressible flows.  _Comput. Methods Appl. Mech. Engrg._ , 195(44-47):6011–6045, 2006.
  * Badia and Codina [2007] S. Badia and R. Codina.  Algebraic pressure segregation methods for the incompressible Navier-Stokes equations.  _Arch. Comput. Methods Eng._ , 15:1–52, 2007.
  * Díaz et al. [2023] F. Díaz, E. Castillo, R.C. Cabrales, and N.O. Moraga.  Third-order in time non-relaxed finite volume fractional step schemes for unsteady incompressible flows.  _Comput. Math. Appl._ , 2023.  Accepted.
  * Boscarino et al. [2013] S. Boscarino, L. Pareschi, and G. Russo.  Implicit-explicit runge–kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit.  _SIAM J. Sci. Comput._ , 35(1):A22–A51, 2013\.
  * Hochbruck and Leibold [2021] M. Hochbruck and J. Leibold.  An implicit-explicit time discretization scheme for second-order semilinear wave equations with application to dynamic boundary conditions.  _Numer. Math._ , 147(4):869–899, 2021.
  * Burman and Guzmán [2022] E. Burman and J. Guzmán.  Implicit-explicit multistep formulations for finite element discretisations using continuous interior penalty.  _ESAIM Math. Model. Numer. Anal._ , 56:349–383, 1 2022\.
  * Guesmi et al. [2023] M. Guesmi, M. Grotteschi, and J. Stiller.  Assessment of high‐order IMEX methods for incompressible flow.  _Internat. J. Numer. Methods Fluids_ , 2 2023.
  * Burman et al. [2023] E. Burman, D. Garg, and J. Guzman.  Implicit-explicit time discretization for oseen’s equation at high reynolds number with application to fractional step methods.  _SIAM J. Numer. Anal._ , In press, 2023.
  * John [2016] V. John.  _Finite element methods for incompressible flow problems_.  Springer, Basel, 2016.
  * Stiller [2020] J. Stiller.  A spectral deferred correction method for incompressible flow with variable viscosity.  _J. Comput. Phys._ , 423:109840, 12 2020.
  * Olshanskii et al. [2009] M. Olshanskii, G. Lube, T. Heister, and J. Löwe.  Grad–div stabilization and subgrid pressure models for the incompressible navier–stokes equations.  _Comput. Methods Appl. Mech. Engrg._ , 198:3975–3988, 11 2009.
  * Wang, Haijin et al. [2016] Wang, Haijin, Wang, Shiping, Zhang, Qiang, and Shu, Chi-Wang.  Local discontinuous galerkin methods with implicit-explicit time-marching for multi-dimensional convection-diffusion problems.  _ESAIM: M2AN_ , 50(4):1083–1105, 2016.
  * Wang and Shu [2020] Zhang Q. Wang S. Wang, H. and C.-W. Shu.  Local discontinuous galerkin methods with explicit-implicit-null time discretizations for solving nonlinear diffusion problems.  _Sci. China Math._ , 63:183–204, 2020.
  * Tan et al. [2022] Meiqi Tan, Juan Cheng, and Chi-Wang Shu.  Stability of high order finite difference and local discontinuous galerkin schemes with explicit-implicit-null time-marching for high order dissipative and dispersive equations.  _Journal of Computational Physics_ , 464:111314, 2022.
  * Ern and Guermond [2021] Alexandre Ern and Jean-Luc Guermond.  _Finite Elements I_.  Springer, 2021.
  * Heywood and Rannacher [1990] J.G. Heywood and R. Rannacher.  Finite-element approximation of the nonstationary navier–stokes problem. part iv: Error analysis for second-order time discretization.  _SIAM J. Numer. Anal._ , 27:353–384, 4 1990.
  * Deteix and Yakoubi [2019] J. Deteix and D. Yakoubi.  Shear rate projection schemes for non-Newtonian fluids.  _Comput. Methods Appl. Mech. Engrg._ , 354:620–636, 2019\.
  * Abugattas et al. [2020] C. Abugattas, A. Aguirre, E. Castillo, and M. Cruchaga.  Numerical study of bifurcation blood flows using three different non-newtonian constitutive models.  _Appl. Math. Model._ , 88:529–549, 2020.
  * Farías et al. [2023] Catalina Farías, Camilo Bayona-Roa, Ernesto Castillo, Roberto C. Cabrales, and Ricardo Reyes.  Reduced order modeling of parametrized pulsatile blood flows: Hematocrit percentage and heart rate.  _International Journal of Engineering Science_ , 193:103943, 2023.
  * González et al. [2021] A. González, O. Ruz, and E. Castillo.  Numerical study of the fluid dynamics and heat transfer for shear-thinning nanofluids in a micro pin-fin heat sink.  _Case Studies in Thermal Engineering_ , 28:101635, 2021\.
  * Aguirre et al. [2022] I. Aguirre, A. González, and E. Castillo.  Numerical study on the use of shear-thinning nanofluids in a micro pin-fin heat sink including vortex generators and changes in pin shapes.  _Journal of the Taiwan Institute of Chemical Engineers_ , 136:104400, 2022.
  * Galdi et al. [2008] G.P. Galdi, R. Rannacher, A.M. Robertson, and S. Turek.  _Hemodynamical Flows_ , volume 37 of _Oberwolfach Seminars_.  Birkhäuser, Basel, 2008.
  * Neofytou [2005] P. Neofytou.  A 3rd order upwind finite volume method for generalised newtonian fluid flows.  _Adv. Eng. Softw._ , 36:664–680, 10 2005.
  * Li et al. [2014] Q. Li, N. Hong, B. Shi, and Z. Chai.  Simulation of power-Law fluid flows in two-dimensional square cavity using multi-relaxation-time lattice Boltzmann method.  _Commun. Comput. Phys._ , 15:265–284, 2014.
  * Cho and Kensey [1991] Y.I. Cho and K.R. Kensey.  Effects of the non-newtonian viscosity of blood on flows in a diseased arterial vessel. part 1: Steady flows.  _Biorheology_ , 28:241–262, 6 1991.
  * John et al. [2017] L. John, P. Pustějovská, and O. Steinbach.  On the influence of the wall shear stress vector form on hemodynamic indicators.  _Comput. Vis. Sci._ , 18(4-5):113–122, 2017.
  * Thiriet [2008] M. Thiriet.  _Biology and Mechanics of Blood Flows_.  Springer, New York, NY, 2008.
  * Pacheco [2022] D.R.Q. Pacheco.  On the numerical treatment of viscous and convective effects in relative pressure reconstruction methods.  _Internat. J. Numer. Methods Biom. Engrg._ , 38, 3 2022.