1.5exL1.2

# A Barrier Method for Contact Avoiding Particles in Stokes Flow

Anna Broms1∗ and Anna-Karin Tornberg1

1 Department of Mathematics
KTH Royal Institute of Technology
Lindstedtsvägen 25, 114 28 Stockholm, Sweden
∗e-mail<EMAIL_ADDRESS>https://www.kth.se/profile/annabrom

###### Abstract

Rigid particles in a Stokesian fluid can physically not overlap, as a thin
layer of fluid always separates a particle pair, exerting increasingly strong
repulsive forces on the bodies for decreasing separations. Numerically,
resolving these lubrication forces comes at an intractably large cost even for
moderate system sizes. Hence, it can typically not be guaranteed that particle
collisions and overlaps do not occur in a dynamic simulation, independently of
the choice of method to solve the Stokes equations. In this work, non-overlap
constraints, in terms of the Euclidean distance between boundary points on the
particles, are represented via a barrier energy. We solve for the minimum
magnitudes of repelling contact forces between any particle pair in contact to
correct for overlaps by enforcing a zero barrier energy at the next time
level, given a contact-free configuration at a previous instance in time. The
method is tested using a multiblob method to solve the mobility problem in
Stokes flow applied to suspensions of spheres, rods and boomerang shaped
particles. Collision free configurations are obtained at all instances in
time. The effect of the contact forces on the collective order of a set of
rods in a background flow that naturally promote particle interactions is also
illustrated.

Key words:  Stokes flow, contact problem, rigid particles, barrier method
Highlights

  * •

Strategy presented to avoid numerically introduced particle contacts in Stokes
flow.

  * •

Non-colliding bodies guaranteed by enforcing zero barrier energy at each time-
step.

  * •

Contact force magnitudes are minimised for reduced effect on the physics of
the system.

  * •

Contact forces and contact distances satisfy discrete complementarity.

  * •

Both convex and non-convex rigid particles can be handled robustly.

##  1 Introduction

We present an algorithm for contact forces between rigid particles immersed in
an unbounded viscous fluid in 3D, introduced only when needed to avoid
unphysical particle collisions and overlaps. Such collisions are caused from
hard-to-avoid numerical artifacts and can also be the results of Brownian
increments in a stochastic setting where thermal fluctuations in the fluid are
considered. In this work, full hydrodynamic interaction is encountered for,
meaning that there is a global coupling between all the particles in the
system.

We compute these hydrodynamic interactions with a cheap and fast technique:
the so called rigid multiblob method, as carefully described in [1, 2]. Note
that this is only one of many numerical methods for Stokes flows and the
contact avoiding strategy that we develop can be used also in combination with
other techniques. The idea in the multiblob method is to model a rigid body by
a collection of spheres or “blobs” and has been used in a very large number of
works, see e.g. [3, 4, 5, 6, 7] and references therein. Each blob interacts
hydrodynamically with all other blobs in the system in a pairwise manner and
blobs belonging to the same particle are constrained to move as a rigid body
via forces applied at each blob center, such that the blob forces sum to the
net force and torque on the particle. Mathematically, one could view this as a
regularised single layer boundary integral formulation, with the blob radius
the regularisation parameter. Example multiblob geometries are displayed in
Figure 1. For each instance in time, we solve the Stokes mobility problem,
that is, given assigned external forces and torques on the particles, such as
e.g. gravity or some electrostatic forcing, the resulting particle
translational and rotational velocities are computed. The inertia of the
particles is typically negligible in the Stokes regime, and hence, the
computed rigid body velocities can be used to update the particle positions,
applying some suitable time-stepping scheme, and the configuration of
particles can in this way be studied dynamically. By choosing the
regularisation parameter and the surface where blobs are placed in relation to
the true surface of the particle as the solution to a small off-line
optimisation problem for each (axisymmetric) particle type, good accuracy in
the particle velocities can be obtained for moderately separated particles
even with coarse grids of the particle surfaces [2].

(a) A suspension of slender rods.

(b) A multiblob boomerang with its center of mass indicated with a cross.

Figure 1: Multiblob particles in a Stokesian fluid. Colors indicate the depth
in the figures.

As particles get closer however, their interactions become increasingly
difficult to accurately resolve. In the multiblob framework, accuracy is
suffering for particles close to being in contact, and this is a challenge
independently of the choice of Stokes solver. Insufficient treatment of close
interactions, caused either by numerical errors due to insufficient spatial
resolution of the particle grid or by accumulated errors from the time-
stepping scheme of the moving particles, can in the worst case lead to non-
physical particle overlaps. Even if close interactions are resolved, and an
adaptive time-stepping scheme is used, another unwanted effect is stalling
when the adaptive time-step size becomes increasingly small as the time-step
is adjusted to the stiff problem of particles in close proximity with
increasingly strong repulsive forces. Lubrication forces between particles,
which physically guarantees no-collision, can be resolved only using high
fidelity methods with high spatial discretisation of the surfaces of the
particles [8]. Both high temporal and spatial discretisation hence comes at a
large cost, especially in dense suspensions. Keeping the surface grid
resolution moderate, simulations become cheaper, but the accuracy is worsened
and there is a large need for a contact avoiding strategy.

Repelling contact forces should be introduced in such a way that the physical
properties of the system are preserved. Contact forces are artificial to the
system and introduced only to decrease the impact from the equally artificial
numerical errors difficult to avoid for particles in close proximity. Hence,
we want to find the smallest possible contact forces for particles to stay
apart. At the same time, we also want to avoid introducing stiffness in the
system, and eliminate the risk of having to take very small time-steps. Such
stiffness is often the draw back of applying a strongly repelling potential to
avoid particle overlaps, as an alternative to contact resolution algorithms.
Potentials of this type include Lennard-Jones or e.g. a potential based on
Hard Gaussian Overlap (HGO), which is designed for ellipsoids [9] and
discussed for rods in [10]111For multiblob particles of a more general
geometry, one could construct a potential from an HGO-potential centered on
each individual blob building up the particle.. Due to the problem with
stiffness, it is difficult to guarantee non-overlapping particles with a
potential-based method, and there is a risk that particles become “soft” [11].

Contact problems with particles _not_ immersed in a fluid has a richer
literature [12, 13, 14, 15, 16] with more benchmarks available. An overview of
techniques for contact dynamic problems is found in [17]. There is also a body
of work in computer graphics [18, 19], and of most relevance to this work, the
recently introduced method of Incremental Potential Contact (IPC), which is a
penalty method where colliding and overlapping configurations are penalised by
a barrier energy [20, 21, 22].

Different contact resolution techniques have been suggested to better resolve
particle collisions in Stokes flow, see e.g. the works [23, 24, 25, 26, 27],
where no-slip boundary conditions are imposed at the point of contact,
constraining the colliding particles to move with equal speed at the collision
point. A complementarity formulation is favored in [28, 29, 30, 31, 32, 33],
where the idea is to solve a nonlinear complementarity problem for contact
force magnitudes and some definition of separation. The difficult-to-solve
nonlinear problem is in turn approximated by one or a sequence of _linear_
complementarity problems (LCPs) that can be solved more easily. One method in
this class is presented in the works by Yan et al., [29, 30, 32], and has the
advantages that it is easily applicable to any hydrodynamic solver and is
based on a geometric formulation that fulfills Newton’s third law (forces on
every pair of particles are balanced). One drawback, however, is that the
method cannot guarantee non-overlapping configurations at the end of each
time-step, as the solution of a single LCP does not imply a solution to the
original nonlinear problem. A second drawback is that non-convex particles
cannot be handled. Another method of the same flavour for Stokesian fluids is
based on so called Space-Time Interference Volumes (STIV) by Lu et al. [28,
34] and Bystricky et al. [31]. The STIV technique is inspired by the work of
Harmon et al. [19] in computer graphics. The general idea is to consider
Stokes equation in variational form and after a candidate time-step, compute
the volume in space-time swept out by the trajectories of the particles in
contact. If the volume is negative, particles overlap. Enforcing this volume
to be zero, by determining appropriate repulsion forces, gives rise to a
complementarity problem. In an STIV approach, no-collision is obtained even
for large time-steps by solving sequences of LCPs and particles that pass
through each other during one time-step can be detected and the corresponding
time-step corrected. However, each STIV-formulation is strongly linked to a
specific Stokes solver, the method is not easily applicable to general
geometries, and in 3D, the STIV volume has to be efficiently computed in four
dimensions [34, 28]. Another difficulty is that using the STIV, contact forces
are not automatically balanced and balanced forces are difficult to obtain
(Newton’s third law is _not_ satisfied) [35].

We choose to approach the contact problem by introducing repelling contact
forces to particles that come too close to each other in terms of the
Euclidean distance. Key features of the method is that all contact forces are
balanced via Newton’s third law and that the next time-step in a sequence of
time-steps is guaranteed to be contact free (in fact, by construction, a
minimum separation distance is guaranteed between particles). The method is
strongly inspired by the work of Yan et al. [29, 30, 32] in how contact forces
are geometrically motivated and by the work of Zorin and coauthors in the work
on IPC, [20, 22, 21], in how non-overlapping constraints are set up. However,
instead of penalising contact by a barrier energy as in [20, 22, 21], a
barrier energy is rather used to represent a large number of non-overlap
constraints and the complementarity condition with the associated contact
force. In the work on IPC, [20, 22, 21], a tetrahedral or triangular mesh of
the surface of each particle is considered, where distances are computed
robustly between all pairs of edges and points to triangles respectively. In
our work, particles are rigid and smooth and we can utilise the known
parameterisation of the particle surfaces to compute distances robustly; given
a set of points defining the grid of a particle surface, we instead flag the
closest point of contact on the neighbouring particle to be part in the
collision handling for each surface grid node within some set threshold of the
other particle. This is possible as we are not dependent on second order
derivatives of the distances in solving the optimisation problem with our
formulation, in contrast to IPC; gradients of the distances with respect to
particle coordinates suffice. The handling of the geometry and the constraints
is also in contrast to the geometric approach by Yan et al, where only a
single point of contact or ”the most overlapping point” has to be determined,
and more similar to the STIV technique, where multiple segments of the
boundary can be flagged for the same particle contact pair. This is a
beneficial property especially for non-convex particle geometries, or if
surfaces are close to parallel, where the computed contact torque otherwise
becomes very sensitive to the choice of contact point. Details are outlined in
Section 2.1. Another difference in our work compared to [20, 22, 21] is that
particles are immersed in a Stokesian fluid, where intertia is negligible, and
we solve for the force magnitudes rather than the particle coordinates in an
implicit time-step as in IPC for rigid bodies [21]. The dimension of the
optimisation variable in the resulting optimisation problem hence becomes
smaller. Keeping the rigidity of the particles is a non-issue (no additional
constraints have to be enforced as discussed in [21]). We utilise the
linearity of the Stokes equations and formulate a minimisation problem for
contact force magnitudes to solve in every time-step where contact occurs.
Note that we do not minimise a combined, weighted, energy formulation as in
IPC [20, 22, 21]. Hence, there are no parameters to tune in our formulation
for the importance of the collision constraints relative to the minimisation
of other contributions to the total energy of the system.

The time-stepping schemes used both for STIV and the geometric approach by Yan
et al. is explicit (to avoid a large cost at every time-step) and of low order
[36, 29, 30, 32]. The repulsion force is assumed to be constant over the
course of one time-step and the basis for both types of contact algorithms in
their vanilla version is explicit Euler. This is the time-stepping method that
will be used for demonstration also in this work. An additional motivation to
this choice is in settings where the contact resolution algorithm is coupled
to Brownian motion, modelled by a stochastic differential equation, where it
is difficult to obtain anything better than first order accuracy in time. Note
however that there is nothing that prevents a higher order time-stepping
scheme from being used if contact avoiding is applied in a deterministic
setting. The most straight-forward idea is then to use a higher order time-
stepping method for the particles not involved in any contacts. For particles
to which repulsive forces are added in a certain time-step, there is no
expected gain in using a high order method.

###  1.1 The Stokes mobility problem

Before introducing the optimisation problem, we start with some preliminaries.
Each 3D particle in the fluid suspension can be described by its center
coordinates $\boldsymbol{x}_{i}$ and rotation quaternion,
$\boldsymbol{q}_{i}$. We collect these generalized coordinates for all the $N$
particles in the system in the vector $\boldsymbol{\mathcal{Q}}$ such that

$\boldsymbol{\mathcal{Q}}=\left[\boldsymbol{x}_{1}^{T},\boldsymbol{q}_{1}^{T},\boldsymbol{x}_{2}^{T},\boldsymbol{q}_{2}^{T},\dots,\boldsymbol{x}_{N}^{T},\boldsymbol{q}_{N}^{T}\right]^{T}.$
(1)

Let $\boldsymbol{\mathcal{U}}$ be a vector of all rigid body velocities of the
particles in the system, where $\boldsymbol{u}_{i}\in\mathbb{R}^{3}$ is the
translational velocity and $\boldsymbol{\omega}_{i}\in\mathbb{R}^{3}$ the
rotational velocity of particle $i$. Similarly, let
$\boldsymbol{\mathcal{F}}_{\text{ext}}$ be a vector of all the externally
applied forces $\boldsymbol{f}_{i}\in\mathbb{R}^{3}$ and torques
$\boldsymbol{t}_{i}\in\mathbb{R}^{3}$ on the particles in the system:

$\boldsymbol{\mathcal{U}}=\begin{bmatrix}\boldsymbol{u}_{1}^{T}&\boldsymbol{\omega}_{1}^{T}&\boldsymbol{u}_{2}^{T}&\boldsymbol{\omega}_{2}^{T}&\dots&\boldsymbol{u}_{N}^{T}&\boldsymbol{\omega}_{N}^{T}\end{bmatrix}^{T},\quad\boldsymbol{\mathcal{F}}_{\text{ext}}=\begin{bmatrix}\boldsymbol{f}_{1}^{T}&\boldsymbol{t}_{1}^{T}&\boldsymbol{f}_{2}^{T}&\boldsymbol{t}_{2}^{T}&\dots&\boldsymbol{f}_{N}^{T}&\boldsymbol{t}_{N}^{T}\end{bmatrix}^{T}.$
(2)

Now, the Stokes mobility problem can be stated as

$\boldsymbol{\mathcal{U}}=\boldsymbol{\mathcal{U}}_{\text{bg}}+\boldsymbol{\mathcal{U}}_{\text{Brownian}}+\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{F}}_{\text{ext}},$
(3)

where $\boldsymbol{\mathcal{M}}$ is the mobility matrix of size $6N\times 6N$,
with each $6\times 6$ block corresponding to the interaction between a
specific pair of particles. The vector $\boldsymbol{\mathcal{U}}_{\text{bg}}$
is the velocity contribution on the particles from an eventual background flow
and $\boldsymbol{\mathcal{U}}_{\text{Brownian}}$ is a vector of stochastic
velocities included if thermal fluctuations are considered. Note that
implementation-wise, $\boldsymbol{\mathcal{M}}$ is only to be interpreted
symbolically and
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{F}}_{\text{ext}}$ corresponds to
solving Stokes equations along with no-slip boundary conditions using the
multiblob method as described in [2], given forces and torques for all
particles, stacked in $\boldsymbol{\mathcal{F}}_{\text{ext}}$.

With no contact forces present, $\boldsymbol{\mathcal{Q}}$ is related to the
velocities in $\boldsymbol{\mathcal{U}}$ as

$\dot{\boldsymbol{\mathcal{Q}}}=\boldsymbol{\Psi}\boldsymbol{\mathcal{U}},$
(4)

with $\boldsymbol{\Psi}$ a geometry-dependent matrix relating velocities to
particle positions and quaternions.

When particles are in contact, the velocities are corrected using the action
of a vector of contact forces $\boldsymbol{\mathcal{F}}_{c}$, such that

$\dot{\boldsymbol{\mathcal{Q}}}=\boldsymbol{\Psi}\boldsymbol{\mathcal{U}}+\boldsymbol{\Psi}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{F}}_{c},$
(5)

with

$\boldsymbol{\mathcal{F}}_{c}=\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}.$
(6)

The matrix
$\boldsymbol{\mathcal{D}}\in\mathbb{R}^{6N}\times\mathbb{R}^{N_{c}}$
determines the direction of contact forces and torques and depends on the
particle configuration and geometries and the vector
$\boldsymbol{\lambda}\in\mathbb{R}^{N_{c}}$ determines the contact force
magnitudes, with $N_{c}$ the number of particle pairs in contact. Ideally, we
would like to determine contact forces such that if contact forces are needed
for particles to stay apart during the time-step, the force should be active
and the corresponding force magnitude positive. On the other hand, if the
particles stay separate without contact forces, the contact force magnitude
should be zero and the contact force passive. The nature of the contact
problem is hence a complementarity problem. It is also a nonlinear problem, as
the distances at the next instance of time depends nonlinearly of the contact
forces.

##  2 A barrier method for particles in contact

In practice, it is not advisable to let particles come so close to each other
that they eventually touch, due to limitations in all numerical solvers for
particles in Stokes flow – the accuracy for almost touching particles is
suffering if the resolution of the particles is not excessively high. For this
reason, we define a contact to occur if the distance $d$ between points on
particles in close proximity is smaller than a set buffer distance $\hat{d}$.
The parameter $\hat{d}$ is typically set to ensure a prescribed accepted
accuracy from the multiblob method222One could also consider to set $\hat{d}$
in relation to any inhomogenities on the surfaces of the physical particles
modelled by our rigid particles.. We would like to find contact force
magnitudes $\boldsymbol{\lambda}$ such that this separation distance is
guaranteed at the next time-step, i.e.

$d_{i}\left(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}(\boldsymbol{\lambda})\right)>\hat{d},\quad\text{for all $i$ pairs of points
on distinct particles},$ (7)

with $d_{i}$ depending non-linearly on $\boldsymbol{\mathcal{Q}}_{t+\Delta
t}$, the particle configuration in the next time step, and the set of relevant
$d_{i}$ for each contact pair properly defined in Section 2.1. The
complementarity condition for the particle distances and contact force
magnitudes take the form

$\lambda_{k}\max\left(0,\min\limits_{i}\left\\{d_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t})-\hat{d}\right\\}\right)=0,\quad\text{$i$ associated with particle contact
pair $k$ and }\lambda_{k}\geq 0.$ (8)

Instead of solving this (modified) nonlinear complementarity problem sharply,
we will set up a contact optimisation problem to fulfill all constraints in
(7) and still approximately solve (8). Note that depending on the number and
concentration of particles in the system and their geometries, the number of
constraints in (7) might be very large. This section explains how such a
minimisation problem can be formulated and discuss modelling choices for the
objective function, the constraints and the geometry of the contact forces.
The constraints in (7) can be represented via a minimisation of a sum of
indicator functions [37]:

$\min\limits_{\boldsymbol{\lambda}\geq\boldsymbol{0}}\quad\sum_{i}I(d_{i}\left(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}(\boldsymbol{\lambda})\right)-\hat{d}),\quad\text{with
}I(s)=\left\\{\begin{aligned} 0,\quad s\geq 0,\\\ \infty,\quad
s<0.\end{aligned}\right.$ (9)

This representation is however difficult to work with and we therefore replace
(9) with a smoother counterpart. Introduce a so-called barrier function, $b$,
that is zero for sufficiently large distances $d$ and increasingly large for
distances smaller than the set threshold $\hat{d}$:

$b(d,\hat{d})=\left\\{\begin{aligned}
-(d-\hat{d})^{2}\ln\left(d/\hat{d}\right),\quad&0<d<\hat{d},\\\
0,\quad&d\geq\hat{d}.\end{aligned}\right.\\\ $ (10)

We can collect all the non-overlapping constraints in (7) into what we define
as a _barrier energy_ , mimicking the sum of indicators in (9). A non-
overlapping configuration is one with non-negative contact force magnitudes
$\boldsymbol{\lambda}$ applied for $t\in[t,t+\Delta t]$ that minimises the
barrier energy

$\min\limits_{\boldsymbol{\lambda}\geq\boldsymbol{0}}\quad\sum_{i}b\left(d_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}),\hat{d}\right),$ (11)

with the summation being over all geometries in contact (corresponding to all
constraints in (7)). In the Stokesian fluid, the particle coordinates evolve
with time according to

$\dot{\boldsymbol{\mathcal{Q}}}=\boldsymbol{\Psi}\left(\boldsymbol{\mathcal{U}}+\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}\right).$
(12)

If we discretise this equation in time, we can express
$\boldsymbol{\mathcal{Q}}_{t+\Delta t}$ in terms of previous, known,
coordinate vectors. With the forward Euler method, the minimisation problem
becomes

$\displaystyle\min\limits_{\boldsymbol{\lambda}\geq\boldsymbol{0}}\quad$
$\displaystyle\sum_{i}b\left(d_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}),\hat{d}\right),$ (13) s.t.
$\displaystyle\boldsymbol{\mathcal{Q}}_{t+\Delta
t}=\boldsymbol{\mathcal{Q}}_{t}+\Delta
t\boldsymbol{\Psi}\left(\boldsymbol{\mathcal{U}}+\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}\right).$

This is a constrained minimisation problem to solve for
$\boldsymbol{\lambda}\in\mathbb{R}^{N_{c}}$, i.e. with a scalar contact force
magnitude per particle pair in contact. Note however that for some particle
configurations, any vector $\boldsymbol{\lambda}$ with sufficiently large
magnitude would solve (13). We hence need to constrain the solution by
penalising large $\boldsymbol{\lambda}$, with a natural choice being

$\min\limits_{\boldsymbol{\lambda}\geq\boldsymbol{0}}\quad\|\boldsymbol{\lambda}\|+\alpha\sum\limits_{i}b\left(d_{i}\left(\boldsymbol{\mathcal{Q}}_{t}+\Delta
t\boldsymbol{\Psi}\left(\boldsymbol{\mathcal{U}}+\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}\right)\right),\hat{d}\right),$
(14)

where $\alpha$ is a parameter that we have to set to balance the minimisation
of the barrier energy at the next time-step and the minimisation of contact
forces. The parameter $\alpha$ also has the role of enforcing the
complementarity condition for the contact forces and contact distances
(relative to the buffer region) in (8). A different way of saying the same
thing is that we would like to pick $\lambda_{k}$ sufficiently small so that
non-overlap is obtained, but where the force is active, the minimum contact
distance should not be larger than $\hat{d}$ after applying the contact
forces.

For the choice of norm in (14), we have two options: A 1-norm penalty of the
contact force magnitudes, $\sum_{i}\lambda_{i}$, penalises the forces at all
contacts to an equal extent and not only those where the magnitude is large.
One could also choose the norm
$\left(\boldsymbol{\lambda}^{T}\boldsymbol{\mathcal{D}}^{T}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}\right)^{1/2}$,
or its square, corresponding to the dissipated energy induced by the contact
forces. As a main goal is to minimise the impact of the artificial repulsive
forces on the system, the dissipated energy due to the introduced contact
forces could be a reasonable choice of objective function. For our Stokes
solver, the rigid multiblob method, it is a known problem that the accuracy in
$\boldsymbol{\mathcal{M}}$ is suffering a lot for closely interacting
particles [2]. Hence,
$\boldsymbol{\lambda}\boldsymbol{\mathcal{D}}^{T}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$
is a very crude approximation to the true dissipated energy and a large error
is associated with this quantity. For this reason we will use
$\sum_{i}\lambda_{i}$ as our objective function of choice. See section 3.1.2
for a numerical comparison of the two choices.

Note that it is not trivial to find the best choice of $\alpha$ in (14) for
general applicability as the magnitude of both $\boldsymbol{\lambda}$ and the
barrier energy depends on the number of pairs of particles in contact and how
much overlap a non-corrected time-step would yield. A suitable level of
$\alpha$ also depends on the tolerance chosen as stopping criteria when
solving the problem in (14). This is a similar problem as the one encountered
when minimising the energy formulation in the work on IPC, where a penalty has
to be chosen in an adaptive manner to give the barrier energy the proper
weight, see the supplementary in [20]. Here, we would like to avoid such
hyperparameter tuning. What _is_ known, however, is that the barrier energy is
identically zero for a non-overlapping configuration. An alternative
formulation of the minimisation problem in (14) is therefore

$\displaystyle\min_{\boldsymbol{\lambda}\geq\boldsymbol{0}}\quad\sum_{i}\lambda_{i},$
(15) $\displaystyle\text{s.t.
}\sum\limits_{i}b\left(d_{i}\left(\boldsymbol{\mathcal{Q}}_{t}+\Delta
t\boldsymbol{\Psi}\left(\boldsymbol{\mathcal{U}}+\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}\right)\right),\hat{d}\right)=0,$

where contact force magnitudes are minimised while enforcing a zero barrier
energy. The formulation in (15) is simpler than (14) in the sense that one
hyperparameter is reduced, but possibly more challenging as a nonlinear
equality constraint has been added. Note that the Lagrangians of (14) and (15)
are identical, but in the case of (15), the parameter $\alpha$ is a Lagrange
multiplier to be solved for instead of a set parameter. In the remaining of
this work, we will consider the formulation in (15) for determining contact
forces.

######  Remark 1.

In contrast to in the work on IPC by Zorin et al., the value of the barrier
energy itself is not used explicitly in the formulation presented in this
paper but necessary for

  1. 1.

Imposing a large number of constraints of minimum distances, as we seek only
the force magnitudes where the barrier energy is zero. What is gained is that
we do not have to deal with the large number of constraints explicitly and are
hence able to solve smaller optimisation problems in every time-step where
contact occurs.

  2. 2.

Obtaining a direction of the contact force and torque that takes more
information into account: by using the barrier energy, a large number of pairs
of contact points are weighted by their distance to give a net force and
torque (more on this in Section 2.1).

An alternative formulation for the same set of constraints is
$\min_{i}d_{i}\geq\hat{d}$. For a large class of optimisation methods, the
interior point methods, such a problem requires a feasible starting guess,
$\boldsymbol{\lambda}\geq\boldsymbol{0}$ such that all non-overlap constraints
hold, which might be very hard to obtain even for systems containing a few
particles. Rewriting the non-overlap constraints in (7) as a barrier energy is
hence necessary for robustness.

###  2.1 Defining the contact distance and contact force

Let us now focus on the Euclidean contact distances $d_{i}$ between points on
the surfaces of two particles in close proximity. We make a distinction
between $d^{\text{p}}$ and $d^{\text{s}}$, with $d^{\text{p}}$ the shortest
distance between particles and $d^{\text{s}}$ the shortest distance between a
pair of surface points on two particles. In principle, $d_{i}$ in the
constrained optimisation problem (15) can denote either $d_{i}^{\text{p}}$ or
$d_{i}^{\text{s}}$. Practically, we compute $d^{\text{p}}$ directly using
known parameterisations of the particles, while the surface-point-to-surface-
point distance $d^{\text{s}}$ is based on discretisations of the particle
surfaces, described by a distribution of points: For each node on the surface
of one particle in the pair sufficiently close to contact, we then use the
known parameterisation of the center line of the other particle to determine
the closest point on this line. The particle surfaces considered in this work
are all one radius away from the particle center line and the surface-point-
to-surface-point distance can hence easily be computed. Geometric
considerations determine if we choose to use $d^{\text{p}}$ or $d^{\text{s}}$:

  1. 1.

The particle-particle distance $d^{\text{p}}$ may be a good choice if the
particle shapes are simple enough, such as e.g. for spheres or axisymmetric
rods with semi-spherical caps. In the latter case, we solve an equation to
determine the closest distance between two line segments, following [38, 39],
and subtract $2R_{\text{rod}}$ to determine the closest distance between
particles, see Figure 2(c) for an illustration.

  2. 2.

The surface-point-to-surface-point distance, $d^{\text{s}}$, is especially
beneficial in two different cases: a.) For close to parallel surfaces, such
that it is hard to define a single point of contact and b.) For non-convex
particles. We will focus most of our attention on this case.

(a) Contact distances $d_{ik}^{\text{s}}$ are the distances from each surface
grid node on one particle (black dots) to the closest point on the other
particle, less than a set threshold $\hat{d}_{\text{try}}$ (buffer region
indicated in yellow). Closest points are determined from the shortest distance
to the center line segment of the other particle (blue line). Pairs of grid
nodes and computed points are marked with red dots, with red lines drawn in
between displaying the distance.

(b) The corresponding directions of the contact force, representing the terms
in (23), are visualised with red arrows scaled correlated with their
contribution to the total force. There is one red arrow for each red distance
line $d_{ik}^{\text{s}}$ marked in (a). The resulting contact force direction
on each particle in (23) is indicated with black thick arrows.

(c) The shortest particle-particle distance, $d_{i}^{\text{p}}$, is given by
the length of the red line drawn between the closest points on the two
surfaces, computed from the closest distance between the two center line
segments, and subtracting $2R_{\text{rod}}$. The resulting contact force
direction on each particle is indicated with black thick arrows and have the
same direction as the normal drawn between the contact points.

Figure 2: Illustration of the two contact distances $d_{i}^{\text{p}}$ and
$d_{ik}^{\text{s}}$ and the corresponding contact force directions on the
particles for two fat rods close to contact in pair $i$.

We think of the contact forces as a (hopefully small) correction to the
_trial_ time-step $\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*}$, i.e. the time-
step without contact forces, given by

$\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*}=\boldsymbol{\mathcal{Q}}_{t}+\Delta
t\boldsymbol{\Psi}\boldsymbol{\mathcal{U}}.$ (16)

Let $B_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*},\hat{d}_{\text{try}})$ be
a modified barrier energy associated with the particle contact pair $i$
aggregated over all relevant particle-particle distances at the trial time-
step, with the threshold $\hat{d}_{\text{try}}$ chosen so that
$\hat{d}_{\text{try}}>\hat{d}$ ($B_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}_{\text{try}})$ is the sum in (11), but for a single particle
pair and with modified threshold). Then, we define the _contact force
potential_ ,
$\sum\limits_{i}\lambda_{i}B_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}_{\text{try}})$, a weighted sum of the barrier energy in the
trial time-step, and let the contact forces be given by the gradient of this
contact force potential, i.e.

$\boldsymbol{\mathcal{F}}_{c}=-\nabla_{\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}}\left(\sum\limits_{i}\lambda_{i}B_{i}\left(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}_{\text{try}}\right)\right).$ (17)

As stated in (5)-(6), we may write the contact forces on the form
$\boldsymbol{\mathcal{F}}_{c}=\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$.
Let $\boldsymbol{B}$ be a vector of all the modified barrier energies for the
$N_{c}$ particle pairs in contact. The sparse matrix
$\boldsymbol{\mathcal{D}}\in\mathbb{R}^{6N\times N_{c}}$, representing contact
force and torque directions, depends on the geometry of the particles and
their locations at the trial time-step. From (17) we identify that
$\boldsymbol{\mathcal{D}}$ is given by

$\boldsymbol{\mathcal{D}}=-\nabla_{\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}}\boldsymbol{B}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}_{\text{try}}).$ (18)

The matrix has the structure

$\boldsymbol{\mathcal{D}}=\begin{bmatrix}\boldsymbol{D}_{1}&\boldsymbol{D}_{2}&\dots&\boldsymbol{D}_{N_{c}}\end{bmatrix},$
(19)

where $\boldsymbol{D}_{i}$ represents contact $i$, between particle $j$ and
$k$. In general, $\boldsymbol{D}_{i}$ takes the form

$\boldsymbol{D}_{i}=\begin{bmatrix}0\leavevmode\nobreak\
\dots\leavevmode\nobreak\
0&\dfrac{{\boldsymbol{\hat{f}}^{c}_{i}}^{T}}{\|\boldsymbol{\hat{f}}^{c}_{i}\|}&\dfrac{{\boldsymbol{\hat{t}}^{c}_{ij}}^{T}}{\|\boldsymbol{\hat{f}}^{c}_{i}\|}&0\leavevmode\nobreak\
\dots\leavevmode\nobreak\
0&-\dfrac{{\boldsymbol{\hat{f}}^{c}_{i}}^{T}}{\|\boldsymbol{\hat{f}}^{c}_{i}\|}&\dfrac{{\boldsymbol{\hat{t}}^{c}_{ik}}^{T}}{\|\boldsymbol{\hat{f}}^{c}_{i}\|}&0\leavevmode\nobreak\
\dots\leavevmode\nobreak\ 0\end{bmatrix}^{T},$ (20)

see [15] for a detailed motivation. From (20) and the relation
$\boldsymbol{\mathcal{F}}_{c}=\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$,
it is clear that the forces applied to the two particles in pair $i$ are equal
in magnitude (with the magnitude given by $\lambda_{i}$), and differ only in
sign, due to Newton’s third law. In the contact force potential, all particles
or grid nodes closer to each other than a tolerance $\hat{d}_{\text{try}}$ are
flagged from the trial configuration that _might_ come into contact in the
corrected time step. All these flagged points are considered when setting up
$\boldsymbol{\mathcal{D}}$ and determines the direction of contact forces and
torques. The threshold $\hat{d}_{\text{try}}$ is chosen to set up a buffer,
not to miss any colliding particles and sets at the same time the dimension of
the optimisation problem, i.e. the number of colliding particle pairs $N_{c}$.
The contact force potential and the force definition in (17) allows for
contact force complementarity _with respect to_ $\hat{d}_{\text{try}}$ at the
_trial_ time-step, meaning that we allow for a non-zero contact force
component in a certain direction only if the associated contact distance at
the trial time-step is such that $d_{i}(\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta
t})<\hat{d}_{\text{try}}$.

######  Remark 2.

We could also choose to construct the contact force potential at the previous,
contact free time-step, with $\boldsymbol{\mathcal{Q}}_{t}$. Note that this
would require a larger $\hat{d}_{\text{try}}$, as particles are expected to
move more in a full time-step than with only the correction from contact
forces.

######  Remark 3.

Even if we here write the forces and torques as the gradient of a barrier
energy, it is not the gradient of a conservative potential as also
$\boldsymbol{\lambda}$ depends on $\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*}$,
implicitly.

Next, we will compare two different strategies of determining the contact
force potential and hence assembling the matrix $\boldsymbol{\mathcal{D}}$,
using only the very closest two points on two distinct particles or all pairs
of grid nodes sufficiently close to each other. These strategies correspond to
using $d^{\text{p}}(\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta t})$ or
$d^{\text{s}}(\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta t})$ in the contact force
potential. The differences are also outlined in Figure 2.

In the case of a contact force direction determined for the _particle pair_ in
contact $i$, corresponding to using $d_{i}^{\text{p}}$, we identify that
$B_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}_{\text{try}})=b(d_{i}^{p}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}),\hat{d}_{\text{try}})$. Let $\boldsymbol{n}_{i}$ be the outward unit
normal from the particle in the contact pair with the lowest index (pointing
away from the contact) in the trial time-step. We let

$\displaystyle\boldsymbol{\hat{f}}^{c}_{i}$ $\displaystyle=-\frac{\partial
b(x,\hat{d}_{\text{try}})}{\partial
d_{i}}\Big{|}_{d_{i}^{\text{p}}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*})}\boldsymbol{n}_{i},$ (21) $\displaystyle\boldsymbol{\hat{t}}^{c}_{ij}$
$\displaystyle=-\frac{\partial b(x,\hat{d}_{\text{try}})}{\partial
d_{i}}\Big{|}_{d_{i}^{\text{p}}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*})}\left(\boldsymbol{n}_{i}\times\boldsymbol{s}_{ij}\right)$
$\displaystyle{\boldsymbol{\hat{t}}^{c}_{ik}}$ $\displaystyle=\frac{\partial
b(x,\hat{d}_{\text{try}})}{\partial
d_{i}}\Big{|}_{d_{i}^{\text{p}}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*})}(\boldsymbol{n}_{i}\times\boldsymbol{s}_{ik}),$

where $\boldsymbol{s}_{ij}$ and $\boldsymbol{s}_{ik}$ are the vectors from the
center of each particle to the point of contact, in the global reference
frame.

We can also choose to build $\boldsymbol{\mathcal{D}}$ for all _pairs of
surface points_ in contact, with $d^{\text{s}}$ the definition of separation
used in the contact force potential. Let $l_{j}$ be the index of a grid node
on particle $j$, $l_{k}$ an index of a grid node on particle $k$ and
$\boldsymbol{r}_{l}$ denote a grid node. Then, let $\mathcal{C}_{i}$ denote
the set of grid node pairs on different particles that are within a distance
$\hat{d}_{\text{try}}$ from each other. Here, the modified barrier energy
$B_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*},\hat{d}_{\text{try}})$
denotes the barrier energy for all grid nodes involved in the collision for
the contact pair $i$,

$B_{i}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*})=\sum_{l\in\mathcal{C}_{i}}b\left(d^{\text{s}}_{l}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}),\hat{d}_{\text{try}}\right).$ (22)

Following (18), this means that the direction of the force will be given by
the sum of the normal directions for each pair of grid nodes close to contact,
weighted by the corresponding derivative of the barrier function. All in all,
the components of $\boldsymbol{D}_{i}$ take the form

$\displaystyle\boldsymbol{\hat{f}}^{c}_{i}$
$\displaystyle=\sum\limits_{(l_{j},l_{k})\in\mathcal{C}_{i}}\left(\frac{\partial
b(x,\hat{d}_{\text{try}})}{\partial
x}\Big{|}_{x=\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|}\right)\frac{\boldsymbol{r}_{l_{k}}-\boldsymbol{r}_{l_{j}}}{\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|},$
(23) $\displaystyle\boldsymbol{\hat{t}}^{c}_{ij}$
$\displaystyle=\sum\limits_{(l_{j},l_{k})\in\mathcal{C}_{i}}\left(\frac{\partial
b(x,\hat{d}_{\text{try}})}{\partial
x}\Big{|}_{x=\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|}\right)\left(\frac{\left(\boldsymbol{r}_{l_{k}}-\boldsymbol{r}_{l_{j}}\right)}{\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|}\times\left(\boldsymbol{r}_{l_{j}}-\boldsymbol{x}_{j}\right)\right),$
$\displaystyle\boldsymbol{\hat{t}}^{c}_{ik}$
$\displaystyle=\sum\limits_{(l_{j},l_{k})\in\mathcal{C}_{i}}\left(\frac{\partial
b(x,\hat{d}_{\text{try}})}{\partial
x}\Big{|}_{x=\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|}\right)\left(\frac{\left(\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\right)}{\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|}\times\left(\boldsymbol{r}_{l_{k}}-\boldsymbol{x}_{k}\right)\right).$

meaning that a certain force direction
$\left(\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\right)/\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|$
has a larger weight if the corresponding distance
$\|\boldsymbol{r}_{l_{j}}-\boldsymbol{r}_{l_{k}}\|$ is smaller and the overlap
relative to the minimum allowed distance $\hat{d}$ larger. Distances very
close to $\hat{d}$ give only very small contributions to the sums representing
force and torque directions in (23).

Note that alternatives could be considered for the barrier function, defining
it by

$\quad b_{0}(d,\hat{d})=\left\\{\begin{aligned}
-\ln\left(d/\hat{d}\right),\quad&0<d<\hat{d},\\\
0,\quad&d\geq\hat{d},\end{aligned}\right.\quad\text{or}\quad
b_{1}(d,\hat{d})=\left\\{\begin{aligned}
(d-\hat{d})\ln\left(d/\hat{d}\right),\quad&0<d<\hat{d},\\\
0,\quad&d\geq\hat{d},\end{aligned}\right.\\\ $ (24)

with regularity $C_{0}$ and $C_{1}$ respectively and discussed also in [20].
We choose the barrier function $b$ in (10) due to its smooth transition at
$d=\hat{d}$ ($b$ has regularity $C_{2}$) which makes $b$ well-suited for the
optimisation problem in (15), especially for computing derivatives of the
barrier function which has to be done in any iterative method to solve (15).
For determining the contact force direction at the trial time-step, we have no
such regularity demand and choose the $C^{0}$ function $b^{0}$, for which the
weight in (23) from the derivative of the barrier function is the reciprocal
of the distance.

######  Remark 4.

Note that particles might not only violate the threshold $\hat{d}$ but also
overlap in the iterative optimisation procedure before a contact force of the
right magnitude is applied, such that $d_{i}<0$. We cannot evaluate the
barrier function $b$ in (10) with a negative argument. One solution is to map
$d_{i}\to d_{i}+\epsilon_{\text{reg}}$ and
$\hat{d}\to\hat{d}+\epsilon_{\text{reg}}$, with $\epsilon_{\text{reg}}$ a
parameter to be chosen. Note that if $\epsilon_{\text{reg}}$ is chosen large,
the gradient of the barrier function is small also for small $d_{i}$ relative
to $\hat{d}$, which might have an impact on the performance of solving the
optimisation problem (15). Note also that in a Brownian setting, we may need
to choose a large $\epsilon_{\text{reg}}$ if a large time-step size is used
(depending on how repelling the conservative potential is that gives rise to
external forces and torques). A different idea to avoid a large
$\epsilon_{\text{reg}}$ is to map all $d_{i}<\epsilon_{\text{cap}}$ linearly
such that the barrier function at $d_{i}=\epsilon_{\text{cap}}$ is $C^{1}$ and
extended to a linear function, with a typical choice of
$\epsilon_{\text{cap}}$ being a small fraction of $\hat{d}$, e.g.
$\epsilon_{\text{cap}}=10^{-2}\hat{d}$. Such an extension is visualised in
Figure 3(b). The three alternative barrier functions in (10) and (24) and
their regularisations are visualised in Figure 3. In numerical experiments in
Section 3, we employ the strategy with $\epsilon_{\text{reg}}$, as we
numerically have observed a reduced number of iterations for solving the
optimisation problem with this choice, as gradients of $b$ are smaller and
hence easier to handle.

(a) Modification of the barrier function $b$, where negative inter-particle
distances are allowed corresponding to a maximum overlap of size
$\epsilon_{\text{reg}}$.

(b) Alternative modification of the barrier function $b$ and its cousins
$b_{0}$ and $b_{1}$, where a particle overlap with $d<\epsilon_{\text{cap}}$
is mapped to a linear continuation with $C^{1}$ regularity at
$d=\epsilon_{\text{cap}}$.

Figure 3: Barrier functions as defined in (10) and (24) are smooth
approximations of the indicator function. Here, the barrier function is
modified to be defined for negative inter-particle distances, as particles
might overlap in the trial time-step and in iterations of the optimisation
method before contact forces of the right magnitudes are found.

###  2.2 Connection to a complementarity formulation

For comparison to complementarity techniques in the literature, as in [28, 29,
30, 31, 32, 33], the accumulated non-overlap constraints in
$\sum\limits_{i}B_{i}(\boldsymbol{\lambda})=0$ can equivalently be expressed
as a complementarity problem for the barrier energy of each contact pair and
the corresponding contact force magnitude:

$\boldsymbol{0}\leq-\boldsymbol{B}(\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta
t}+\Delta
t\boldsymbol{\Psi}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda},\hat{d})\perp\boldsymbol{\lambda}\geq\boldsymbol{0}.$
(25)

(despite the barrier function $b$ never taking negative values). The problem
in (25) can be solved by first linearising about the trial time-step so that

$\boldsymbol{0}\leq-\boldsymbol{B}(\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta
t},\hat{d})+\left(\nabla^{T}_{\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}}\boldsymbol{B}\left(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}\right)\Delta
t\boldsymbol{\Psi}\boldsymbol{\mathcal{M}}\nabla_{\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}}\boldsymbol{B}\left(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*},\hat{d}_{\text{try}}\right)\right)\boldsymbol{\lambda}\perp\boldsymbol{\lambda}\geq\boldsymbol{0}$
(26)

and techniques for LCPs may be used as in [40]. To guarantee a solution to the
original nonlinear problem in (25), a sequence of linear problems of the form
in (26) have to be solved, each with an update of the trial time-step,
$\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta t}$, containing already computed
position updates by an accumulated contact force.

The equation in (25) is reduced to the formulation by Yan et al. in [29, 30,
32] if only one point is involved in the contact per particle per contact pair
so that

$\boldsymbol{0}\leq\boldsymbol{d}^{p}(\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}(\boldsymbol{\lambda}))-\hat{d}\perp\boldsymbol{\lambda}\geq\boldsymbol{0}.$
(27)

The only difference is then that we view contact forces as a correction to a
trial time-step at $t+\Delta t$ and linearise about
$\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*}$ instead of linearising about the
previous coordinate vector $\boldsymbol{\mathcal{Q}}_{t}$, which is done in
[29, 30, 32]. At $\boldsymbol{\mathcal{Q}}_{t}$, the configuration is contact
free. That is however not guaranteed at $\boldsymbol{\mathcal{Q}}_{t+\Delta
t}^{*}$, which likely affect the speed of convergence of the optimisation
problem, as we are likely to start with an infeasible initial guess. In this
work, in contrast to the work by Yan et al, we find a solution to the
nonlinear complementarity problem in (25) by computing a solution vector
$\boldsymbol{\lambda}$ satisfying
$\sum\limits_{i}B_{i}(\boldsymbol{\lambda})=0$, instead of considering a
linearised problem. Note that it is not possible to strictly solve (25), as it
is possible to have a configuration with
$\boldsymbol{B}(\boldsymbol{\mathcal{Q}}^{*}_{t+\Delta t}+\Delta
t\boldsymbol{\Psi}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda},\hat{d})=\boldsymbol{0}$
and $\boldsymbol{\lambda}=\boldsymbol{0}$ and hence we allow for some
relaxation in the complementarity condition. We leave it to later work for a
more thorough comparison between the nonlinear barrier method presented in
this paper and versions of solving linear complementarity problems if the form
in (26) and as discussed in [28, 29, 30, 31, 32, 33].

######  Remark 5.

With the choice $\hat{d}_{\text{try}}=\hat{d}$, methods for Quadratic
Programming problems (QPs) are eligible to solve (26) (the matrix in the LCP
in (26) is symmetric positive semi-definite [29]) such as in [41, 42].
However, with such a small choice of the buffer region at the trial time-step
(implied by choosing $\hat{d}_{\text{try}}$ to be as small as the natural
choice of $\hat{d}$), there is a risk of missing potential contacts in the
corrected time-step, leading to a long sequence of linear complementarity
problems to be solved to eventually converge to a solution to the nonlinear
complementarity problem. On the other hand, if $\hat{d}$ is instead increased
to be as large as the natural choice of $\hat{d}_{\text{try}}$, to get
equality between the two parameters, we would have a very large buffer region
around each particle at the corrected time-step, with resulting contact forces
heavily affecting the physics of the suspension.

If the dissipative energy induced by contact forces is chosen as objective
function, the Lagrangian in our nonlinar optimisation problem (15) takes on a
similar form as in the works by Yan et al., in which a QP on the form

$\min\limits_{\boldsymbol{\lambda}\geq\boldsymbol{0}}\quad\boldsymbol{\lambda}^{T}\boldsymbol{\mathcal{D}}^{T}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}+\boldsymbol{G}^{T}\boldsymbol{\lambda},$
(28)

with $\boldsymbol{G}$ a specified vector is considered. For this QP, the first
order optimality conditions is an LCP for the contact distances at the known
time-step $\boldsymbol{\mathcal{Q}}_{t}$ and the force magnitudes
$\boldsymbol{\lambda}$. The difference in this paper is that we do not
consider a linearised complementarity problem and the second term in the
Lagrangian of (15) is nonlinear and corresponds to the nonlinear constraints.
Solving the problem becomes harder, but in contrast to the work by Yan et al.,
non-overlaps can be guaranteed.

###  2.3 Solution strategy

The problem (15) is solved with fmincon, an inbuilt solver for constrained
minimisation in Matlab, employing an interior-point method. In this and any
other iterative method that could be considered for solving the problem, a
full Stokes solve is required for each evaluation of the barrier energy,
requiring as much work as one time-step without contact handling. A large
number of such iterations hence become intractably expensive and it is
important to keep the number of iterations as low as possible. This number is
determined by a set of parameters that have to be carefully set. The time-step
size, $\Delta t$, is set a priori depending on the typical magnitudes of the
computed rigid body velocities in the problem, which depend on the type of
background flow and the magnitude of applied forces and torques, so that the
spatial updates of the particles per time-step are reasonable in magnitude. To
demonstrate the robustness of the method, we will vary the threshold
$\hat{d}$, determining the minimum allowed distance between particles and the
non-zero contribution to the barrier energy at the new time-step. In a general
application, we recommend to pick $\hat{d}$ as a fraction of the particle
radius, e.g. $\hat{d}=10^{-2}R$. As a rule of thumb for the remaining
parameters, we pick:

  * •

The stopping criterion for minimising contact force magnitudes,

$\max\limits_{i}|\lambda_{i}^{k+1}-\lambda_{i}^{k}|<10^{-2}\|\boldsymbol{\lambda}\|_{\infty}.$
(29)

  * •

The threshold $\hat{d}_{\text{try}}$ determining _both_ the $N_{c}$ particle
pairs to potentially be assigned contact forces _and_ the contact force and
torque directions for these pairs assembled in the matrix
$\boldsymbol{\mathcal{D}}$: For rods and boomerangs given $\Delta t$, we set
$\hat{d}_{\text{try}}$ adaptively in each time-step as

$\hat{d}_{\text{try}}=\hat{d}+\Delta
t\max_{i=1,\dots,N}\left(\boldsymbol{u}_{i}+(L/2)\boldsymbol{v}_{i}\times\boldsymbol{\omega}_{i}\right),$
(30)

where the expression in the parenthesis is the tip velocity of a particle,
with $\boldsymbol{v}_{i}$ the unit direction along the axis of particle $i$
and $L$ the particle length. For spheres, we pick

$\hat{d}_{\text{try}}=\hat{d}+\Delta
t\max\limits_{i=1,\dots,N}\boldsymbol{u}_{i}.$ (31)

  * •

The regularisation parameter $\epsilon_{\text{reg}}$ for the barrier function
as $\epsilon_{\text{reg}}=\hat{d}_{\text{try}}$.

Results are reported in Section 3. The influence of these parameter choices
will specifically be considered in Section 3.1.2.

##  3 Numerical results

We perform numerical experiments with spheres, rods and boomerangs. It is
difficult to study the performance of the contact algorithm dynamically, as a
small error in the time-discretisation or small differences in the contact
forces might result in completely different particle trajectories after
sufficiently long time. No standard benchmarks are available for systems of
multiple particles. For this purpose, we choose to focus on two different
types of tests:

  1. 1.

Investigating the ability of the contact resolution strategy to find contact
force magnitudes that maintain the allowed distance $\hat{d}$. This is a test
of the complementarity condition stating that a non-zero force should be
applied if and only if the contact is “active”. The test is performed for
geometries where the particles are pushed to come into contact by an external
force, starting from a contact free configuration.

  2. 2.

Investigating the ability to preserve symmetries in a background flow that
naturally enhance particle interactions to quantify the impact of different
choices of $\hat{d}$ and $\Delta t$.

In some examples, we will exaggerate $\Delta t$ and/or force magnitudes and
magnitudes of background flows to trigger the contact resolution algorithm
with the purpose of demonstrating its robustness.

###  3.1 Spheres

####  3.1.1 A random suspension of spheres

We perform an experiment with configurations of 500 non-overlapping unit
spheres randomly distributed in a cube of length $L$. The geometry is
exemplified in Figures 4(a) and 4(b) for the packing densities $24\%$ and
$12\%$. In each of 200 configurations for each density, $\hat{d}$ is set to be
the minimum separation distance. The spheres are then assigned external forces
uniformly sampled from a sphere and scaled so that
$\|\boldsymbol{f}_{i}\|=100$. A trial time-step with $\Delta t=0.01$ is taken
and for spheres that have come too close, contact forces are computed with
$d_{i}$ determined at the particle level ($d_{i}^{\text{p}}$ is used). Among
all the corrected particle pairs, the smallest distance at the next time-level
is reported vs $\hat{d}$ in Figure 4(d). In Figure 4(c) and 4(e), normalised
distances are displayed for the particles flagged to potentially be in contact
during the time-step, before and after applying the contact forces and scaled
relative to $\hat{d}_{\text{try}}$ and $\hat{d}$ respectively. In conclusion,
the forces do not artificially push the network of particles apart, but
approximately keep the inter-particle distances, except for the particle-pairs
that have come too close, for which the “overlap” is avoided by applying a
contact force.

(a) An example configuration with a packing density of $24\%$, with arrows
indicating randomly sampled external forces pushing the particles together.

(b) As in (a), but with a packing density of $12\%$.

(c) Inter-particle distances for the $N_{c}$ contact pairs before and after
applying contact forces, with statistics collected from the 200 configurations
with packing density $24\%$. For the size of $\hat{d}_{\text{try}}$ relative
to $\hat{d}$, see panel (d).

(d) Minimum distances upon applying contact forces,
$d_{\min}(\boldsymbol{\lambda})$, respect the set distance $\hat{d}$ both for
the denser and coarser configurations. For each configuration, the
corresponding $\hat{d}_{\text{try}}$ is also displayed, determining the
$N_{c}$ particle pairs flagged to be part of the contact force optimisation.

(e) Same as in (c) but with distances scaled relative to $\hat{d}$ instead of
$\hat{d}_{\text{try}}$.

Figure 4: Contact forces are computed for all particle pairs that violate
$\hat{d}$ among 500 randomly positioned spheres in a cube, where the packing
density is $24\%$ and $12\%$ respectively. The spheres are affected by
external forces with directions randomly drawn from the unit sphere and
$\|\boldsymbol{f}_{i}\|=100$ and a single time-step is taken with $\Delta
t=0.01$. At the trial time-step, particles are in contact if they are closer
to each other than $\hat{d}$, here picked for each random configuration as the
minimum separation distance at the previous, contact-free time-step. Minimum
separation distances $\hat{d}$ are respected with the contact forces, and
moreover, the forces do not drastically change the inter-particle distances
for the closest particles, meaning that the forces are not unnecessarily
large.

####  3.1.2 Hyperparameter robustness test

We perform the same test as for the dense suspension in Section 3.1.1 but vary
two parameters: $\hat{d}_{\text{try}}$, that determines the buffer region
around each particle in the trial time-step and the number of particle pairs
flagged to potentially come in contact during the time-step, and the stopping
criterion in the contact force optimisation, TOL, where the iterative
optimisation method is stopped if
$\max\limits_{i}|\lambda_{i}^{k+1}-\lambda_{i}^{k}|<\text{TOL}\|\boldsymbol{\lambda}\|_{\infty}$.
The experiment is repeated with
$\hat{d}_{\text{try}}\in\\{\hat{d}_{\text{try}}^{*},1.1\hat{d}_{\text{try}}^{*},1.2\hat{d}_{\text{try}}^{*}\\}$,
with $\hat{d}_{\text{try}}^{*}$ defined in (31) and
$\text{TOL}\in\\{10^{-2},5\cdot 10^{-3},5\cdot 10^{-6}\\}$. By viewing the
statistics in Figure 5(a), one can conclude that only the elements of
$\boldsymbol{\lambda}$ corresponding to particle pairs that would overlap
without a contact force correction get assigned a larger contact force and all
the other magnitudes are small, for all tested combinations of
hyperparameters. The choice of TOL will affect the magnitudes of the computed
contact forces to a very small extent as long as TOL is small enough. One
would expect that with a more restrictive TOL, the number of non-negligible
contact forces is reduced, but on the other hand, the number of iterations
required to solve (15) is increased. In practice, the dependence of TOL for
the magnitude of the contact foce is however very small, see Figure 5(a). In
Figure 5(b), the histogram shows the contact force magnitudes only for the
pairs flagged with $\hat{d}_{\text{try}}>\hat{d}_{\text{try}}^{*}$ that are
not flagged with $\hat{d}_{\text{try}}=\hat{d}_{\text{try}}^{*}$. Even if more
particle pairs are flagged within the time-step with a larger
$\hat{d}_{\text{try}}$, we show that all the extra contact forces computed
with a larger buffer region are small in magnitude. In Table 1 the mean
squared deviation and max relative difference of the three components of the
contact force for any particle in 200 configurations are presented versus the
setting with $\hat{d}_{\text{try}}=\hat{d}_{\text{try}}^{*}$ and
$\text{TOL}=5\cdot 10^{-6}$.

We end this section with a discussion on the choice of objective function. If
we redo the experiment with $\text{TOL}=5\cdot 10^{-6}$ and
$\hat{d}_{\text{try}}\in\\{\hat{d}_{\text{try}}^{*},1.1\hat{d}_{\text{try}}^{*},1.2\hat{d}_{\text{try}}^{*}\\}$,
but change the objective function to the approximated dissipative energy
$\boldsymbol{\lambda}^{T}\boldsymbol{\mathcal{D}}^{T}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$,
the “extra” flagged particle pairs for larger $\hat{d}_{\text{try}}$ will get
assigned a larger contact force magnitude, relative to choosing the 1-norm
objective function. The force magnitudes for these particle pairs are reported
in Figure 6 and can be compared to Figure 5(b).

(a) Bars display statistics for the contact force magnitudes with different
buffer regions as given by $\hat{d}_{\text{try}}$ and stopping criteria
$\max\limits_{i}|\lambda_{i}^{k+1}-\lambda_{i}^{k}|<\text{TOL}\|\boldsymbol{\lambda}\|_{\infty}$.
Whiskers show the minimum and maximum contact force magnitudes relative to
$\|\boldsymbol{f}_{i}\|$, the upper box edge the top $2.5\%$ of force
magnitudes, the lower box edge the bottom $1\%$ and the line inside the box,
the median.

(b) Histogram of the contact force magnitudes for the contact pairs with
$\hat{d}_{\text{try}}>\hat{d}_{\text{try}}^{*}$ not flagged for collision with
$\hat{d}_{\text{try}}=\hat{d}_{\text{try}}^{*}$. All are small in magnitude.

Figure 5: Almost all contact forces are small in magnitude relative to the
external force, which implies that only the pair or pairs that actually
violate the set minimum allowed distance $\hat{d}$ are assigned a significant
contact force, even if more particle pairs are flagged to be part of the
contact optimisation. For each hyperparameter combination, statistics is
collected from all flagged contacts in 200 random configurations of spheres
with packing density $24\%$. The same behaviour is noted for all
hyperparameter combinations.

| (a) MSD of contact force
---
$\hat{d}_{\text{try}}$ / TOL | $1\cdot 10^{-2}$ | $5\cdot 10^{-3}$ | $5\cdot 10^{-6}$
$\hat{d}_{\text{try}}^{*}$ | $2.45\cdot 10^{-6}$ | $6.63\cdot 10^{-7}$ | Reference
$1.1\hat{d}_{\text{try}}^{*}$ | $3.76\cdot 10^{-3}$ | $9.60\cdot 10^{-6}$ | $8.84\cdot 10^{-6}$
$1.2\hat{d}_{\text{try}}^{*}$ | $3.86\cdot 10^{-3}$ | $1.03\cdot 10^{-5}$ | $9.72\cdot 10^{-6}$
| (b) Max (relative) difference of contact force
---
$\hat{d}_{\text{try}}$ / TOL | $1\cdot 10^{-2}$ | $5\cdot 10^{-3}$ | $5\cdot 10^{-6}$
$\hat{d}_{\text{try}}^{*}$ | $0.14$ | $0.14$ | Reference
$1.1\hat{d}_{\text{try}}^{*}$ | $0.15$ | $0.16$ | $0.16$
$1.2\hat{d}_{\text{try}}^{*}$ | $0.16$ | $0.17$ | $0.17$
Table 1: Mean squared deviation and max difference of contact forces relative
to the reference with $\hat{d}_{\text{try}}=\hat{d}_{\text{try}}^{*}$ and
$\text{TOL}=5\cdot 10^{-6}$ with different hyperparameter combinations
($\hat{d}_{\text{try}}$, TOL) computed over all three components of the
contact force in for all spheres in 200 random configurations.  Figure 6:
Histogram of the contact force magnitudes for the contact pairs with
$\hat{d}_{\text{try}}>\hat{d}_{\text{try}}^{*}$ not flagged for collision with
$\hat{d}_{\text{try}}=\hat{d}_{\text{try}}^{*}$ using the objective function
$\boldsymbol{\lambda}^{T}\boldsymbol{\mathcal{D}}\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$
and the strice stopping criterion $\text{TOL}=5\cdot 10^{-6}$. A lot of these
extra contact forces are not small in magnitude, which would be expected, and
we therefore conclude that this choice of objective function is not as robust
as $\sum_{i}\lambda_{i}$.

###  3.2 Rods

####  3.2.1 Geometric considerations

Rods constituting of a central cylindrical part and semi-spherical caps at
both ends are considered. The multiblob grid is set via a matching strategy
for the mobility coefficients of a single particle as described in [2]. The
surface of the rod is everywhere one radius from its center line segment and
inter-particle distances $d^{\text{p}}_{i}$ can hence be determined by
computing the shortest distance between two line segments and subtracting
$2R_{\text{rod}}$, see [38, 39]. The surface-node-to-surface-node distances
$d^{\text{s}}_{ik}$ are determined by computing the shortest distance between
each source grid node belonging to an upsampled grid on particle 1 and the
line segment of particle 2. A new target surface node is introduced on
particle 2 where the vector of shortest distance cuts the surface of the
particle. This source/target pair is added to the list of surface points that
will be used to form the elements in the matrix $\boldsymbol{\mathcal{D}}$ in
(23). The procedure is repeated for all points sufficiently close to the other
particle on the upsampled surfaces of both particles. Finally, the closest
points of contact as defined by $d_{i}^{\text{p}}$ are added to the list not
to underestimate the closest distance between the particles, as for all grid
node pairs $k$, $d^{\text{s}}_{ik}\leq d_{i}^{\text{p}}$.

####  3.2.2 A chain of rods

Chains of rods of length $L_{\text{rod}}=0.5$ and radius
$R_{\text{rod}}=L_{\text{rod}}/20$ are considered, where for each chain, a
unit direction vector $\boldsymbol{t}\in\mathbb{R}^{3}$ and rotation
$(\theta,\phi)$ are drawn at random from the first orthant. The chain is
constructed with the first rod placed at the origin with orientation
coinciding with the $z$-axis. The consecutive rods are obtained from the
previous by rotating the particle by $(\theta,\phi)$ and then translating the
center coordinate by $\beta(\hat{d})\boldsymbol{t}$ in the coordinate frame of
the previous particle, with the constant $\beta(\hat{d})$ determining the
magnitude of the translation such that the smallest distance between a pair of
particles is $\hat{d}$. Example geometries are visualised in Figure 7(a).

Every rod in each chain is assigned a force in the direction of the next
particle so that a collision is caused in the next time-step ($\hat{d}$ is
violated) if no contact forces are applied. The given set of forces on the
chain configuration is then corrected with contact forces. Statistics for the
resulting particle-particle distances $d(\boldsymbol{\lambda})$, upon applying
contact forces, are reported versus $\hat{d}$ for a large range of $\hat{d}$
in Figure 7(b). Note that the largest $\hat{d}$ do not represent a realistic
choice for a dynamic simulation, but rather demonstrate the robustness of the
method and that the choices of $\hat{d}_{\text{try}}$ depend on $\Delta t$ and
therefore do not follow the curve for $\hat{d}$. By comparing
$d(\boldsymbol{\lambda})$ to the reference distance $\hat{d}$, it can be
concluded that the minimum allowed distance $\hat{d}$ is respected for almost
all particles in all chains. Statistics for
$\left(d(\boldsymbol{\lambda})-\hat{d}\right)/\hat{d}$ is illustrated for all
chains before and after contact in Figure 7(c). In Figure 7(d), histograms
display probable contact force magnitudes on any rod in any chain for fixed
$\hat{d}$ relative to the magnitude of the external force triggering the
contact avoiding algorithm. Contact forces are comparable in magnitude to the
external forces.

For reproducibility, the time-step size is set to $\Delta t=0.05$ and the
force magnitude is $\|\boldsymbol{f}_{i}\|=2\hat{d}^{1/3}$ (this particular
choice was made to promote collisions for the entire range of $\hat{d}$
considered in the test).

(a) Example rod chains with the same relative translation and rotation between
each consecutive pair of rods so that the inter-particle distance $\hat{d}$ is
equal between every pair. Red arrows indicate external forces towards the
center of the next particle that deliberately cause collisions during the next
time-step. Particle colors indicate the depth in the suspension.

(b) Distributions of the computed distance $d(\boldsymbol{\lambda})$ with
statistics collected from all rods in 100 chains for each $\hat{d}$, as
compared to the reference minimum allowed distance $\hat{d}$. For each
$\hat{d}$, the corresponding distance $\hat{d}_{\text{try}}$ is also
indicated, for which particles are flagged to be part of the collision
algorithm.

(c) Bars show the relative difference in computed distances
$d(\boldsymbol{\lambda})$ to $\hat{d}$ at the trial time-step without contact
forces and upon applying contact forces, for each choice of $\hat{d}$.
Whiskers display the minimum and maximum relative distance difference, box
edges the 10th and 90th percentile of the relative difference and the box
center line display the median.

(d) Magnitude of contact forces relative to the magnitude of the external
forces $\|\boldsymbol{f}_{i}\|$. All computed contact forces are comparable in
size to $\|\boldsymbol{f}_{i}\|$ or smaller.

Figure 7: Chains of 40 rods are forced to come into contact by externally
applied forces towards the next particle in the chain. By finding optimal
magnitudes of the contact forces, the minimum separation distance $\hat{d}$
can be kept. For each $\hat{d}$, the experiment is repeated with 100 randomly
generated chains.

####  3.2.3 Rods in a biaxial compression flow

Consider a 2D grid of 30 rods of length $L_{\text{rod}}=2$ and radius
$R_{\text{rod}}=L_{\text{rod}}/4$ arranged as in Figures 8(a)-8(b) and
affected by a background biaxial compression flow
$\boldsymbol{\mathcal{U}}_{\text{bg}}$ as indicated in the figure, given by

$\boldsymbol{\mathcal{U}}_{\text{bg}}=\begin{bmatrix}(\boldsymbol{E}\boldsymbol{x}_{1})^{T},\boldsymbol{0},(\boldsymbol{E}\boldsymbol{x}_{2})^{T},\boldsymbol{0},\dots,(\boldsymbol{E}\boldsymbol{x}_{n})^{T},\boldsymbol{0}\end{bmatrix}^{T}\quad\text{with
}\boldsymbol{E}=\begin{bmatrix}\dot{\gamma}&0&0\\\ 0&-\dot{\gamma}/2&0\\\
0&0&-\dot{\gamma}/2\end{bmatrix},\medspace\gamma=1.$ (32)

We choose different minimum allowed distances $\hat{d}$ and a hierarchy of
time-step sizes $\Delta t$ and discretise the dynamics of the system with
forward Euler with one hyperparameter setting $(\hat{d},\Delta t)$ at a time.
Whenever the smallest allowed distance $\hat{d}$ is violated at the end of the
time-step, contact forces are computed for particles considered sufficiently
close to contact at the trial time-step, as given by $\hat{d}_{\text{try}}$.
Due to the contractile nature of the background flow, the contact avoiding
algorithm will be triggered in almost every time-step, see Table 2. In Figure
8(c), the change in coordinate position per time-step is displayed for all
particles until the particles start to diverge along the $x$-axis (where the
background flow is diverging) and the simulation is stopped. The contact
forces cause no drastic jumps in the particle trajectories and a vast majority
of the coordinate updates are very close to zero. Table 2 indicates that
contact forces are robustly computed for a long sequence of time-steps. The
maximum deviation along the $x$-axis for the particles is displayed as
function of time in Figure 8(e). A smaller time-step allows for a slightly
smaller deviation than a larger time-step. Note however that for all
hyperparameters, particles stay in the $yz$-plane.

If all symmetries of the problem were kept, the particles should not only stay
in the same plane, but also keep the alignment with the $z$-axis. The
alignment can be quantified with the Onsager order parameter $S$ defined by
[43]

$S=\frac{1}{N}\sum_{i}^{N}\left\\{3/2\left(\boldsymbol{e}_{z}\cdot\boldsymbol{u}_{i}\right)^{2}-1/2\right\\}.$
(33)

If $S=1$, the particles are perfectly aligned with the $z$-axis and if
$S=-1/2$, particles are all perpendicular to the $z$-axis. The order parameter
is visualised for two choices of the parameter $\hat{d}$ in Figure 8(d). We
expect the order parameter to be $S(t)=1$ throughout the simulation, but
numerically this only holds up until some point in time when the alignment is
broken for some particles and $S(t)$ decreases. With varying choices of
$\hat{d}$ and $\Delta t$, the same qualitative behaviour can be noted for the
order parameter; The particles are initially ordered symmetrically and the
background flow is symmetric, then, particles first start to come in contact
in the vertical direction and later both vertically and horizontally. Due to
contact forces not being perfectly symmetric, the alignment is broken. This
divergence of the order-parameter however happens after a large number of
time-steps for all ($\hat{d},\Delta t$), from which we can conclude that the
contact handling is robust.

(a) Initial particle geometry in the $yz$-plane. Arrows indicate flow
direction.

(b) Initial particle geometry in the $xz$-plane. Arrows indicate flow
direction.

(c) Binned particle coordinate change per time-step with the large time-step
size $\Delta t=0.02$. The contact forces cause no drastic jumps in the
particle trajectories and all coordinate changes are small.

(d) The order parameter $S$ as defined in (33), using various time-step sizes
and two different choices of the minimum allowed distance $\hat{d}$. When all
particles are aligned with the $z$-axis, $S=1$. The alignment of the particles
is broken after some finite time. With the largest $\Delta t$, this happens
after approximately 100 time-steps using $\hat{d}=0.01R_{\text{rod}}$, at
$t\approx 3$.

(e) Maximal deviation in the $x$-coordinate: Despite the large number of time-
steps with the contact algorithm active, particles stay in the $yz$-plane for
a long time. The number of time-steps taken are reported in Table 2.

Figure 8: A grid of 30 rods affected by a biaxial compression flow. Contact forces are robustly computed whenever the minimum allowed distance $\hat{d}$ is violated. The symmetry of the configuration is broken but this happens after a large number of time-steps. Despite the challenging setting with the background flow pushing particles together, the particles stay in the original 2D-grid for a relatively long time. Time-step $\Delta t$ as fraction of $\Delta t^{*}=0.02$ |  # contact time-steps $\hat{d}=0.01R_{\text{rod}}$ $T_{\text{end}}=9$ ($T_{S}=4.36$) |  # contact time-steps $\hat{d}=0.1R_{\text{rod}}$ $T_{\text{end}}$ ($T_{S}$) |  total # steps $T_{\text{end}}$ ($T_{S}$)
---|---|---|---
$\Delta t^{*}$ | 410 (178) | 412 (180) | 450 (218)
$\Delta t^{*}/2$ | 819 (357) | 824 (362) | 900 (436)
$\Delta t^{*}/4$ | 1639 (715) | 1647 (723) | 1800 (872)
Table 2: Rods in a biaxial compression flow: Number of time-steps where the
contact algorithm is triggered for the different time-step sizes $\Delta t$,
together with the total number of steps taken, with $T_{S}$ the time up until
which the order parameter $S$ is displayed in Figure 8(d) and $T_{\text{end}}$
the time when particles start to diverge along the $x$-axis (where the
background flow is diverging).

###  3.3 Boomerangs

####  3.3.1 Geometric considerations

We consider boomerang particles of the type displayed in Figure 1(b). Such a
boomerang is constructed by dividing a rod of aspect ratio
$L_{\text{rod}}/R_{\text{rod}}=8$ and radius $R_{\text{rod}}=0.5$ as described
in Section 3.2.1 into two equal pieces, each with a cylindrical part and a
semi-spherical cap at one end. One of these half-rods is placed vertically and
one placed horisontally so that the corners of the cut area of the two pieces
touch. A circle is rotated around this point so that the two rod pieces are
joined into one particle. The center coordinate $\boldsymbol{x}$ is defined to
be the mass center of the particle.

For the barrier energy and the direction of contact forces and torques in
$\boldsymbol{\mathcal{D}}$, the contact distance $d^{\text{s}}$, based on
surface-node-to-surface-node distances, is used, as the particle geometry is
non-convex. For this purpose, the center curve of each particle is considered,
given by two line segments and a quarter of a circle joining the two segments.
For each point on the discretised surface of particle 1, we may easily compute
the shortest distance to the center curve of particle 2, subtracting one
particle radius to give the shortest distance between the surfaces. The
corresponding source/target pair of surface points is added to the list of
surface nodes used to construct the matrix $\boldsymbol{\mathcal{D}}$ in (23).
The procedure is repeated with reversed numbering of the two particles in the
pair.

####  3.3.2 A random suspension of boomerangs

Systems of 40 randomly positioned boomerangs in a cube of length $L=12$ are
considered, as exemplified in Figure 9(a). The boomerangs are generated one at
the time with a uniformly sampled center coordinate $\boldsymbol{x}$ and
quaternion $\boldsymbol{q}$ such that the distance to any other particle is
larger than $\delta=10^{-2}$. We perform the same type of test as for the
spheres in Section 3.1: For each generated configuration of 40 boomerangs, the
minimum separation distance is computed and $\hat{d}$ is chosen to be this
distance. The boomerangs are then assigned forces and torques randomly sampled
from a sphere with $\|\boldsymbol{\mathcal{F}}_{\text{ext}}\|=100$ so that in
a trial time-step with $\Delta t=0.01$, some particles violate the constraints
for surface-node-to-surface-node distances. For these particle pairs,
correcting contact forces are computed. Among all the particle pairs flagged
to potentially be assigned a contact force, the smallest distance obtained
with contact forces at time $t+\Delta t$ is reported vs $\hat{d}$ in Figure
9(b) and all distances for the flagged pairs are reported before and after the
correction with contact forces in Figure 9(c). It can here be concluded that
contact forces do not push particles apart unnecessarily far (note that the
parameter $\hat{d}_{\text{try}}$ is chosen large relative to $\hat{d}$ not to
miss any collisions). Despite the non-convex particle shape and the risk of
particles getting stuck in locked configurations, $\hat{d}$ can be maintained
for all the boomerangs in all of the 200 configurations.

(a) An example configuration, with colors indicating depth in the suspension.

(b) The minimum distances upon applying contact forces all respect the minimum
allowed distance $\hat{d}$. For each configuration, the corresponding
$\hat{d}_{\text{try}}$ is also displayed, determining the $N_{c}$ particle
pairs flagged to be part of the contact force optimisation.

(c) Inter-particle distances for the $N_{c}$ contact pairs before and after
applying contact forces, with statistics collected from all 200
configurations. Note that $\hat{d}_{\text{try}}$ is large relative to
$\hat{d}$, see (b).

Figure 9: Contact avoidance for 200 random configuration of 40 boomerangs.

##  4 A discussion on the computational cost

Despite the increased cost in any time-step with contacts, a lot can be gained
in terms of computational cost by applying the contact avoiding algorithm
presented in this paper, as a much larger time-step can be considered for
dynamical simulations than what would be allowed in a simulation without
contact forces – all time-step sizes used for the numerical experiments in
this work are excessively large, but the contact algorithm is still robust.

The main cost of finding optimal contact force magnitudes is the need for
determining the matrix vector product
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}$ in the computation of the
gradient of the barrier energy with respect to the vector of contact force
magnittudes, which is needed in the interior point method used for solving the
optimisation problem in (15). The matrix-matrix product
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}$ has to be computed and
stored once per time-step where contact occurs, as $\boldsymbol{\mathcal{M}}$
depends on $\boldsymbol{\mathcal{Q}}_{t}$ and $\boldsymbol{\mathcal{D}}$
depends on $\boldsymbol{\mathcal{Q}}_{t+\Delta t}^{*}$ (no dependence on
$\boldsymbol{\lambda}$). For the multiblob method, the mobility matrix can
easily be computed explicitly if the total number of particles is small. For
larger particle systems, evaluating
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}$ would amount to solving
$N_{c}$ Stokes mobility problems, the cost of which is determined by
$\hat{d}_{\text{try}}$ that will set $N_{c}$, which of course also depends on
the particle density in the system. One Stokes mobility problem also has to be
solved for every evaluation of the action of the contact forces,
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$, that
enters in the barrier energy. The number of such solves depends on the number
of iterations to solve the optimisation problem (15) with the interior-point
method, which, in turn, depends on the particle type, the time-step size and a
few hyperparameters as discussed in Section 2.3: Choices that need to be made
is e.g. how to regularise the barrier function, i.e. how to handle negative
distances between boundaries, and what stopping criteria to pick when solving
for minimum contact force magnitudes in (15). From numerical experiments in
this paper, it is expected that the number of mobility solves
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{D}}\boldsymbol{\lambda}$ at
least is in the range 10-20. The number of iterations is larger if the
stopping criteria is very strict and smaller if the stopping criteria is less
strict, which could be allowed by choosing a slightly larger buffer region
$\hat{d}$ than required for accuracy in $\boldsymbol{\mathcal{M}}$. We may
compare the number of iterations to what is reported in the literature for the
LCP technique in the work of Yan and collaborators, but at the same time
emphasise that the iteration count will vary depending on the setting: Yan et
al. report $\approx 20$ extra Stokes solves per time-step for spheres in [30]
and 5-10 iterations for dilute suspensions of spherocylinders in [29], but
$\mathcal{O}(1000)$ close to the random-close-packing limit.

The first-order interior-point method combined with the constrained minisation
problem considered in this work also has benefits over second order methods,
e.g. based on Newton’s method as applied to the non-constrained IPC-
formulation [20]. In contrast to an IPC-formulation, we are not dependent on
second order derivatives of the distances with respect to the particle
coordinates and also avoid solving large linear systems – other than the
Stokes mobility problem – while iteratively updating the solution vector in
the optimisation problem. The matrix in the linear system in an IPC-
formulation of the problem would be a function of $\boldsymbol{\mathcal{M}}$,
but also depend on the gradient of external forces and torques and on the
Hessian of the distances with respect to both particle center coordinates and
quaternions. The structure of this matrix is not trivially investigated for
the general case a priori and the matrix has to be projected onto the cone of
nonnegative matrices for Newton’s method to converge. Moreover, fast matrix-
vector techniques do not apply as for the Stokes mobility problem (applying
$\boldsymbol{M}\boldsymbol{\mathcal{F}}$, for some vector
$\boldsymbol{\mathcal{F}}$) [1]. Because of these difficulties, the current
formulation coupled to a first order method, such as the interior-point
method, is especially beneficial in very large particle systems, where we
neither want to compute $\boldsymbol{\mathcal{M}}$ explicitly nor want to
repeatedly solve a large linear system that is a non-trivial function of
$\boldsymbol{\mathcal{M}}$333An option for larger particle systems could
potentially be to solve the problem by employing the sparse nonlinear
optimiser provided in the SNOPT package [44]. Care however has to be taken for
how second order derivatives are computed..

##  5 Conclusions

We have presented an optimisation procedure to guarantee non-overlapping
configurations of 3D particles in an unbounded Stokesian fluid. The method is
based on a barrier formulation where non-overlapping constraints are rewritten
as a barrier energy, constrained to be zero for non-overlapping
configurations. Numerical examples are provided for spheres, rod-like
particles with semi-spherical caps and boomerangs. The numerical examples show
the performance of the contact force algorithm for dense suspensions of
particles deliberately pushed together by external forces or a background
flow. In all examples, the highlighted property of the proposed method is its
ability to keep a set minimum separation distance $\hat{d}$. Discrete
complementarity is obtained at the trial time-step between
$\boldsymbol{\lambda}$ and the modified barrier energy (with parameter
$\hat{d}_{\text{try}}$). We have numerically shown that the method is robust
for different choices of $\hat{d}$, particle shapes and flow scenarios. The
magnitudes of external forces, background flows and time-step sizes are in
this work deliberately chosen to force particles to come too close to each
other and in a realistic simulation, contact will occur less frequently and
$\hat{d}$ is supposed to be small. Strengths of the method is its simple
formulation, its ability of handling non-convex particles, the no-collision
guarantee of the formulation and the minimum impact on the system by contact
forces as these are balanced and their magnitudes are minimised. By solving
for force magnitudes explicitly, as we do here, a penalisation parameter can
be omitted, which otherwise has to be carefully tuned. Moreover, the size of
the optimisation problem is much more moderate (at least for moderatly crowded
systems) than if all particle coordinates are solved for implicitly, as in
Incremental Potential Contact (IPC) [20].

For future work we identify two directions:

  1. 1.

Approximation of the action of the mobility matrix: There is a global coupling
between all particles in the suspension in the mobility matrix,
$\boldsymbol{\mathcal{M}}$. Motivated by the fact that close interactions are
dominating, one option is to build an approximation to the action of
$\boldsymbol{\mathcal{M}}$ for matrix-vector multiplies
$\boldsymbol{\mathcal{M}}\boldsymbol{\mathcal{F}}$ constructed only from all
pairs of particles in contact, instead of applying the global mobility matrix.
An issue is however that a contact force implying no-collision determined with
such an approximation does not necessarily have to lead to a non-overlapping
configuration with $\boldsymbol{\mathcal{M}}$. A workaround could be to
introduce a small buffer for the minimum allowed distance so that
$\hat{d}\to\hat{d}(1+\delta)$ or to use the approximated solution vector for
contact force magnitudes as an initial guess for iterations with the full
mobility matrix. Along similar lines, separated clusters of particles could be
considered to compute contact forces locally.

  2. 2.

A continuous contact force potential: In this work, a discrete version of the
contact force potential is considered, computed from all points on the
particle surfaces sufficiently close to each other. One could also consider a
continuously defined potential for two particles in contact, where the barrier
energy at the trial time-step is expressed in terms of integrals over the
particle surfaces close to contact. The asymmetry of the contact forces noted
in the experiment with a biaxial compression flow in Section 3.2.3 has two
explanations: To start with, the discrete surface grid on the particles
results in different barrier energies for seemingly similar relative rotations
and particle-particle distances for different colliding pairs – depending on
the rotations of the particles around their own axis, a different number of
grid nodes may be flagged for collision. Furthermore, the shortest particle-
particle distance $d^{\text{p}}$ is sensitive to slight differences in
orientation when particles are close to parallel. One benefit of a continuous
contact force potential could be to better preserve symmetries. We leave this
approach for future work starting in 2D.

We end with a note on the optimisation algorithm: In this work, the inbuilt
optimisation procedure fmincon in Matlab has been used to solve the
optimisation problem for the contact force magnitude vector
$\boldsymbol{\lambda}$, employing a highly tuned interior point method. The
fact that a standard solver in Matlab can be used off-the-shelf is a strength
of the method. We are satisfied with a possibly _local_ optimum as long as a
zero barrier energy is reached, implying sufficiently separated particles.

## Acknowledgements

The authors thank Anders Forsgren for discussions on the choice of
optimisation algorithm, Georg Stadler for discussing alternative ways of
setting up the non-overlap constraints and Mattias Sandberg for constructive
input on the numerical examples provided in this work. We acknowledge the
support from the Swedish Research Council: grant no. 2019-05206 and the
research environment grant INTERFACE (biomaterials), no. 2016-06119.

## References

  * [1] F. B. Usabiaga, B. Kallemov, B. Delmotte, A. P. S. Bhalla, A. Donev, and B. E. Griffith, “Hydrodynamics of Suspensions of Passive and Active Rigid Particles : A Rigid Multiblob Approach,” _Comm. App. Math. Comp. Sci._ , vol. 11, no. 2, 2016 doi: 10.2140/camcos.2016.11.217
  * [2] A. Broms, M. Sandberg, and A.-K. Tornberg, “A locally corrected multiblob method with hydrodynamically matched grids for the stokes mobility problem,” 2022 doi: 10.48550/ARXIV.2207.11210
  * [3] S. Delong, F. Balboa Usabiaga, and A. Donev, “Brownian dynamics of confined rigid bodies,” _J. Chem. Phys_ , vol. 143, no. 14, 2015 doi: 10.1063/1.4932062
  * [4] B. Sprinkle, F. Balboa Usabiaga, N. A. Patankar, and A. Donev, “Large scale Brownian dynamics of confined suspensions of rigid particles,” _J. Chem. Phys_ , vol. 147, no. 24, 2017 doi: 10.1063/1.5003833
  * [5] Q. Brosseau, F. B. Usabiaga, E. Lushi, Y. Wu, L. Ristroph, J. Zhang, M. Ward, and M. J. Shelley, “Relating rheotaxis and hydrodynamic actuation using asymmetric gold-platinum phoretic rods,” _Phys. Rev. Lett_ , vol. 123, no. 17, p. 178004, 2019 doi: 10.1103/PhysRevLett.123.178004
  * [6] Q. Brosseau, F. B. Usabiaga, E. Lushi, Y. Wu, L. Ristroph, M. D. Ward, M. J. Shelley, and J. Zhang, “Metallic microswimmers driven up the wall by gravity,” _Soft Matter_ , vol. 17, no. 27, pp. 6597–6602, 2021 doi: 10.1039/d1sm00554e
  * [7] A. M. Fiore and J. W. Swan, “Fast Stokesian dynamics,” _J. Fluid Mech._ , vol. 878, pp. 544–597, 2019 doi: 10.1017/jfm.2019.640
  * [8] L. af Klinteberg and A. K. Tornberg, “Fast Ewald summation for Stokesian particle suspensions,” _Int. J. Numer. Methods Fluids_ , vol. 76, no. 10, pp. 669–698, 2014 doi: 10.1002/fld.3953
  * [9] B. J. Berne and P. Pechukas, “Gaussian model potentials for molecular interactions,” _J. Chem. Phys_ , vol. 56, no. 8, pp. 4195–4205, 1972 doi: 10.1063/1.1677837
  * [10] S. Varga and G. Jackson, “Study of the pitch of fluids of electrostatically chiral anisotropic molecules: Mean-field theory and simulation,” _Mol. Phys._ , vol. 104, no. 22-24, pp. 3681–3691, 2006 doi: 10.1080/00268970601058556
  * [11] Y. G. Tao, W. K. Den Otter, J. T. Padding, J. K. Dhont, and W. J. Briels, “Brownian dynamics simulations of the self- and collective rotational diffusion coefficients of rigid long thin rods,” _J. Chem. Phys_ , vol. 122, no. 24, 2005 doi: 10.1063/1.1940031
  * [12] D. Baraff, “Issues in computing contact forces for non-penetrating rigid bodies,” _Algorithmica_ , vol. 10, no. 2-4, pp. 292–352, 1993 doi: 10.1007/BF01891843
  * [13] Anitescu, Cremer, and Potra, “Formulating 3D Contact Dynamics Problems,” _Mechanics of Structures and Machines_ , vol. 24, no. 4, pp. 405–437, 1996 doi: 10.1080/08905459608905271
  * [14] A. Tasora and M. Anitescu, “A Fast NCP Solver for Large Rigid-Body Problems with Contacts, Friction, and Joints,” _Multibody Dynamics_ , pp. 45–55, 2008 doi: 10.1007/978-1-4020-8829-2_3
  * [15] A. Tasora, D. Negrut, and M. Anitescu, “Large-scale parallel multi-body dynamics with frictional contact on the graphical processing unit,” _Proc. Inst. Mech. Eng., Part K: J. Multi-Body Dyn_ , vol. 222, no. 4, pp. 315–326, 2008 doi: 10.1243/14644193JMBD154
  * [16] M. Tur, F. J. Fuenmayor, and P. Wriggers, “A mortar-based frictional contact formulation for large deformations using Lagrange multipliers,” _Comput. Methods Appl. Mech. Eng._ , vol. 198, no. 37-40, pp. 2860–2873, 2009 doi: 10.1016/j.cma.2009.04.007
  * [17] P. Wriggers, _Computational contact mechanics, second ed., Springer_ , 2006\. ISBN 9783540326083 doi: 10.1007/978-3-540-32609-0
  * [18] J. M. Snyder, “Interactive tool for placing curved surfaces without interpenetration,” _Proceedings of the ACM SIGGRAPH Conference on Computer Graphics_ , pp. 209–218, 1995 doi: 10.1145/218380.218444
  * [19] D. Harmon, D. Panozzo, O. Sorkine, and D. Zorin, “Interference-aware geometric modeling,” _ACM Trans. Graph._ , vol. 30, no. 6, p. 1, 2011 doi: 10.1145/2070781.2024171
  * [20] M. Li, Z. Ferguson, T. Schneider, T. Langlois, D. Zorin, D. Panozzo, C. Jiang, and D. M. Kaufman, “Incremental Potential Contact: Intersection- and Inversion-free, Large-Deformation Dynamics,” _ACM Trans. Graph._ , vol. 39, no. 4, 2020 doi: 10.1145/3386569.3392425
  * [21] Z. Ferguson, M. Li, T. Schneider, F. Gil-Ureta, T. Langlois, C. Jiang, D. Zorin, D. M. Kaufman, and D. Panozzo, “Intersection-free rigid body dynamics,” _ACM Trans. Graph._ , vol. 40, no. 4, pp. 1–16, 2021 doi: 10.1145/3476576.3476773
  * [22] M. Li, D. M. Kaufman, and C. Jiang, “Codimensional incremental potential contact,” _ACM Trans. Graph._ , vol. 40, no. 4, 2021 doi: 10.1145/3450626.3459767
  * [23] S. Yamamoto and T. Matsuoka, “A method for dynamic simulation of rigid and flexible fibers in a flow field,” _J. Chem. Phys_ , vol. 98, no. 1, pp. 644–650, 1993 doi: 10.1063/1.464607
  * [24] ——, “Dynamic simulation of fiber suspensions in shear flow,” _J. Chem. Phys_ , vol. 102, no. 5, pp. 2254–2260, 1995 doi: 10.1063/1.468746
  * [25] D. Das and D. Saintillan, “Electrohydrodynamic interaction of spherical particles under Quincke rotation,” _Phys. Rev. E_ , vol. 87, no. 4, pp. 1–14, 2013 doi: 10.1103/PhysRevE.87.043014
  * [26] B. Delmotte, E. Climent, and F. Plouraboué, “A general formulation of Bead Models applied to flexible fibers and active filaments at low Reynolds number,” _J. Comput. Phys_ , vol. 286, pp. 14–37, 2015 doi: 10.1016/j.jcp.2015.01.026
  * [27] E. Corona, L. Greengard, M. Rachh, and S. Veerapaneni, “An integral equation formulation for rigid bodies in Stokes flow in three dimensions,” _J. Comput. Phys_ , vol. 332, pp. 504–519, 2017 doi: 10.1016/j.jcp.2016.12.018
  * [28] L. Lu, “Parallel contact-aware algorithms for large-scale direct blood flow simulations,” Ph.D. dissertation, Courant Institute of Mathematical Sciences, New York University, 2019.
  * [29] W. Yan, H. Zhang, and M. J. Shelley, “Computing collision stress in assemblies of active spherocylinders: Applications of a fast and generic geometric method,” _J. Chem. Phys_ , vol. 150, no. 6, 2019 doi: 10.1063/1.5080433
  * [30] W. Yan, E. Corona, D. Malhotra, S. Veerapaneni, and M. Shelley, “A scalable computational platform for particulate Stokes suspensions,” _J. Comput. Phys_ , vol. 416, p. 109524, 2020 doi: 10.1016/j.jcp.2020.109524
  * [31] L. Bystricky, S. Shanbhag, and B. Quaife, “Stable and contact-free time stepping for dense rigid particle suspensions,” _Int. J. Numer. Methods Fluids_ , vol. 92, no. 2, pp. 94–113, 2020 doi: 10.1002/fld.4774
  * [32] W. Yan, S. Ansari, A. Lamson, M. A. Glaser, R. Blackwell, M. D. Betterton, and M. Shelley, “Toward the cellular-  scale simulation of motor-  driven cytoskeletal assemblies,” _eLife_ , no. 11, p. e74160, 2022 doi: 10.7554/eLife.74160
  * [33] R. Kohl, E. Corona, V. Cheruvu, and S. Veerapaneni, “Fast and accurate solvers for simulating Janus particle suspensions in Stokes flow,” pp. 1–24, 2021 doi: 10.48550/arXiv.2104.14068
  * [34] L. Lu, M. J. Morse, A. Rahimian, G. Stadler, and D. Zorin, “Scalable simulation of realistic volume fraction red blood cell flows through vascular networks,” _Int. Conf. High Perform. Comput. Netw. Storage Anal._ , 2019 doi: 10.1145/3295500.3356203
  * [35] L. Bystricky, “Contact-free Simulations of Rigid Particle Suspensions Using Boundary Integral Equations,” Ph.D. dissertation, Florida State University. ISBN 9780438447769 2018.
  * [36] L. Lu, A. Rahimian, and D. Zorin, “Contact-aware simulations of particulate Stokesian suspensions,” _J. Comput. Phys_ , vol. 347, pp. 160–182, 2017 doi: 10.1016/j.jcp.2017.06.039
  * [37] S. Boyd and L. Vandenberghe, _Convex optimization_ , 2009 doi: 10.1142/9789814412520_0010
  * [38] V. J. Lumelsky, “On fast computation of distance between line segments,” _Inf Process Lett_ , vol. 21, no. 2, pp. 55–61, 1985 doi: 10.1016/0020-0190(85)90032-8
  * [39] Ondrej, “Fast shortest distance between two line segments (in n dimensions),” MATLAB Central File Exchange, Retrieved May 4, 2022, https://www.mathworks.com/matlabcentral/fileexchange/49502-fast-shortest-distance-between-two-line-segments-in-n-dimensions.
  * [40] A. Fischer, “A special newton-type optimization method,” _Optimization_ , vol. 24, no. 3-4, pp. 269–284, 1992 doi: 10.1080/02331939208843795
  * [41] Y. H. Dai and R. Fletcher, “Projected Barzilai-Borwein methods for large-scale box-constrained quadratic programming,” _Numerische Mathematik_ , vol. 100, no. 1, pp. 21–47, 2005 doi: 10.1007/s00211-004-0569-y
  * [42] R. Fletcher, “On the barzilai-borwein method,” in _Optimization and Control with Applications_ , L. Qi, K. Teo, and X. Yang, Eds. Boston, MA: Springer US, 2005 doi: 10.1007/b104943. ISBN 978-0-387-24255-2 pp. 235–256.
  * [43] H. Doi, K. Z. Takahashi, K. Tagashira, J. ichi Fukuda, and T. Aoyagi, “Machine learning-aided analysis for complex local structure of liquid crystal polymers,” _Sci. Rep._ , vol. 9, no. 1, pp. 1–12, 2019 doi: 10.1038/s41598-019-51238-1
  * [44] P. E. Gill, W. Murray, and M. A. Saunders, “Snopt: An sqp algorithm for large-scale constrained optimization,” _SIAM Rev._ , vol. 47, no. 1, pp. 99–131, 2005 doi: 10.1137/S0036144504446096