# Discretisation of an Oldroyd-B viscoelastic fluid flow using a Lie
derivative formulation

Ben S. Ashby  Institute for Mathematical Innovation,
University of Bath, Bath, UK.  Tristan Pryer  Institute for Mathematical
Innovation,
University of Bath, Bath, UK.

###### Abstract

In this article we present a numerical method for the Stokes flow of an
Oldroyd-B fluid. The viscoelastic stress evolves according to a constitutive
law formulated in terms of the upper convected time derivative. A finite
difference method is used to discretise along fluid trajectories to
approximate the advection and deformation terms of the upper convected
derivative in a simple, cheap and cohesive manner, as well as ensuring that
the discrete conformation tensor is positive definite. A full implementation
with coupling to the fluid flow is presented, along with detailed discussion
of the issues that arise with such schemes. We demonstrate the performance of
this method with detailed numerical experiments in a lid-driven cavity setup.
Numerical results are benchmarked against published data, and the method is
shown to perform well in this challenging case.

##  1 Introduction

Viscoelastic fluids produce a number of phenomena distinct from the familiar
behaviour of Newtonian fluids. A broad class of such fluids is that of polymer
solutions, for which these phenomena are the result of long chain molecules in
the fluid that are large enough to have a meaningful effect on the fluid flow.
Such fluids are found in industry in the form of molten plastics or machine
lubricants, as well as several biological fluids such as blood. Polymeric
fluids are modelled by augmenting the usual fluid flow equations for
conservation of mass and momentum with a constitutive law for the polymeric
stress, which results in a system of either partial differential or integro-
differential equations. This article will focus on the former, with the
constitutive law given by a hyperbolic partial differential equation (PDE).

In Newtonian fluid dynamics, the _Reynolds number_ indicates the relative
scaling of inertial terms to viscous terms in the Navier-Stokes equations, and
can be used to predict flow behaviour and transition between different regimes
such as laminar and turbulent flow. In viscoelastic flows, there are two
dimensionless groups which are used to predict flow behaviour, the
_Weissenberg number_ , $\operatorname{Wi}$ and the _Deborah number_ ,
$\operatorname{De}$. The Weissenberg number quantifies the effect of
nonlinearities that arise due to non-Newtonian normal stress differences [11,
33]. The Deborah number has been interpreted as the ratio of the magnitudes of
elastic and viscous forces [4]. In the limit as $\operatorname{De}\to 0$,
Newtonian flow behaviour is recovered. Experimental results for dilute polymer
solutions reported in [32] suggest that a transition to turbulent flow
behaviour may occur when elastic effects become dominant over viscous, an
effect known as elastic instability, even for small Reynolds number. In
complex flow geometries and time-dependent flows, it may be difficult to
define these dimensionless groups appropriately, and indeed one may be more
relevant than the other. However, in keeping with much of the literature on
numerical methods for non-Newtonian flows, we will mainly consider the
Weissenberg number $\operatorname{Wi}$, defined as

$\operatorname{Wi}:=\frac{\lambda U}{L},$

where $\lambda$ is the relaxation time and $U$ and $L$ are characteristic
velocity and length scales of the problem.

Numerical methods posed for the solution of viscoelastic fluid problems,
especially using the Oldroyd-B model, have been observed to suffer from the
so-called high Weissenberg number problem, where numerical simulations fail to
converge at moderate values of the Weissenberg number. There are several
potential causes of the high Weissenberg number problem such as the fact that
many viscoelastic fluid flow problems exhibit strong, spatially exponential
layers which require fine mesh resolution to adequately resolve. In [34], the
authors observe a relationship between the Weissenberg number and the
principle eigenvalue, determined by numerical simulation and curve fitting.
They find that this maximum value increases exponentially with respect to
$\operatorname{Wi}$. It is therefore a significant computational cost to
resolve this layer for larger Weissenberg number. It has also been suggested
that the Oldroyd-B model itself is insufficient to properly describe
viscoelastic behaviour. Indeed, extensional viscosity of the Oldroyd-B model
blows up at a finite strain rate [21]. Despite this, the Oldroyd-B model
remains a popular test case for numerical modelling. This point is discussed
in greater detail in §2.

It is also well known that for the Oldroyd-B model, if the initial data for
the conformation tensor is positive definite, then the conformation tensor
remains so for all time [13, 40]. Loss of positive definiteness at the
discrete level is often a precursor to numerical blowup, and indeed it was
shown in [42] that the Upper-convected Maxwell model (which is the model that
results if the solvent contribution is neglected in the Oldroyd-B model)
becomes ill-posed. It therefore appears that designing schemes which preserve
the positive definiteness of the discrete conformation tensor is important not
only to maintain the physical structures of the model, but also crucial for
numerical stability.

These phenomena motivated what seems to be the dominant method for solving
viscoelastic constitutive laws in differential form: the log-conformation
representation, introduced in [15]. In this approach, a differential equation
is derived for the tensor logarithm of the (positive definite) conformation
tensor, thus guaranteeing that the numerical approximation to the conformation
tensor will remain positive definite. It is argued that this method is more
appropriate for resolving exponential spatial layers, and indeed has provided
numerical simulations with relatively large Weissenberg numbers [16, 23] which
are stable, though questions of accuracy remain. It is however suggested by
the results of [45] that careful application of the method is required to
preserve physical characteristics of the model. In particular, they report
reduced turbulent drag reduction. In addition, solving for the logarithm of
the conformation tensor introduces additional nonlinearities into the
numerical model that are not present in the physical model, increasing
computational expense. Other schemes which preserve positivity have been
derived using for example the square root of the conformation tensor [28] and
the use of the hyperbolic tangent transformation [24].

Schemes based on a Lie derivative formulation offer a potential solution to
the issues described above. While characteristic schemes are motivated by the
idea of discretising the material derivative as a single object after
expressing it as a rate along fluid trajectories, the upper-convected
derivative is similarly expressible as a rate along fluid paths, this time
including information on deformation of the moving fluid. Indeed, when
expressed in a convected coordinate system, the Oldroyd-B constitutive law is
a tensor-valued ordinary differential equation. The method presented in this
work can be thought of as resulting from an approximation of a rate of change
of stress that also incorporates rate of change of the convected coordinate
system used to express it. Discretising the Oldroyd-B constitutive law with an
appropriately chosen finite difference methods provides a natural way to
preserve positive definiteness of the conformation tensor. To our knowledge,
the idea of conducting computations for viscoelastic fluids by discretising
the upper convected derivative expressed as a Lie derivative was first
proposed in [35]. Theory was developed in [25, 26] for a finite element method
in the framework of Ricatti equations. A finite difference scheme for the
Oldroyd-B constitutive law with prescribed velocity field was presented with
truncation analysis in [30]. These methods can be viewed as an extension of
characteristic schemes, which are well established in fluid dynamics due to
their favourable stability properties, see [37] for an early reference, with
analysis of a finite element method given in [44]. For example, when applied
to the transport equation, stable solutions are obtained even for large CFL
numbers. These methods have also been used to discretise the advective part of
the upper convected derivative in computations involving viscoelastic fluids
[8].

The discretisation of this problem in the Lie derivative framework has
distinct advantages and unique challenges. One of the goals of the paper is to
compare these features with other common approaches, and to test the limits of
Lie derivative discretisation methods. We present a scheme based upon the Lie
derivative where the conformation tensor is evolved pointwise, and spatial
discretisation enters only through the approximation of the deformation
gradient tensor and particle trajectories and interpolation. We note that
similar schemes have appeared in the literature [3, 25, 26, 30], but our aim
in this work is to present results from a full implementation in a challenging
test case that can be compared with other published data.

In this article, we provide an accessible overview of viscoelastic fluid
dynamics, particularly focused on polymer solutions modeled by the Oldroyd-B
constitutive law. We develop a numerically tractable scheme for the model and
prove the well-posedness. Our developed scheme also performs favorably in
comparison to more computationally intensive alternatives that typically
require high performance computing to realise. A thorough comparative study is
presented, demonstrating the robustness of our approach. Furthermore, we
discuss appropriate mesh design critical for capturing the complex dynamics of
viscoelastic flows. This article will focus on the lid-driven cavity problem:
a challenging test problem for non-Newtonian flows due to strong shearing that
occurs near the moving lid, as well as strong stress gradients at the corners.
It is a well-studied benchmark test case in the numerical literature on
simulation of the Oldroyd-B model [7, 10, 16, 19, 34, 43], and serves as a
useful example to illustrate key differences between Newtonian and non-
Newtonian flows.

The outline of the rest of the article is as follows. In §2, we describe the
Oldroyd-B model. Key concepts from continuum mechanics are introduced in §3,
where a reformulation of the upper convected time derivative is described. In
§4, a semidiscrete upper convected time derivative is introduced, and its
consistency and approximation order is given. We discuss fully discrete
numerical methods in §5. We work with a decoupled formulation, and so we
discuss solution of the Stokes and Oldroyd-B systems separately. Descriptions
of computational aspects and implementation details, as well as numerical
experiments are presented in §6.

##  2 The Oldroyd-B model

Let $\Omega\subset\mathbb{R}^{d}$ be a convex domain with polygonal boundary,
and let $\bm{u}:\Omega\times[0,T]\to\mathbb{R}^{d}$ be the spatial velocity
field, with $p:\Omega\times[0,T]\to\mathbb{R}$ the pressure. The (symmetric)
_conformation tensor_ $\bm{\sigma}:\Omega\to\mathbb{R}^{d^{2}}$ describes the
non-Newtonian component of the stress. We let $\bm{\varepsilon}(\bm{u})$ be
the rate of strain tensor, or symmetrised velocity gradient:

$\bm{\varepsilon}(\bm{u})=\frac{1}{2}\left(\nabla\bm{u}+(\nabla\bm{u})^{T}\right).$

The Oldroyd-B model for the flow of a fluid in $\Omega$ consists of the
Navier-Stokes equations for incompressible flow given by

$\begin{split}\operatorname{Re}\left(\frac{\partial\bm{u}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{u}\right)&=-\nabla
p+2\beta\operatorname{div}\\!\left({\bm{\varepsilon}(\bm{u})}\right)+\frac{1-\beta}{\operatorname{Wi}}\operatorname{div}\\!\left({\bm{\sigma}}\right)\\\
\operatorname{div}\bm{u}&=0.\end{split}$ (1)

This is coupled to a hyperbolic constitutive law for the evolution of the
stress, and is given below in non-dimensional form.

$\frac{\partial\bm{\sigma}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{\sigma}-(\nabla\bm{u})\bm{\sigma}-\bm{\sigma}(\nabla\bm{u})^{T}=-\frac{1}{\operatorname{Wi}}(\bm{\sigma}-\bm{I}).$
(2)

The system has three parameters, namely the familiar Reynolds number
$\operatorname{Re}$, the elasticity ratio $\beta$ and the _Weissenberg number_
, $\operatorname{Wi}$. The regime of interest in this work is low
$\operatorname{Re}$ and high $\operatorname{Wi}$. We will assume that
$\beta>0$ to avoid the complications that can arise from the viscous term
being removed and the system becoming fully hyperbolic (these difficulties are
discussed in more detail in [6]). We note that the conformation tensor is
related to the polymeric extra-stress tensor $\bm{\tau}$ by

$\bm{\sigma}=\bm{I}+\frac{\operatorname{Wi}}{1-\beta}\bm{\tau}.$

We note that in the context in this work (i.e. a lid-driven cavity), the
Weissenberg and Deborah numbers are related to each other via the aspect ratio
of the cavity [33], and are in fact equal for a cavity of square cross-
section.

We prescribe Dirichlet boundary conditions $\bm{w}$ on $\partial\Omega$ for
the velocity such that $\bm{w}\cdot\bm{n}=0$ on $\partial\Omega$, where
$\bm{n}$ is the outward facing normal vector to the boundary, i.e. the problem
is driven by tangential boundary values for the velocity. We note that since
Dirichlet boundary conditions are prescribed on the whole boundary, to ensure
that there exists a uniquely determined pressure field we include the
constraint that the pressure has mean zero. Initial conditions for the fluid
velocity and conformation tensor are given by
$\bm{u}(\bm{x},0)=\bm{u}^{0}(\bm{x})$ and
$\bm{\sigma}(\bm{x},0)=\bm{\sigma}^{0}(\bm{x})$ respectively.

The left hand side of Equation (2) is the _upper convected time derivative_ ,
introduced by Oldroyd [31] and discussed in greater detail in §3. See also
[21] for an overview. It gives the rate of change of a tensor quantity in a
convected coordinate system which moves and deforms with the fluid. Such a
rate is appropriate to derive constitutive laws which exhibit material frame
indifference that is, the principle that the physics of the constitutive law
should be independent of the frame of the observer.

######  Remark 2.1 (Boundary conditions for the conformation tensor).

The problem described above is an enclosed flow. Since there is no inflow
boundary, and since the constitutive relation (2) is hyperbolic, no boundary
condition is required for the conformation tensor. The case of inflow boundary
conditions for the velocity, which we do not consider here, requires
additional boundary conditions for the conformation tensor. This is not a well
studied problem and another potential source of instabilities in these flows.

To complete the exposition, it is appropriate to discuss some of the
assumptions made in the derivation of the Oldroyd-B model. The introduction of
the upper convected time derivative by Oldroyd [31] assumes that the physics
of the constitutive relation are formulated in terms of the evolution of
contravariant components of a second order tensor in the convected frame
moving and deforming with the fluid. This is a choice, justified by
experiments exhibiting rod-climbing (a phenomenon that the Oldroyd-B model
replicates), rather than physics, and other choices may be more appropriate
(see the discussion in [21]).

The model has also been derived from kinetics using a dumbbell and spring
model in which no limit is placed upon the extension of the spring. However,
if the flow is strong (that is, if the largest eigenvalue of $\nabla\bm{u}$
exceeds $\tfrac{1}{2\lambda}$, where $\lambda$ is the relaxation time [21]),
the molecule length increases without limit. It has been shown that, as a
result of this, infinite stress can occur in the interior of a steady flow. In
many regimes however, the Oldroyd-B fluid is a good approximation of a Boger
fluid. It is the simplest differential viscoelastic model available which
enjoys material frame indifference (simpler models do exist such as the linear
Maxwell model, see [39, Chapter 2] for an overview, but do not have this
property). It is also popular, with a wide range of benchmark numerical data
available in the literature. It therefore represents a prototypical model for
the derivation and testing of numerical methods, see for example [16, 24, 26,
29, 36] and many others.

######  Remark 2.2.

Taking the kinetic viewpoint, the Oldroyd-B model can be improved by
considering for example variants of the _Finitely Extensible Nonlinear
Elastic_ (FENE) model for the spring force in the polymer molecules, which
limits molecules to be finitely extensible and avoiding the issue of unbounded
stress growth. This does however result in a more complicated nonlinearity in
the constitutive law. For example, the FENE-P model has constitutive law given
by

$\frac{\partial\bm{\sigma}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{\sigma}-(\nabla\bm{u})\bm{\sigma}-\bm{\sigma}(\nabla\bm{u})^{T}=\frac{1}{\operatorname{Wi}}\left(\bm{I}-\frac{1}{1-\frac{1}{b^{2}}\operatorname{tr}(\bm{\sigma})}\bm{\sigma}\right),$

where the additional parameter $b$ represents the maximum extensibility of the
polymer molecules. The extension of the current work to more realistic models
will be the subject of further research.

##  3 A Lie derivative framework for the upper convected time derivative

In his seminal paper [31], Oldroyd described the importance of formulating a
time derivative with which to construct constitutive laws for non-Newtonian
stresses that does not depend on the choice of reference frame. Such
derivatives have to take account of the moving fluid flow and its deformation.
The _upper convected time derivative_ , one of two derivatives introduced by
Oldroyd, is one of (infinitely) many such time derivatives. It is now well
known that these rates are special cases of Lie derivatives (see eg [22,
chapter II]), and formulating the constitutive law in a way that respects this
can have structure-preserving advantages. In this section we give the concepts
necessary to formulate the constitutive law for a non-Newtonian fluid as a Lie
derivative.

###  3.1 Flow map and deformation gradient

Let $\bm{y}$ be the flow map for the fluid motion. That is, the fluid particle
which at time $t$ has position $\bm{x}$ has position $\bm{y}(\bm{x},t;s)$ at
time $s$. The flow map satisfies the following ordinary differential equation

$\begin{split}\partial_{s}\bm{y}(\bm{x},t;s)&=\bm{u}(\bm{y}(\bm{x},t;s),s)\\\
\bm{y}(\bm{x},t;t)&=\bm{x}.\end{split}$

One can express the material derivative in terms of the flow map by examining
the infinitesimal derivative, that is

$\left.\frac{\partial}{\partial
s}\bm{u}(\bm{y}(\bm{x},t;s),s)\right|_{s=t}=\left(\partial_{t}\bm{+}\bm{u}\cdot\nabla\right)\bm{u}(x,t).$

We define the velocity gradient with the convention that it has components

$(\nabla\bm{u})_{ij}=\frac{\partial u_{i}}{\partial x_{j}},$

and make note that some authors define the velocity gradient as the transpose
of the above tensor. The deformation gradient tensor describes the deformation
experienced by a fluid parcel in a neighbourhood of $(\bm{x},t)$ between times
$s_{1}$ and $s_{2}$. In a Cartesian coordinate system, it has components

$\bm{F}_{ij}(\bm{x},t;s_{1},s_{2})=\frac{\partial\bm{y}_{i}(\bm{x},t;s_{2})}{\partial\bm{y}_{j}(\bm{x},t;s_{1})}.$
(3)

It can alternatively be expressed in the following way:

$\bm{F}(\bm{x},t;t,s)=\nabla\bm{y}(\bm{x},t;s).$

Where no confusion can occur, we will omit the particle label and for brevity
write $\bm{F}(s_{1},s_{2})$ for $\bm{F}(\bm{x},t;s_{1},s_{2})$ and $\bm{y}(s)$
for $\bm{y}(\bm{x},t;s)$.

It follows immediately from the definition that $\bm{F}(s,s)=\bm{I}$, where
$\bm{I}$ is the identity tensor, and that
$\bm{F}(s_{1},s_{2})=\bm{F}(s_{2},s_{1})^{-1}$. The deformation gradient is a
two point tensor, and the following lemma describes its evolution with respect
to the two faux time indices:

######  Lemma 3.1 (Evolution of the deformation gradient along
characteristics).

The deformation gradient tensor satisfies the following differential
equations.

$\frac{\partial}{\partial
s_{1}}\bm{F}(\bm{x},t;s_{1},s_{2})=-\bm{F}(\bm{x},t;s_{1},s_{2})\nabla\bm{u}(\bm{y}(\bm{x},t;s_{1}),s_{1})$
(4) $\frac{\partial}{\partial
s_{2}}\bm{F}(\bm{x},t;s_{1},s_{2})=\nabla\bm{u}(\bm{y}(\bm{x},t;s_{2}),s_{2})\bm{F}(\bm{x},t;s_{1},s_{2}).$
(5)

###### Proof.

We begin with Equation (3) and use the chain rule to compute the components of
$\partial_{s_{2}}\bm{F}(s_{1},s_{2})$:

$\begin{split}\frac{\partial}{\partial
s_{2}}\bm{F}(s_{1},s_{2})_{ij}&=\frac{\partial}{\partial
s_{2}}\left(\frac{\partial\bm{y}_{i}(s_{2})}{\partial\bm{y}_{j}(s_{1})}\right)=\frac{\partial\bm{u}_{i}(\bm{y}(s_{2}),s_{2})}{\partial\bm{y}_{j}(s_{1})}=\frac{\partial\bm{u}_{i}(\bm{y}(s_{2}),s_{2})}{\partial\bm{y}_{k}(s_{2})}\frac{\partial\bm{y}_{k}(s_{2})}{\partial\bm{y}_{j}(s_{1})}\\\
&=\nabla\bm{u}(\bm{y}(s_{2}),s_{2})_{ik}\bm{F}(s_{1},s_{2})_{kj},\end{split}$

showing (5). The proof of (4) follows in a similar manner. ∎

######  Proposition 3.2.

The determinant of the deformation gradient $\bm{F}$ is an invariant of
Equations (4) and (5). In particular,
$\det\\!\left({\bm{F}(s_{1},s_{2})}\right)=1$ for all $s_{1},s_{2}$.

###### Proof.

For any $s$, we have
$\det\\!\left({\bm{F}(s,s)}\right)=\det\\!\left({\bm{I}}\right)=1$. Now, (5)
is a tensor-valued differential equation with coefficient matrix
$\nabla\bm{u}(\bm{y}(s_{2}),s_{2})$. Since the fluid is incompressible, the
trace of the velocity gradient is zero. We can therefore invoke [20, §IV.3,
Lemma 3.1] which guarantees $\det\\!\left({\bm{F}(s_{1},s_{2})}\right)$ is
invariant. ∎

###  3.2 The upper convected time derivative

For a generic sufficiently smooth tensor function $\bm{\zeta}$, the upper
convected derivative is defined through

$\overset{\triangledown}{\bm{\zeta}}:=\frac{\partial\bm{\zeta}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{\zeta}-(\nabla\bm{u})\bm{\zeta}-\bm{\zeta}(\nabla\bm{u})^{T}.$

This has a relationship to the notion of a _Lie derivative_ of the tensor
quantity $\bm{\zeta}$ which we will exploit in this work. Indeed, one may
introduce the Lie derivative which follows the fluid with deformation gradient
tensor $\bm{F}$ by

$(\mathcal{L}_{u}\bm{\zeta})(\bm{y}(\bm{x},t;s),s):=\bm{F}(t,s)\frac{\partial}{\partial
s}\left[\bm{F}(s,t)\bm{\zeta}(\bm{y}(s),s)\bm{F}(s,t)^{T}\right]\bm{F}(t,s)^{T}.$

The upper-convected derivative can be inferred from the fundamental properties
of $\bm{F}$ given in §3.1 [25, 30]

$\overset{\triangledown}{\bm{\zeta}}=\lim_{s\to
t}(\mathcal{L}_{u}\bm{\zeta})(\bm{y}(\bm{x},t;s),s).$

The Oldroyd-B constitutive relation for the conformation tensor can therefore
be formulated in terms of the Lie derivative resulting in an equation in one
(characteristic) time variable as follows

$\overset{\triangledown}{\bm{\sigma}}=\lim_{s\to
t}(\mathcal{L}_{u}\bm{\sigma})(\bm{y}(\bm{x},t;s),s).=-\frac{1}{\operatorname{Wi}}\left(\bm{\sigma}-\bm{I}\right).$
(6)

We see that the constitutive law in Oldroyd-B is a linear ODE along the
characteristics induced by the deformation gradient $\bm{F}$. This is why
Oldroyd-B is often referred to as a linear viscoelastic model. It also offers
the opportunity to utilise the a wide variety of timestepping methods designed
for problems of this kind.

##  4 Discretisation of the upper convected time derivative

In this section we introduce a discretisation of the upper convected time
derivative motivated by its formulation as a Lie derivative given in §3. The
method will be similar to characteristic schemes, where the material
derivative is discretised along the characteristics of the flow [37, 41, 44].
Adapting this methodology to account for the deformation of the fluid provides
a natural way to preserve the positive definiteness of the conformation tensor
of any given discretisation. This idea has been explored previously in the
context of several different numerical frameworks [8, 26, 30, 35].

We divide the time domain $[0,T]$ into a uniform partition of $N$ subintervals
of length $\Delta t$ (so that $T=N\Delta t$). We write $t^{n}=n\Delta t$ for
$n=0,1,...,N$, non uniform timestepping is certainly possible however for
clarity of exposition we will restrict to uniform. We also use the notation
$\bm{u}^{n}(\bm{x})=\bm{u}(\bm{x},t^{n})$. Starting from

$(\mathcal{L}_{u}\bm{\zeta})(\bm{y}(\bm{x},t;s),s)=\bm{F}(t,s)\frac{\partial}{\partial
s}\left[\bm{F}(s,t)\bm{\zeta}(\bm{y}(s),s)\bm{F}(s,t)^{T}\right]\bm{F}(t,s)^{T},$
(7)

we will utilise a finite difference approach which motivates the following
definitions:

$\bm{y}^{t}_{\Delta t}(\bm{x})$$(\bm{x},t)$$\bm{y}(\bm{x},t;t-\Delta
t)$$\bm{u}(\bm{x},t)$

Figure 1: Graphical illustration of the departure point
$\bm{y}(\bm{x},t;t-\Delta t)$ of $(\bm{x},t)$ under the flow map, and the
approximate departure point $\bm{y}^{t}_{\Delta t}(\bm{x})$. The curve
represents the characteristic that the particle $(\bm{x},t)$ belongs to.

######  Definition 4.1 (Semidiscrete upper convected derivative).

For $\bm{x}\in\Omega$, $0<\Delta t<T$ and $t\in[\Delta t,T]$, define

$\bm{y}^{t}_{\Delta t}(\bm{x})=\bm{x}-\Delta t\,\bm{u}(\bm{x},t),$ (8)

as an approximation of the departure point of $\bm{x}$ under the fluid
velocity at time $t$. Further, let

$\bm{F}^{t}_{\Delta t}(\bm{x})=\bm{I}+\Delta t\,\nabla\bm{u}(\bm{x},t),$ (9)

denote an approximation of the deformation gradient tensor between times
$t-\Delta t$ and $t$. We define a discrete upper-convected time derivative of
a generic, sufficiently smooth, tensor function $\bm{\zeta}$ by the velocity
field $\bm{u}(\bm{x},t)$ as

$D_{(\bm{u},\,\Delta t)}\bm{\zeta}(\bm{x},t)=\frac{1}{\Delta
t}\Big{(}\bm{\zeta}(\bm{x},t)-\bm{F}^{t}_{\Delta
t}(\bm{x})\bm{\zeta}(\bm{y}^{t}_{\Delta t},t-\Delta t)\bm{F}^{t}_{\Delta
t}(\bm{x})^{T}\Big{)}.$

######  Remark 4.2.

Different choices are available to make this approximation depending on the
desired properties of the scheme. Definition 4.1 results from choosing a
backward Euler discretisation, whereby the right hand side of Equation (7) is
approximated by

$\bm{F}(t,s)\frac{1}{\Delta
t}\Big{\\{}\bm{F}(s,t)\bm{\zeta}(\bm{y}(s),s)\bm{F}(s,t)^{T}-\bm{F}(s-\Delta
t,t)\bm{\zeta}(\bm{y}(s-\Delta t),s-\Delta t)\bm{F}(s-\Delta
t,t)^{T}\Big{\\}}\bm{F}(t,s)^{T},$

which after evaluating at $s=t$ becomes

$\begin{split}\frac{1}{\Delta t}\left\\{\bm{\zeta}(\bm{x},t)-\bm{F}(t-\Delta
t,t)\bm{\zeta}(\bm{y}(t-\Delta t),t-\Delta t)\bm{F}(t-\Delta
t,t)^{T}\right\\}.\end{split}$ (10)

It is not always practical (or possible) to exactly evaluate the deformation
gradient and flow map, so they are replaced with approximations
$\bm{F}^{t}_{\Delta t}$ and $\bm{y}^{t}_{\Delta t}$ respectively in the
definition of $D_{(\bm{u},\,\Delta t)}$.

The structure associated to (6) is dissipative in nature hence a BDF
discretisation is quite natural. A detailed discussion of approximations of
the upper convected time derivative in the manner described above is given in
[30], where the authors discuss first and second order approximations of BDF
type.

We now turn our attention to well-posedness and properties of the scheme. We
first show in Lemma 4.3 that $D_{(\bm{u},\,\Delta t)}$ is well-defined as long
as $\Delta t$ is chosen to be sufficiently small. For a velocity $\bm{u}$
satisfying a homogeneous Dirichlet boundary condition, this result is proved
in [41]. We provide an extension of that result appropriate for nonzero
tangential velocity boundary conditions such as the lid-driven cavity.

######  Lemma 4.3.

Let $\Omega$ be the unit square or the unit cube. Assume that
$\bm{u}\cdot\bm{n}=0$ on $\partial\Omega$, and that $\Delta t$ is chosen so
that $\Delta
t\left\|\bm{u}\right\|_{C^{0}([0,T];\operatorname{W}^{1,\infty}(\Omega))}<1$.
Then for any $\bm{x}\in\Omega$, and any $t\in[\Delta t,T]$, we have
$\bm{y}^{t}_{\Delta t}(\bm{x})\in\Omega$.

###### Proof.

Let $\bm{x}\in\Omega$. Then, by (8)

$|\bm{x}-\bm{y}^{t}_{\Delta t}(\bm{x})|=\Delta t|\bm{u}(\bm{x},t)|.$

The boundary of $\Omega$ has the form $\partial\Omega=\cup_{i=1}\Gamma_{i}$
where $\Gamma_{i}$ is a line segment when $d=2$ and a plane segment when
$d=3$. Let $\bm{x}\in\Omega$ and for some $i$, denote by $\bm{x}_{\Gamma_{i}}$
the unique point on $\Gamma_{i}$ that satisfies
$|\bm{x}-\bm{x}_{\Gamma_{i}}|=\text{dist}(\bm{x},\Gamma_{i})$. As $\Omega$ is
convex this point clearly lies in the interior of $\Gamma_{i}$, so that the
unit outward normal to the boundary at $\bm{x}_{\Gamma_{i}}$ is well-defined,
and the line segment joining $\bm{x}$ and $\bm{x}_{\Gamma_{i}}$ is
perpendicular to $\Gamma_{i}$. We consider the component of the velocity in
the direction of $\bm{n}(\bm{x}_{\Gamma_{i}})$. Since
$\bm{u}(\bm{x}_{\Gamma_{i}},t)\cdot\bm{n}(\bm{x}_{\Gamma_{i}})=0$,

$\begin{split}\bm{u}(\bm{x},t)\cdot\bm{n}(\bm{x}_{\Gamma_{i}})&=\int_{0}^{1}\frac{\partial}{\partial
s}\bm{u}(\bm{x}_{\Gamma_{i}}+s(\bm{x}-\bm{x}_{\Gamma_{i}}),t)\cdot\bm{n}(\bm{x}_{\Gamma_{i}})\mathop{}\\!\mathrm{d}s,\\\
&=(\bm{x}-\bm{x}_{\Gamma_{i}})\cdot\int_{0}^{1}\nabla\bm{u}(\bm{x}_{\Gamma_{i}}+s(\bm{x}-\bm{x}_{\Gamma_{i}}),t)\cdot\bm{n}(\bm{x}_{\Gamma_{i}})\mathop{}\\!\mathrm{d}s\\\
&\leqslant\text{dist}(\bm{x},\Gamma_{i})\left\|\bm{u}\right\|_{C^{0}([0,T];\operatorname{W}^{1,\infty}(\Omega))},\end{split}$

where $\bm{n}(\bm{x}_{\Gamma_{i}})$ is the unit normal to $\partial\Omega$ at
$\bm{x}_{\Gamma_{i}}$. Again, using (8) we have shown that

$|(\bm{x}-\bm{y}^{t}_{\Delta
t}(\bm{x}))\cdot\bm{n}(\bm{x}_{\Gamma_{i}})|\leqslant\text{dist}(\bm{x},\Gamma_{i})\Delta
t\left\|\bm{u}\right\|_{C^{0}([0,T];\operatorname{W}^{1,\infty}(\Omega))}\leqslant\text{dist}(\bm{x},\Gamma_{i}),$

that is, the velocity at $\bm{x}$ in the direction of
$\bm{n}(\bm{x}_{\Gamma_{i}})$ is not large enough for the line segment joining
$\bm{x}$ and $\bm{y}^{t}_{\Delta t}(\bm{x})$ to cross $\Gamma_{i}$. Since the
choice of $i$ was arbitrary, we must have that $\bm{y}^{t}_{\Delta
t}(\bm{x})\in\Omega$. ∎

In the following lemmata, we discuss the approximation properties of
$\bm{y}^{t}_{\Delta t}$, $\bm{F}^{t}_{\Delta t}$ and finally
$D_{(\bm{u},\,\Delta t)}$. We will utilise certain results which were proved
in [30]. The first, Lemma 4.4, gives a bound on the error arising from
approximating the departure point using (8) (see Figure 1 for an
illustration), while Lemma 4.5 bounds the error in evaluating a tensor
function at this approximate departure point as opposed to the true point. The
error in approximating the deformation gradient is bounded in Lemma 4.6, and
finally a bound showing the consistency and approximation order of the
discrete upper convected derivative $D_{(\bm{u},\,\Delta t)}$ is given in
Theorem 4.8.

######  Lemma 4.4 (Approximation of characteristics [30, Lemma 2]).

Suppose that $\bm{u}\cdot\bm{n}=0$ on $\partial\Omega$ and that
$\Delta
t\left\|\bm{u}\right\|_{C^{0}([0,T];\operatorname{W}^{1,\infty}(\Omega))}<1$.
Then for any $\bm{x}\in\Omega$ and $t\geqslant\Delta t$,

$\bm{y}(\bm{x},t;t-\Delta t)-\bm{y}^{t}_{\Delta t}(\bm{x})=\frac{\Delta
t^{2}}{2}\frac{D\bm{u}}{Dt}+\mathcal{O}(\Delta t^{3}).$ (11)

######  Lemma 4.5.

In addition to the hypotheses of Lemma 4.4, let $\bm{\zeta}$ be an
appropriately differentiable tensor function. Then for any $\bm{x}\in\Omega$
and $t\geqslant\Delta t$,

$\bm{\zeta}(\bm{y}(\bm{x},t;t-\Delta t),t-\Delta
t)-\bm{\zeta}(\bm{y}^{t}_{\Delta t}(\bm{x}),t-\Delta t)=\frac{\Delta
t^{2}}{2}\left(\frac{D\bm{u}}{Dt}\cdot\nabla\right)\bm{\zeta}+\mathcal{O}(\Delta
t^{3}).$

######  Lemma 4.6 (Approximation of deformation gradient [30, Lemma 1]).

Suppose that the hypotheses of Lemma 4.4 are satisfied. Then for any
$\bm{x}\in\Omega$ and $t\geqslant\Delta t$,

$\bm{F}(t-\Delta t,t)-\bm{F}^{t}_{\Delta t}(\bm{x})=\Delta
t^{2}\left[(\nabla\bm{u})^{2}-\frac{D(\nabla\bm{u})}{Dt}\right]+\mathcal{O}(\Delta
t^{3}).$

######  Remark 4.7.

The bounds of Lemma 4.4 and Lemma 4.6 require spatial and temporal regularity
of $\bm{u}$ due for example to the appearance of the material derivative of
$\bm{u}$ on the right hand side of Equation (11). In addition, the result of
Lemma 4.5 requires spatial regularity of $\bm{\zeta}$.

We finally arrive at the approximation of the upper convected derivative (see
[30, Remark 7]).

######  Theorem 4.8.

Suppose that the hypotheses of Lemma 4.4 are satisfied, and that $\bm{\zeta}$
is an appropriately differentiable tensor function. Then for any
$\bm{x}\in\Omega$ and $t\geqslant\Delta t$,

$\overset{\triangledown}{\bm{\zeta}}(\bm{x},t)=D_{(\bm{u},\,\Delta
t)}\bm{\zeta}(\bm{x},t)+\mathcal{O}(\Delta t).$

Theorem 4.8 may be interpreted as an expansion around $\bm{x},t$ in that if
follows directly that we have, at $\bm{y}^{t}_{\Delta t}(\bm{x})$,

$\bm{F}^{t}_{\Delta t}(\bm{x})\bm{\zeta}(\bm{y}^{t}_{\Delta t},t-\Delta
t)\bm{F}^{t}_{\Delta t}(\bm{x})^{T}=\bm{\zeta}(\bm{x},t)-\Delta
t\overset{\triangledown}{\bm{\zeta}}(\bm{x},t)+\mathcal{O}(\Delta t^{2}).$

We are now in a position to state the semidiscrete scheme for the constitutive
law (2).

###  4.1 Semidiscretisation of the Oldroyd-B constitutive law

Suppose that $\bm{u}(\bm{x},t)$ is given for $\bm{x}\in\Omega$ and
$t\in[0,T]$. Suppose also that an initial condition $\bm{\sigma}^{0}(\bm{x})$
is given. Then for $n\geqslant 1$ we define $\bm{\Sigma}^{n}$ by

$\begin{split}D_{(\bm{u},\,\Delta
t)}\bm{\Sigma}^{n}(\bm{x})&=-\frac{1}{\operatorname{Wi}}\left(\bm{\Sigma}^{n}(\bm{x})-\bm{I}\right)\\\
\bm{\Sigma}^{0}(\bm{x})&=\bm{\sigma}^{0}(\bm{x}).\end{split}$ (12)

######  Proposition 4.9 (Consistency).

The semidiscrete scheme given in (12) is a consistent, first order
approximation.

###### Proof.

By Theorem 4.8, we have the following representation:

$\bm{\sigma}(\bm{x},t)=\Delta
t\,\overset{\triangledown}{\bm{\sigma}}(\bm{x},t)+\bm{F}^{t}_{\Delta
t}(\bm{x})\bm{\sigma}(\bm{y}^{t}_{\Delta t}(\bm{x}),t-\Delta
t)\bm{F}^{t}_{\Delta t}(\bm{x})^{T}+\mathcal{O}(\Delta t^{2}).$

Then, using Equation (2), we obtain

$\left(1+\frac{\Delta
t}{\operatorname{Wi}}\right)\bm{\sigma}(\bm{x},t)=\bm{F}^{t}_{\Delta
t}(\bm{x})\bm{\sigma}(\bm{y}^{t}_{\Delta t}(\bm{x}),t-\Delta
t)\bm{F}^{t}_{\Delta t}(\bm{x})^{T}-\frac{\Delta
t}{\operatorname{Wi}}\bm{I}+\mathcal{O}(\Delta t^{2}).$ (13)

Assuming that $\bm{\sigma}(\bm{y}^{t}_{\Delta t}(\bm{x}),t-\Delta t)$ is
given, one step of the semidiscrete scheme gives that $\bm{\Sigma}^{n}$
satisfies

$\left(1+\frac{\Delta
t}{\operatorname{Wi}}\right)\bm{\Sigma}^{n}(\bm{x})=\bm{F}^{t}_{\Delta
t}(\bm{x})\bm{\sigma}(\bm{y}^{t}_{\Delta t}(\bm{x}),t-\Delta
t)\bm{F}^{t}_{\Delta t}(\bm{x})^{T}-\frac{\Delta t}{\operatorname{Wi}}\bm{I}.$
(14)

Finally, noting that for any $\Delta t>0$, we have
$0<\tfrac{\operatorname{Wi}}{\operatorname{Wi}+\Delta t}<1$, so that
subtracting (14) from (13) gives

$\bm{\sigma}(\bm{x},t)-\bm{\Sigma}^{n}(\bm{x})=\mathcal{O}(\Delta t^{2}),$

and therefore the semidiscrete scheme has local truncation error proportional
to $\Delta t^{2}$. ∎

######  Proposition 4.10 (The positive definiteness of the conformation tensor
is preserved).

Suppose that

$\bm{\Sigma}^{n-1}$ is given, and is positive definite, and assume that
$\bm{u}$ is given. Then there exists $\Delta t>0$ such that $\bm{\Sigma}^{n}$
is positive definite.

###### Proof.

We can write (12) as

$\left(1+\frac{\Delta
t}{\operatorname{Wi}}\right)\bm{\Sigma}^{n}(\bm{x})=\bm{F}^{t}_{\Delta
t}(\bm{x})\bm{\Sigma}^{n-1}(\bm{y}^{t}_{\Delta t}(\bm{x}))\bm{F}^{t}_{\Delta
t}(\bm{x})^{T}+\frac{\Delta t}{\operatorname{Wi}}\bm{I}.$ (15)

Thus, as long as $\Delta t$ is chosen small enough, a sufficient condition
would be $\Delta t\leqslant\tfrac{1}{2}\left\|\bm{u}\right\|_{1,\infty}$ for
example, $\bm{F}^{t}_{\Delta t}$ is invertible by the Gershgorin circle
theorem [18, Thm 7.2.1]. Let $\bm{0}\neq\bm{\xi}\in\mathbb{R}^{d}$. Then since
$\bm{F}^{t}_{\Delta t}$ is invertible, the product $(\bm{F}^{t}_{\Delta
t})^{T}\bm{\xi}$ is a nonzero vector in $\mathbb{R}^{d}$, so that by
hypothesis

$0<\bm{\xi}^{T}\bm{F}^{t}_{\Delta t}\bm{\Sigma}^{n-1}(\bm{y}^{t}_{\Delta
t}(\bm{x}))(\bm{F}^{t}_{\Delta t})^{T}\bm{\xi}.$

Then from (15), $\bm{\Sigma}^{n}(\bm{x})$ is a sum of two positive definite
terms. ∎

##  5 Spatial discretisation

In this section we present numerical methods for the Oldroyd-B fluid in the
case $\operatorname{Re}=0$. The problem is solved in a decoupled manner, in
which at each time step a Stokes problem is solved for the fluid velocity,
which is then used as a forcing in the time-dependent constitutive relation.
This avoids the need to solve a very large monolithic system and therefore
keeps computational costs down. We first give a finite element method to solve
the Stokes problem and discuss well-posedness. We then present the fully
discrete constitutive law and discuss the solution procedure for the full
problem.

###  5.1 Finite element discretisation for the Stokes problem

We consider $\mathscr{T}{}$ to be a conforming triangulation of $\Omega$,
namely, $\mathscr{T}{}$ is a finite family of sets such that

  1. 1.

$K\in\mathscr{T}{}$ implies $K$ is an open box (i.e. a quadrilateral or
hexahedron),

  2. 2.

for any $K,J\in\mathscr{T}{}$ we have that $\overline{K}\cap\overline{J}$ is a
full lower-dimensional box (i.e., it is either $\emptyset$, a vertex, an edge
(or face) or the whole of $\overline{K}$ and $\overline{J}$) of both
$\overline{K}$ and $\overline{J}$ and

  3. 3.

$\bigcup_{K\in\mathscr{T}{}}\overline{K}=\overline{\Omega}$.

Further, we define $h:\Omega\to\mathbb{R}$ to be the piecewise constant
_meshsize function_ of $\mathscr{T}{}$ given by

$h(\bm{x}):=\max_{\overline{K}\ni\bm{x}}\text{diam}(K).$

With $\mathbb{Q}_{1}(K)$ the space of bilinear polynomials over a
quadrilateral (or trilinear over a hexahedron), we introduce the _finite
element space_

$\mathbb{Q}_{1}:=\\{\phi\in\operatorname{H}^{1}(\Omega):\phi|_{K}\in\mathbb{Q}_{1}(K)\\},$

to be the usual space of continuous piecewise bilinear functions. Let the
vertices of the triangulation $\mathscr{T}{}$ be denoted $\bm{x}_{i}$. For a
vertex $\bm{x}_{i}$, we define

$\hat{\bm{x}}_{i}:=\\{K\in\mathscr{T}{}:\bm{x}_{i}\in\bar{K}\\}.$

To facilitate strong imposition of boundary conditions for the velocity, we
introduce finite element spaces

$\mathbb{V}_{\bm{w}}:=\\{\bm{v}\in\mathbb{Q}_{1}^{d}:\bm{v}\mid_{\partial\Omega}=\Pi_{1}\bm{w}\\}$
$\mathbb{V}_{\bm{0}}:=\\{\bm{v}\in\mathbb{Q}_{1}^{d}:\bm{v}\mid_{\partial\Omega}=\bm{0}\\}.$

To solve the Stokes problem, we employ the stabilised equal-order
approximation introduced in [12]. The stabilisation makes use of the
$\operatorname{L}^{2}$ projection operator
$\Pi_{0}:\mathbb{Q}_{1}\to\mathbb{P}_{0}$ that maps the finite element space
to a space of discontinuous functions which are constant on each element. The
construction of the stabilisation term is therefore local to elements. In
addition, this scheme performs well on graded meshes, which are a useful tool
due to the sharp boundary layers that occur in viscoelastic fluid flows (see
§6.2). A more detailed discussion and analysis for this scheme can be found in
[5, 12, 14]. Let $\bm{\Sigma}_{h}^{0}$ be a given interpolant of
$\bm{\Sigma}^{0}$. Then the discrete problem for the velocity and pressure is
to find $\bm{u}_{h}\in\mathbb{V}_{\bm{w}},p_{h}\in\mathbb{Q}_{1}$ such that

$2\beta\int_{\Omega}\bm{\varepsilon}(\bm{u}_{h}):\bm{\varepsilon}(\bm{v}_{h})-\int_{\Omega}p_{h}\operatorname{div}\bm{v}_{h}-\int_{\Omega}q_{h}\operatorname{div}\bm{u}_{h}-\int_{\Omega}(p_{h}-\Pi_{0}p_{h})(q_{h}-\Pi_{0}q_{h})=\frac{\beta-1}{\operatorname{Wi}}\int_{\Omega}{\bm{\Sigma}_{h}^{0}:\bm{\varepsilon}(\bm{v}_{h})}$
(16)

for all $\bm{v}_{h}\in\mathbb{V}_{\bm{0}},q_{h}\in\mathbb{Q}_{1}$.

######  Proposition 5.1 (Inf-sup stability [5]).

Let

$B_{h}\\!\left({(\bm{u}_{h},p_{h}),(\bm{v}_{h},q_{h})}\right):=2\beta\int_{\Omega}\bm{\varepsilon}(\bm{u}_{h}):\bm{\varepsilon}(\bm{v}_{h})-\int_{\Omega}p_{h}\operatorname{div}\bm{v}_{h}-\int_{\Omega}q_{h}\operatorname{div}\bm{u}_{h}-\int_{\Omega}(p_{h}-\Pi_{0}p_{h})(q_{h}-\Pi_{0}q_{h}).$
(17)

Then, there exists a $C>0$ independent of $h$ such that

$\sup_{\bm{v}_{h}\in\mathbb{Q}_{1}^{d},q_{h}\in\mathbb{Q}_{1}}\frac{B_{h}\\!\left({(\bm{u}_{h},p_{h}),(\bm{v}_{h},q_{h})}\right)}{\left\|\bm{v}_{h}\right\|_{\operatorname{H}^{1}(\Omega)}+\left\|q_{h}\right\|_{\operatorname{L}^{2}(\Omega)}}\geqslant
C\\!\left({\left\|\bm{u}_{h}\right\|_{\operatorname{H}^{1}(\Omega)}+\left\|p_{h}\right\|_{\operatorname{L}^{2}(\Omega)}}\right)\quad\>\forall\>\bm{u}_{h}\in\mathbb{Q}_{1}^{d},p_{h}\in\mathbb{Q}_{1}.$
(18)

The inf-sup condition guarantees that problem (16) is well-posed.

###  5.2 A finite difference approximation of the upper convected time
derivative

The discretisation of the constitutive law requires several definitions
relating to approximations of characteristics and deformation gradients, which
we now give. We write $\bm{u}_{h}^{n}$ for the discrete solution at time
$t^{n}$.

######  Remark 5.2.

If we assume
$\bm{u}\in\operatorname{W}^{1,\infty}(\Omega)\cap\operatorname{H}^{2}(\Omega)$
then standard arguments can be used to infer the discrete velocity
$\bm{u}_{h}^{n}\in\operatorname{W}^{1,\infty}(\Omega)$. Also since
$\bm{u}_{h}^{n}\cdot\bm{n}=0$ on $\partial\Omega$ we have, in light of Lemma
4.3, the approximation (20) is also well-defined.

######  Definition 5.3.

Let $\bm{u}_{h}^{n}\in\mathbb{Q}_{1}^{d}$ be given. For any $\bm{x}\in\Omega$
and $n\geqslant 1$, define the _discrete departure point_ of $\bm{x}$
according to the velocity $\bm{u}_{h}^{n}$ as

$\bm{y}_{h,\Delta t}^{n}(\bm{x})=\bm{x}-\Delta t\bm{u}_{h}^{n}(\bm{x}).$ (19)

We also define the discrete deformation gradient between time $(n-1)\Delta t$
and $n\Delta t$ as

$\bm{F}_{h,\Delta t}^{n}=\bm{I}+\Delta t\nabla\bm{u}_{h}^{n}(\bm{x}).$

######  Definition 5.4 (Fully discrete approximate upper convected time
derivative).

Let $\bm{u}_{h}^{n}$ be given, and suppose that
$\bm{\zeta}^{n-1},\bm{\zeta}^{n}$ are tensor valued functions representing the
the value of a time dependent tensor field at the time levels $n-1,n$. Then
the fully discrete approximation to the upper convected time derivative of
$\bm{\zeta}$ at time level $n$ is given by

$D^{h}_{(\bm{u}_{h}^{n},\,\Delta t)}\bm{\zeta}^{n}(\bm{x}):=\frac{1}{\Delta
t}\Big{(}\bm{\zeta}^{n}(\bm{x})-\bm{F}_{h,\Delta
t}^{n}(\bm{x})\bm{\zeta}^{n-1}(\bm{y}^{n}_{h,\Delta t})\bm{F}_{h,\Delta
t}^{n}(\bm{x})^{T}\Big{)}.$ (20)

###  5.3 A fully discrete scheme for the Stokes flow of an Oldroyd-B fluid

We are now in a position to define the full finite difference discretisation
of the constitutive law. Values of the discrete conformation tensor are stored
at the vertices $\bm{x}_{i}$. At any time level $k$, these values are
sufficient to uniquely define a function
$\bm{\Sigma}^{k}_{h}\in\mathbb{Q}_{1}^{d(d+1)/2}$ by bilinear (or trilinear)
interpolation. This representation is used to evaluate the conformation tensor
at points which are not vertices, and is the function on the right hand side
of the discrete Stokes equation (16). By a slight abuse of notation, we will
write $\bm{\Sigma}_{h}^{k}$ for both the set of point values at vertices, and
the function in $\mathbb{Q}_{1}^{d(d+1)/2}$ obtained by interpolation using
these values.

We note that the approximation of the deformation gradient needed to evaluate
(20) requires the value of the discrete velocity gradient, which we wish to
use at the vertices of the mesh. Since the discrete velocity lies in
$\mathbb{Q}_{1}^{d}$, the point values of its gradient are not well defined.
For a vertex $\bm{x}_{i}$ of the mesh and a function
$\bm{v}_{h}\in\mathbb{Q}_{1}^{d}$, in place of the velocity gradient, we
therefore use the local average over $\hat{\bm{x}}_{i}$:

$\nabla\bm{v}_{h}(\bm{x}_{i})=\frac{1}{|\hat{\bm{x}}_{i}|}\sum_{K\subseteq\hat{\bm{x}}_{i}}\int_{K}\nabla\bm{v}_{h}.$

The fully discrete problem is given as follows. Let $\bm{\Sigma}^{0}_{h}$ be
an interpolant of the initial condition. Then for $1\leqslant n\leqslant N$,
find $\bm{u}^{n}_{h}\in\mathbb{V}_{\bm{w}}$, $p_{h}^{n}\in\mathbb{Q}_{1}$ such
that

$B_{h}\\!\left({(\bm{u}^{n}_{h},p^{n}_{h}),(\bm{v}_{h},q_{h})}\right)=\frac{\beta-1}{\operatorname{Wi}}\int_{\Omega}{\bm{\Sigma}^{n-1}_{h}:\bm{\varepsilon}(\bm{v}_{h})}$
(21)

for all $\bm{v}_{h}\in\mathbb{V}_{\bm{0}},q_{h}\in\mathbb{Q}_{1}$, then
$\bm{\Sigma}_{h}^{n}$ is given by

$D^{h}_{(\bm{u}_{h}^{n},\,\Delta
t)}\bm{\Sigma}_{h}^{n}(\bm{x}_{i})+\frac{1}{\operatorname{Wi}}(\bm{\Sigma}^{n}_{h}(\bm{x}_{i})-\bm{I})=0.$
(22)

At each mesh point, (22) is a coupled system of ordinary differential
equations forced by the discrete fluid velocity via the deformation gradient,
and approximated in an explicit manner, and therefore the nodal values may be
updated without the need to solve a linear system.

######  Remark 5.5 (Choice of spatial discretisation points for finite
difference scheme).

We note that since all derivatives appearing in the constitutive relation are
approximated via (10), we are free to choose any set of points at which to
approximate the conformation tensor, with the proviso that care must be taken
to adequately resolve the velocity gradient, as this is what forces the
constitutive law. The choice of these points so that they correspond with the
degrees of freedom of a piecewise bilinear finite element space is a natural
one, as it is then clear how to include the discrete conformation tensor in
the momentum equation. Other possible choices are discussed in §7.

######  Proposition 5.6.

If the initial condition $\bm{\sigma}^{0}$ is positive definite, and bilinear
interpolation is used to obtain $\bm{\Sigma}_{h}^{0}$ and to evaluate
$\bm{\Sigma}_{h}^{k}$ at points which are not vertices of the mesh, then at
each time step there exists $\Delta t>0$ such that the discrete conformation
tensor remains positive definite.

###### Proof.

By assumption, all nodal values of $\bm{\Sigma}_{h}^{n-1}$ are positive
definite. Then since interpolated values are weighted sums of the nodal values
with non-negative weights, the interpolant is also positive definite. By the
same argument as in Proposition 4.10, if we choose $\Delta
t\leqslant\tfrac{1}{2}\left\|\bm{u}_{h}^{n}\right\|_{1,\infty}$, and if
$\bm{\Sigma}_{h}^{n-1}$ is positive definite, then $\bm{\Sigma}_{h}^{n}$
computed from the scheme (22) will also be positive definite. ∎

##  6 Numerical Experiments

In this section, we present numerical results obtained by solving (21)-(22)
for a variety of model and discretisation parameters. All simulations
presented here were conducted using deal.II, an open source C++ software
library providing tools for finite element computations [2].

###  6.1 Description of the test case

The lid-driven cavity consists of a fluid-filled, impenetrable box in which
the flow is driven by a moving lid. The flow of a Newtonian fluid in a lid-
driven cavity is by now a very well studied problem. It is known that the flow
characteristics, see Figure 2, depend upon the Reynolds number. Studies from
the numerical literature largely agree with one another qualitatively and
quantitatively for moderate Reynolds numbers. For larger Reynolds numbers,
numerical solutions become more challenging, and differences appear in
computational results. Detailed numerical studies have provided benchmarks for
Reynolds numbers in the tens of thousands [9, 17].

In experimental studies, for small Reynolds number, flow in the cavity reaches
a steady state, and, in a three dimensional cavity, remains approximately two
dimensional. For the setup in Figure 2, i.e., with lid moving left to right,
as the Reynolds number increases the centre of the main vortex moves to the
right. Above a critical value of the Reynolds number, the flow transitions to
a time-periodic state with variation in the neutral direction [1]. Even for
inertialess non-Newtonian flows, the symmetry of the flow is broken (see
Figure 2, right) and moves towards the upper left corner with increasing
Weissenberg number [33].

Moving LidClockwise Circular Flow

Figure 2: An illustration of a lid driven cavity setup (left). Streamlines of
a Newtonian Stokes flow (middle) with left-right symmetry in the velocity
magnitude. Introducing a polymeric stress breaks this symmetry, with the
centre of circulation moving up and to the left with increasing
$\operatorname{Wi}$ (right).

In this work, the regularised lid-driven cavity problem was solved, for which
$\Omega=[0,1]\times[0,1]$, $\bm{u}=(0,0)^{T}$ on the lateral and lower
boundaries, and on the lid (i.e. $\\{(x_{1},x_{2}):x_{2}=1\\}$), the velocity
boundary condition is

$\bm{u}_{lid}(\bm{x},t)=8\left(1+\tanh\left(8\left(t-\tfrac{1}{2}\right)\right)x_{1}^{2}(1-x_{1}^{2}),0\right)^{T}.$

######  Remark 6.1.

For this problem, a Weissenberg number can be defined as

$\operatorname{Wi}=\frac{\lambda U}{H},$

where $\lambda$ is a characteristic relaxation time of the fluid, $U$ is a
characteristic velocity scale and $H$ is the height of the cavity. Thus, $U/H$
is a measure of the shear rate [33]. For the setup described above, $U=1$,
$H=1$.

###  6.2 Mesh design

Mesh design and selection of discretisation parameters appears to be of
crucial importance for the cavity flow of an Oldroyd-B fluid, and indeed very
few studies have reported satisfactory numerical results for the lid-driven
cavity using a uniform mesh, with the notable exception of schemes which
utilise the log-conformation tensor formulation. In addition, results for some
schemes have been shown to be very sensitive to mesh size, for example in [7],
dissipation properties, positive definiteness of the conformation tensor and
existence of steady state change qualitatively upon varying the mesh size.

For the cavity problem, graded meshes of various types are a common choice.
Meshes constructed in this manner are able to provide very fine resolution
near the boundary (particularly the lid), where it is most needed, in an
efficient manner. The situation may be compared to singularly perturbed
convection-diffusion problems, which exhibit boundary layers that put severe
limitations on the mesh size. To make computations tractable, layer-adapted
meshes have been extensively researched in this community (see [27] for a
review). A particularly effective example is the Bakhvalov-type mesh, where
upper bounds on the derivatives of the solution near the boundary are used to
design mesh grading functions that are sufficient to resolve it.

In this work, we will use a mixture of the following. One construction
consists of decreasing element size towards the boundary with a constant
contraction ratio $\gamma$, i.e. $h_{i+1}=\gamma h_{i}$ where the index $i$ is
such that $i+1$ is the next element in the direction of the boundary [10, 43].
Another is to decrease element size by a constant amount $c$ so that
$h_{i+1}=h_{i}-c$ [34]. For comparison, the latter is equivalent to decreasing
the contraction ratio, so that elements size decreases very rapidly near the
boundary. Finally, in [19], mixtures of uniform grid spacing and refined areas
are used.

Figure 3:  Illustrations of the two types of graded mesh used. Both meshes
consist of 1024 elements. Left: mesh $\mathcal{R}_{32}$, offering aggressive
refinement near the boundary. Right: $\mathscr{T}{}_{32}$. Bottom row: zoom of
upper left corners showing an area of $0.08\times 0.08$ to illustrate the
differences in grading.

Meshes constructed from central square elements with a constant reduction
ratio of 0.96 towards each boundary are denoted by $\mathscr{T}{}_{N}$, where
$N$ is the number of elements in each coordinate direction. We denote by
$\mathcal{R}_{N}$ meshes constructed with linearly decreasing element size as
follows (similar to those used in [34]). We define vertices for an $N\times N$
quadrilateral mesh by their $x$ and $y$ coordinates, given respectively by
$x_{i}=2\left(\tfrac{i}{N}\right)^{2}$ for $i=0,1,...,\tfrac{N}{2}$, with
$x_{i}=1-x_{N-i}$ for $i=\tfrac{N}{2}+1,...,N$, and
$y_{i}=1-\left(1-\tfrac{i}{N}\right)^{2}$ for $i=0,...,N$. On these meshes,
the finest vertical resolution occurs at the lid, with an element height of
$\tfrac{1}{N^{2}}$.

This means the meshes are highly anisotropic away from the corners. Stabilised
finite element methods often have stabilisation parameters that depend badly
on the anisotropy of the elements, however the solver we used (16) is
parameter free.

Table 1: Mesh statistics Mesh |  # elements |  $h_{\text{min}}$ |  $h_{\text{max}}$
---|---|---|---
$\mathscr{T}{}_{90}$ |  8,100 |  0.0039 |  0.024
$\mathscr{T}{}_{120}$ |  14,400 |  0.0020 |  0.022
$\mathscr{T}{}_{150}$ |  22,500 |  0.0010 |  0.021
$\mathscr{T}{}_{180}$ |  32,400 |  0.00054 |  0.021
$\mathcal{R}_{256}$ |  65,536 |  $1.5\times 10^{-5}$ |  $7.8\times 10^{-3}$

###  6.3 Numerical results

Where possible, the results presented here are compared with the literature,
particularly [16, 34, 43] \- see Table 2. Descriptions of the meshes used are
given in §6.2, and a summary of the key statistics of each is given in Table
1. The data for comparison are the logarithm of the maximum value of
$\bm{\sigma}_{1}1$ along the midline, $x=0.5$, the maximum value of
$\bm{\sigma}_{1}1$ over the domain, and the coordinates of the centre of the
main vortex, $x_{c},y_{c}$.

####  6.3.1 $\operatorname{Wi}=0.5$

The solution appears to have reached a steady state by $t=8$. In Figure 4,
components of the conformation tensor are plotted along cross sections of the
domain for $\operatorname{Wi}=0.5$ and for various computational meshes
described in §6.2 and Table 1. The solutions obtained for $\bm{\sigma}_{11}$
and $\bm{\sigma}_{22}$ agree qualitatively very well with Figure 5 of [34],
albeit with differences in magnitude along the lid. We note that the meshes
used there are finer, which could account for the difference. Other metrics
available for comparison from the numerical literature (see Table 4 in [43]
and Table 2) are the centre of the main vortex in the velocity field and the
logarithm of the maximum value of $\bm{\sigma}_{11}$ attained along the line
$x=0.5$. Computed on mesh $\mathcal{R}_{256}$, the former agrees very closely
(within 1%) although not much variation in prediction of this quantity is
observed across the numerical literature. The latter is approximately 2%
larger than the other predictions given in [43]. Resolution of the layer at
the top boundary appears to be of critical importance for the accuracy of the
approximation. Under-resolution results in altered behaviour of
$\bm{\sigma}_{12}$ at the upper boundary. In addition, coarse meshes appear to
over-estimate $\bm{\sigma}_{11}$. These errors do not appear to pollute the
solution away from the boundary.

Figure 4:  Components of the conformation tensor plotted along cross sections of the domain, obtained on meshes $\mathscr{T}{}_{90}$, $\mathscr{T}{}_{120}$, $\mathscr{T}{}_{150}$ and $\mathcal{R}_{256}$ with $\operatorname{Wi}=0.5$ and $t=10$. Top row: plots over the top boundary given by $y=1$. Middle row: plots over the line given by $y=0.75$. Bottom row: plots over the midline given by $x=0.5$.  Table 2: Comparison of results with published data for the (regularised) lid-driven cavity, $\beta=0.5$. Value marked ∗ was estimated from [34] using their predicted relationship between this quantity and $\operatorname{Wi}$ (see [34, §4.3]), while values marked † were estimated using a plot digitisation tool [38]. All other values were reported directly. All values from the current work were taken after a steady state had been reached (approx $t=10$ for $\operatorname{Wi}=0.5$ and $t=30$ for $\operatorname{Wi}=1$). Reference |  $\underset{x=0.5}{\max}\ln(\bm{\sigma}_{11})$ |  $\max\bm{\sigma}_{11}$ |  $x_{c},y_{c}$
---|---|---|---
$\operatorname{Wi}=0.5$ |  |  |
Current work $\mathscr{T}{}_{90}$ |  5.76 |  351.21 |  0.466, 0.799
Current work $\mathscr{T}{}_{120}$ |  5.65 |  309.65 |  0.466, 0.799
Current work $\mathscr{T}{}_{150}$ |  5.60 |  295.04 |  0.467, 0.799
Current work $\mathcal{R}_{256}$ |  5.59 |  290.59 |  0.467, 0.799
Pan et al. [34] |  $5.59^{\dagger}$ |  $289^{*}$ |  0.469, 0.798
Sousa et al. [43] M4 |  5.51 |  - |  0.466, 0.800
$\operatorname{Wi}=1$ |  |  |
Current work $\mathscr{T}{}_{180}$ |  10.75 |  48,038.9 |  0.428, 0.819
Current work $\mathcal{R}_{256}$ |  10.22 |  28,069.2 |  0.431, 0.819
Pan et al. [34] |  $9.35^{\dagger}$ |  11,529.43 |  0.439, 0.816
Sousa et al. [43] M4 |  7.80 |  - |  0.434, 0.816

####  6.3.2 $\operatorname{Wi}=1$

The solution appears to have reached a steady state by $t=27.5$. The
conformation tensor field is plotted component-wise in Figure 5, with
benchmark metrics reported in Table 2. The logarithm of the conformation
tensor was computed to allow comparison with other published works, with
selected cross sections presented in Figure 6. A sharp boundary layer at the
lid is observed in $\bm{\sigma}_{11}$, with all components exhibiting large
gradients near the upper corners of the domain. Our computational results
agree qualitatively with those presented in [16, Figure 4], although in our
simulations a steady state was reached much later than $t=8$, as reported
there. In addition, the maximum value $\bm{\sigma}_{11}$ along $x=0.5$ was
significantly overestimated in our work relative to other published data, but
there is significant variation in this figure, with [34] reporting a value two
logarithmic orders larger than [16]. However, it can be seen from Figure 6
that large discrepancies occur only in a very small neighbourhood of the upper
boundary. Qualitative differences are observed only in the upper right hand
corner.

Figure 5:  $\operatorname{Wi}=1$. Components of the conformation tensor at
$t=30$: $\bm{\sigma}_{11}$, $\bm{\sigma}_{12}$, $\bm{\sigma}_{22}$ (note
$\bm{\sigma}_{11}$ and $\bm{\sigma}_{22}$ are logarithmically coloured and
contoured while $\bm{\sigma}_{12}$ is linearly coloured and contoured).
Different colormaps are used to emphasise the difference in magnitude between
the components. Minimum and maximum values are shown on the colour bar.

Figure 6:  Components of the logarithm of the conformation tensor plotted
along cross sections of the domain obtained on the mesh $\mathcal{R}_{256}$
with $\operatorname{Wi}=1$ (solid line). Dotted line is numerical data from
Fattal & Kupferman (2005) [16], while the dashed line is data from Pan, Hao &
Glowinski (2009) [34].

In Figure 7, components of the conformation tensor are plotted along cross
sections for two Weissenberg numbers, 0.5 and 1, to illustrate the difference
in solutions. Larger magnitudes and exaggerated features are observed in all
components for $\operatorname{Wi}=1$ when compared with
$\operatorname{Wi}=0.5$.

Figure 7:  Components of the conformation tensor plotted along cross sections
of the domain obtained on the mesh $\mathcal{R}_{256}$ with
$\operatorname{Wi}=0.5,1$. Top row: plots over the top boundary given by
$y=1$. Middle row: plots over the line given by $y=0.75$. Bottom row: plots
over the midline given by $x=0.5$.

###  6.4 Positive definiteness of the discrete conformation tensor

A major aim of the numerical scheme presented is to maintain a positive
definite comformation tensor, both to respect the physics of the problem and
avoid numerical instabilities that can arise when this property is lost. We
recall from Proposition 5.6 that this is guaranteed if the time step is chosen
to be sufficiently small. To ensure that an appropriate choice was made, the
eigenvalues of the discrete conformation tensor were computed, and indeed at
all times, the smallest eigenvalue of the conformation tensor was strictly
positive, and therefore the discrete conformation tensor remained positive
definite throughout the computation. The eigenvalues of the conformation
tensor are visualised in Figure 8.

Figure 8:  Left: Smallest eigenvalue of the discrete conformation tensor.
Right: largest eigenvalue. Both plots are logarithmically coloured and
contoured. Minimum and maximum values are shown on the colour bar.

##  7 Conclusions & Discussion

Results from simulations of an Oldroyd-B fluid in the creeping flow regime
using a discretisation of the Lie derivative (7) were presented. The scheme
maintained positive definiteness of the discrete conformation tensor and
achieved good agreement with existing published data. Our proposed finite
difference method circumvents quadrature issues, representing advection
effectively compared to some finite element schemes, which can under-represent
advection for small time steps or coarse quadrature formulas [8].

Qualitatively and quantitatively, the results are consistent with those in the
literature using the log-conformation representation [16, 34, 43]. The
scheme’s update process for the conformation tensor is efficient and
straightforward to implement. However, spatial resolution and mesh design
significantly influence the results, requiring fine meshes for convergence.

Further work includes exploring enhancements to preserve additional fluid
structures, particularly incompressibility constraints and investigating
differential equations on the special linear group for a consistent discrete
deformation gradient. This may involve combining the approach with divergence-
free finite element methods for fluid equations or exploring finer meshes or
different finite difference points for the constitutive law to address
challenges in the creeping flow regime.

##  8 Acknowledgements

This work has been partially supported by the Leverhulme Trust Research
Project Grant No. RPG-2021-238. TP is also partially supported by EPRSC grants
EP/W026899/2, EP/X017206/1 and EP/X030067/1. The authors also want to thank
Gabriel Barrenechea and Emmanuil Geourgoulis for helpful discussions and
suggestions.

## References

  * [1] Cyrus K Aidun, NG Triantafillopoulos and JD Benson  “Global stability of a lid-driven cavity with throughflow: Flow visualization studies”  In _Physics of Fluids A: Fluid Dynamics_ 3.9 American Institute of Physics, 1991, pp. 2081–2091
  * [2] Daniel Arndt et al.  “The deal.II finite element library: Design, features, and insights”  In _Computers & Mathematics with Applications_ 81, 2021, pp. 407–422  DOI: 10.1016/j.camwa.2020.02.022
  * [3] Mohammed Bensaada, Driss Esselaoui and Pierre Saramito  “Error estimate for the characteristic method involving Oldroyd derivative in a tensorial transport problem”  In _Electronic Journal of Differential Equations_ , 2004
  * [4] Robert Byron Bird, Robert Calvin Armstrong and Ole Hassager  “Dynamics of polymeric liquids. Vol. 1: Fluid mechanics”  John WileySons Inc., New York, NY, 1987
  * [5] Pavel B Bochev, Clark R Dohrmann and Max D Gunzburger  “Stabilization of low-order mixed finite elements for the Stokes equations”  In _SIAM Journal on Numerical Analysis_ 44.1 SIAM, 2006, pp. 82–101
  * [6] J Bonvin, M Picasso and R Stenberg  “GLS and EVSS methods for a three-field Stokes problem arising from viscoelastic flows”  In _Computer Methods in Applied Mechanics and Engineering_ 190.29-30 Elsevier, 2001, pp. 3893–3914
  * [7] Sébastien Boyaval  “Lid-driven-cavity simulations of Oldroyd-B models using free-energy-dissipative schemes”  In _Numerical Mathematics and Advanced Applications 2009: Proceedings of ENUMATH 2009, the 8th European Conference on Numerical Mathematics and Advanced Applications, Uppsala, July 2009_ , 2010, pp. 191–198  Springer
  * [8] Sébastien Boyaval, Tony Lelièvre and Claude Mangoubi  “Free-energy-dissipative schemes for the Oldroyd-B model”  In _ESAIM: Mathematical Modelling and Numerical Analysis_ 43.3 EDP Sciences, 2009, pp. 523–561
  * [9] Charles-Henri Bruneau and Mazen Saad  “The 2D lid-driven cavity problem revisited”  In _Computers & fluids_ 35.3 Elsevier, 2006, pp. 326–348
  * [10] Raphaël Comminal, Jon Spangenberg and Jesper Henri Hattel  “Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation”  In _Journal of Non-Newtonian Fluid Mechanics_ 223 Elsevier, 2015, pp. 37–61
  * [11] JM Dealy  “Weissenberg and Deborah numbers—their definition and use”  In _Rheol. Bull_ 79.2, 2010, pp. 14–18
  * [12] Clark R Dohrmann and Pavel B Bochev  “A stabilized finite element method for the Stokes problem based on polynomial pressure projections”  In _International Journal for Numerical Methods in Fluids_ 46.2 Wiley Online Library, 2004, pp. 183–201
  * [13] F Dupret and JM Marchal  “Sur le signe des valeurs propres du tenseur des extra-contraintes dans un écoulement de fluide de Maxwell”  In _Journal de mécanique théorique et appliquée_ 5.3, 1986, pp. 403–427
  * [14] Howard C Elman, David J Silvester and Andrew J Wathen  “Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics”  Oxford university press, 2014
  * [15] Raanan Fattal and Raz Kupferman  “Constitutive laws for the matrix-logarithm of the conformation tensor”  In _Journal of Non-Newtonian Fluid Mechanics_ 123.2-3 Elsevier, 2004, pp. 281–285
  * [16] Raanan Fattal and Raz Kupferman  “Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation”  In _Journal of Non-Newtonian Fluid Mechanics_ 126.1 Elsevier, 2005, pp. 23–37
  * [17] UKNG Ghia, Kirti N Ghia and CT Shin  “High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method”  In _Journal of computational physics_ 48.3 Elsevier, 1982, pp. 387–411
  * [18] Gene H Golub and Charles F Van Loan  “Matrix computations”  JHU press, 2013
  * [19] Florian Habla, Ming Wei Tan, Johannes Haßlberger and Olaf Hinrichsen  “Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM®”  In _Journal of Non-Newtonian Fluid Mechanics_ 212 Elsevier, 2014, pp. 47–62
  * [20] Ernst Hairer, Marlis Hochbruck, Arieh Iserles and Christian Lubich  “Geometric numerical integration”  In _Oberwolfach Reports_ 3.1, 2006, pp. 805–882
  * [21] John Hinch and Oliver Harlen  “Oldroyd B, and not A?”  In _Journal of Non-Newtonian Fluid Mechanics_ 298 Elsevier, 2021, pp. 104668
  * [22] Thomas JR Hughes  “Numerical implementation of constitutive models: rate-independent deviatoric plasticity”  In _Theoretical foundation for large-scale computations for nonlinear material behavior: Proceedings of the Workshop on the Theoretical Foundation for Large-Scale Computations of Nonlinear Material Behavior Evanston, Illinois, October 24, 25, and 26, 1983_ , 1984, pp. 29–63  Springer
  * [23] Martien A Hulsen, Raanan Fattal and Raz Kupferman  “Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms”  In _Journal of Non-Newtonian Fluid Mechanics_ 127.1 Elsevier, 2005, pp. 27–39
  * [24] Azadeh Jafari, Alireza Chitsaz, Reza Nouri and Timothy N Phillips  “Property preserving reformulation of constitutive laws for the conformation tensor”  In _Theoretical and Computational Fluid Dynamics_ 32 Springer, 2018, pp. 789–803
  * [25] Young-Ju Lee and Jinchao Xu  “New formulations, positivity preserving discretizations and stability analysis for non-Newtonian flow models”  In _Computer methods in applied mechanics and engineering_ 195.9-12 Elsevier, 2006, pp. 1180–1206
  * [26] Young-Ju Lee, Jinchao Xu and Chen-Song Zhang  “Stable finite element discretizations for viscoelastic flow models”  In _Handbook of numerical analysis_ 16 Elsevier, 2011, pp. 371–432
  * [27] Torsten Linß  “Layer-adapted meshes for convection–diffusion problems”  In _Computer Methods in Applied Mechanics and Engineering_ 192.9-10 Elsevier, 2003, pp. 1061–1105
  * [28] Alexei Lozinski and Robert G Owens  “An energy estimate for the Oldroyd B model: theory and applications”  In _Journal of non-newtonian fluid mechanics_ 112.2-3 Elsevier, 2003, pp. 161–176
  * [29] Mária Lukáčová-Medvid’ová, Hirofumi Notsu and Bangwei She  “Energy dissipative characteristic schemes for the diffusive Oldroyd-B viscoelastic fluid”  In _International Journal for Numerical Methods in Fluids_ 81.9 Wiley Online Library, 2016, pp. 523–557
  * [30] Debora D Medeiros, Hirofumi Notsu and Cassio M Oishi  “Second-order finite difference approximations of the upper-convected time derivative”  In _SIAM Journal on Numerical Analysis_ 59.6 SIAM, 2021, pp. 2955–2988
  * [31] James G Oldroyd  “On the formulation of rheological equations of state”  In _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 200.1063 The Royal Society London, 1950, pp. 523–541
  * [32] Peyman Pakdel and Gareth H McKinley  “Cavity flows of elastic liquids: purely elastic instabilities”  In _Physics of Fluids_ 10.5 American Institute of Physics, 1998, pp. 1058–1070
  * [33] Peyman Pakdel, Stephen H Spiegelberg and Gareth H McKinley  “Cavity flows of elastic liquids: two-dimensional flows”  In _Physics of Fluids_ 9.11 American Institute of Physics, 1997, pp. 3123–3140
  * [34] Tsorng-Whay Pan, Jian Hao and Roland Glowinski  “On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid”  In _International Journal for Numerical Methods in Fluids_ 60.7 Wiley Online Library, 2009, pp. 791–808
  * [35] Jerzy Petera  “A new finite element scheme using the Lagrangian framework for simulation of viscoelastic fluid flows”  In _Journal of non-newtonian fluid mechanics_ 103.1 Elsevier, 2002, pp. 1–43
  * [36] F Pimenta and MA Alves  “Stabilization of an open-source finite-volume solver for viscoelastic fluid flows”  In _Journal of Non-Newtonian Fluid Mechanics_ 239 Elsevier, 2017, pp. 85–104
  * [37] Olivier Pironneau  “On the transport-diffusion algorithm and its applications to the Navier-Stokes equations”  In _Numerische Mathematik_ 38 Springer, 1982, pp. 309–332
  * [38] “PlotDigitizer: Version 3.1.5”, 2023  URL: https://plotdigitizer.com
  * [39] Michael Renardy  “Mathematical analysis of viscoelastic flows”  SIAM, 2000
  * [40] Michael Renardy and Becca Thomases  “A mathematician’s perspective on the Oldroyd B model: Progress and future challenges”  In _Journal of Non-Newtonian Fluid Mechanics_ 293 Elsevier, 2021, pp. 104573
  * [41] Hongxing Rui and Masahisa Tabata  “A second order characteristic finite element scheme for convection-diffusion problems”  In _Numerische Mathematik_ 92.1 Springer, 2002, pp. 161–177
  * [42] IM Rutkevich  “The propagation of small perturbations in a viscoelastic fluid”  In _Journal of Applied Mathematics and Mechanics_ 34.1 Elsevier, 1970, pp. 35–50
  * [43] RG Sousa, RJ Poole, AM Afonso, FT Pinho, PJ Oliveira, A Morozov and MA Alves  “Lid-driven cavity flow of viscoelastic liquids”  In _Journal of Non-Newtonian Fluid Mechanics_ 234 Elsevier, 2016, pp. 129–138
  * [44] Endre Süli  “Convergence and nonlinear stability of the Lagrange-Galerkin method for the Navier-Stokes equations”  In _Numerische Mathematik_ 53 Springer, 1988, pp. 459–483
  * [45] Wen-Hua Zhang, Jingfa Li, Qiankun Wang, Yu Ma, Hong-Na Zhang, Bo Yu and Fengchen Li  “Comparative study on numerical performances of log-conformation representation and standard conformation representation in the simulation of viscoelastic fluid turbulent drag-reducing channel flow”  In _Physics of Fluids_ 33.2 AIP Publishing, 2021


Engaging Students Through Math Competitions

Béla Bajnok111Béla Bajnok is a Professor of Mathematics at Gettysburg College
and is the Director of the American Mathematics Competitions program of the
MAA. His email address is<EMAIL_ADDRESS>
I was fortunate to grow up in Hungary, a country with a long and distinguished
history in mathematical competitions. Like many of my friends, I was keenly
involved in these contests. We liked training for them, being exposed to one
beautiful problem after another; we enjoyed participating in them, even when
we didn’t end up winning; and we loved being in a community that was welcoming
yet tight-knit. I have remained involved with mathematical competitions ever
since, and thus was excited and honored when, in 2017, I was asked to become
the director of the American Mathematics Competitions (AMC) program of the
Mathematical Association of America (MAA). While the winners of our
competitions – especially those participating in the International
Mathematical Olympiad (IMO) – receive much recognition, our goals are to
encourage students to engage in mathematics nationwide, and to discover,
develop, and nurture talent more broadly.

With the participation of over 300,000 students each year, the AMC is the
largest program of the MAA. The competitions start with the AMC 8, the AMC 10,
and the AMC 12 exams, open to students in grade 8 or below, grade 10 or below,
and grade 12 or below, respectively. Based on their performance on these
multiple-choice competitions, approximately 6,000 students are invited to take
the American Invitational Mathematics Exam (AIME), a challenging three-hour
exam where the answer to each of the 15 problems is a nonnegative integer
under 1000. The competition series culminates with the USA Mathematical
Olympiad (USAMO) and the USA Junior Mathematical Olympiad (USAJMO), offered to
approximately 500 students. These Olympiads follow the style of the IMO: they
ask students to provide rigorous proofs for three problems on each of two
consecutive days, with an allowed time of four and a half hours each day. We
hope that our competitions show students – and their teachers and society at
large – the beauty and power of mathematics, usually well beyond what they see
in typical mathematics classrooms.

The creation of problems suitable for the AMC is a highly challenging task.
Even at the beginning levels where problems stay close to the standard school
curriculum, we aim for problems that ask questions in novel and fun ways.
Coming up with beautiful – yet still elementary – problems at the Olympiad
level is particularly challenging. Yet, year after year, collections of
ingenious and captivating problems are created; everyone with an interest is
encouraged to explore them. Beyond the competitions themselves, we hope that
our problems provide learning opportunities for current and future students,
their teachers, and anyone with an interest in mathematics. We present three
examples here: one from an AMC 8 exam, one from an AIME competition, and one
from a USA Mathematical Olympiad.

2022 AMC 8 Problem 25 (created by Silva Chang): A cricket hops between four
leaves, on each turn hopping to one of the other three leaves with equal
probability. After four hops what is the probability that the cricket returns
to the leaf where it started?

Figure 1: Four leaves and a cricket

There is a large variety of different techniques to solve this problem: the
official solution guide lists four of them, from case counting to recursions,
and additional methods such as dynamic programming or generating functions can
be used as well. Problems being suitable for multiple approaches are of course
helpful for the contestants at the time, but this subsequently also allows
teachers to discuss techniques that are new to their students.

2022 AIME I Problem 1 (created by David Altizio): Quadratic polynomials $P(x)$
and $Q(x)$ have leading coefficients $2$ and $-2,$ respectively. The graphs of
both polynomials pass through the two points $(16,54)$ and $(20,53).$ Find
$P(0)+Q(0).$

This problem can, of course, be solved by finding the two quadratic
polynomials that satisfy the conditions of this problem, but this would take a
while, and the answers are not particularly nice. The clever approach (perhaps
suggested by the question itself) would focus on the polynomial $P(x)+Q(x)$,
which is linear; since its graph goes through the points $(16,108)$ and
$(20,106)$, it equals $-x/2+116$, and thus $P(0)+Q(0)=116$.

One of the challenges we have when assembling our exams is to be able to
select the students who will then move on to subsequent competitions and win
awards and prizes, while also providing approachable problems to all
participants. Three hours may seem long for the AIME competition, but time is
an important factor for many participants, so being able to solve the easiest
problems quickly is a big advantage.

2021 USAMO Problem 3 (created by Shaunak Kishore and Alex Zhai): Let $n\geq 2$
be an integer. An $n$-by-$n$ board is initially empty. Each minute, you may
perform one of two moves:
• If there is an L-shaped tromino region of three cells without stones on the
board (see Figure 2; rotations not allowed), you may place a stone in each of
those cells.

Figure 2: The tromino region

• If all cells in a row or column have a stone, you may remove all stones from
that row or column.
For which $n$ is it possible that, after some nonzero number of moves, the
board has no stones?

This problem proved to be one of the most challenging ones in the history of
the USAMO: only seven students who took the exam were able to solve it. The
solution employs what has recently been called the polynomial method, which is
based on the simple fact that polynomials have a limited number of roots.
(Note that one can prove this fact with some basic algebra.) This technique
has seen a lot of interesting applications lately in both mathematics and
computer science, and we were excited to feature it in this beautiful problem.

It is not hard to see that when $n$ is divisible by $3$, the empty board is
achievable. Divide the board into 3-by-3 sub-boards. In each 3-by-3 sub-board,
follow the procedure shown in Figure 3.

Figure 3: Clearing the board when $n$ is divisible by 3

We can then use the polynomial method to prove that no procedure will succeed
when $n$ is not divisible by $3$. Leaving a bit of mystery, we only show this
here for $n=4$; readers are encouraged to think about how this generalizes to
other $n$ that is not divisible by $3$ (and why the argument fails when $n$ is
divisible by $3$).

We place monomials into the $16$ squares of our 4-by-4 board as shown in
Figure 4.

012301231$y$$y^{2}$$y^{3}$$x$$xy$$xy^{2}$$xy^{3}$$x^{2}$$x^{2}y$$x^{2}y^{2}$$x^{2}y^{3}$$x^{3}$$x^{3}y$$x^{3}y^{2}$$x^{3}y^{3}$

Figure 4: Labelling cells of a 4-by-4 board

We keep track of our moves as follows: when stones are placed on the board, we
add their values, and when stones are removed from the board, we subtract
their values. Assume now, indirectly, that a procedure results in an empty
board; to be specific, suppose that in this procedure, a tromino with its
lower left corner in position $(i,j)$ was added $t_{i,j}$ times; the $i$th
column was cleared $c_{i}$ times; and the $j$th row was cleared $r_{j}$ times.
This then means that we must have

$T(x,y)-C(x,y)-R(x,y)=0,$

where

$T(x,y)=\sum_{i=0}^{2}\sum_{j=0}^{2}t_{i,j}x^{i}y^{j}(1+x+y),$
$C(x,y)=\sum_{i=0}^{3}x^{i}(1+y+y^{2}+y^{3}),$
$R(x,y)=\sum_{j=0}^{3}y^{j}(1+x+x^{2}+x^{3}).$

Consider now the set $D=\\{-1,\mathrm{i},-\mathrm{i}\\}$, and observe that
$C(x,y)$ and $R(x,y)$ both equal zero for every $x\in D$ and $y\in D$, but
$1+x+y$ is never equal to zero for any such $x$ and $y$. We then must have

$P(x,y)=\sum_{i=0}^{2}\sum_{j=0}^{2}t_{i,j}x^{i}y^{j}=0$

for each $x\in D$ and $y\in D$. This is a contradiction, however, since the
polynomial $P(x,y)$ is at most quadratic in both $x$ and $y$, and thus it
cannot have more than two values of each variable for which it equals zero.

The AMC program would not be possible without our four remarkable editorial
boards that create our exams each year: they propose problems, review and rate
all submissions, select the problems for the exams, and meticulously edit
these problems and their solutions. They also are instrumental in helping us
shape general AMC policies and practices. These nearly 180 individuals come
from all across the nation and from abroad; include people from academia, high
schools, business, and industry; have a wide range of ages represented, from
undergraduates to retired mathematicians; and possess an impressive variety of
mathematical and cultural expertise. I am immensely grateful to this talented,
diverse, and accomplished group of people. Anyone with an interest in joining
us is encouraged to contact me.

Credits

Figure 1 is courtesy of Ricardo Conceição.

Figure 2 is courtesy of Shaunak Kishore and Alex Zhai.

Figures 3 and 4 are courtesy of the author.

Photo of Béla Bajnok is courtsey of Steve Berg.

DRAFT September 2, 2023