# An Efficient Phase-field Framework for Contact Dynamics between Deformable
Solids in Fluid Flow

Biswajeet Rath<EMAIL_ADDRESS>Xiaoyu Mao<EMAIL_ADDRESS>Rajeev K. Jaiman<EMAIL_ADDRESS>Department of Mechanical Engineering,
The University of British Columbia, Vancouver, BC V6T 1Z4

###### Abstract

Elastic contact in hydrodynamic environments is a complex multiphysics
phenomenon and can be found in applications ranging from engineering to
biological systems. Understanding the intricacies of this problem requires the
development of a generalized framework capable of handling topological changes
and transitioning implicitly from fluid-structure interaction (FSI) conditions
to solid-solid contact conditions. We propose a mono-field interface advancing
method for handling multibody contact simulations in a submerged environment.
Given the physical demands of the problem, we adopt a phase-field based fully
Eulerian approach to resolve the multiphase and multibody interactions in the
system. We employ a stabilized finite element formulation and a partitioned
iterative procedure to solve the unified momentum equation comprising solid
and fluid dynamics coupled with the Allen-Cahn phase-field equation. The
evolution of solid strain in the Eulerian reference frame is governed by the
transport of the left Cauchy Green deformation tensor. We introduce a contact
force approach to handle smooth elastic-elastic and elastic-rigid contact
based on the overlap of the diffused interfaces of two colliding bodies. The
diffuse interface representation ameliorates the issue of discontinuous
contact forces generally found in conventional approaches. The implemented
model is verified for smooth dry contact and FSI contact scenarios. We obtain
a spatial convergence rate of around $1.8$ for the Hertzian contact case
involving a deformable cylinder on an elastic plane. Relative volume error is
found to be within $5\%$ for the FSI contact case. The developed approach
eliminates the need for solving multiple phase-field equations or
reconstructing the solid boundaries at every time step. This reduces the
overall computational time by nearly $10\%$ compared to an explicit interface
extraction method. Using the proposed framework, we demonstrate the collision
dynamics between multiple bodies submerged in an open tank.

Keywords: Phase-field model, Fully Eulerian FSI, Contact modeling, Multibody
contact

##  1 Introduction

Fluid-Structure interaction (FSI) is an omnipresent and complicated phenomenon
that manifests itself in natural, biological and engineering systems. The
relative motion between a solid and the fluid in both internal and external
flows leads to the application of loads on the solid resulting in deformations
and/or displacements in it. Our primary motivation comes from ship-ice
interaction in Arctic regions during the collision of ships with large ice
sheets. The framework should be able to identify different types of interfaces
and handle various structural level interactions apart from free-surface
motion and fluid-structure interaction. Simulating the process of ship-ice
interactions involves handling drifting, large deformation, collisions,
fracture and fragmentation of the ice sheets [1, 2, 3]. In such scenarios,
accurately capturing contact dynamics with hydrodynamic interactions poses
significant challenges for computational modeling. Typically, solid-solid
contact problems are predominantly modeled within a Lagrangian framework,
where the discretization moves with the deforming objects. By considering an
Eulerian formulation, our aim is to develop an integrated and versatile model
capable of effectively handling contact dynamics with fluid-structure
interaction.

In Lagrangian formulations, numerical modeling of contact is a challenging
task for two important reasons. First, both the contact boundary and the
contact boundary condition i.e. contact force are unknown a priori and are a
part of the solution. The second reason is the discontinuous nature of the
contact force which introduces an inherent irregularity and non-smoothness in
the problem. Mathematically, these issues translate into a set of constraints
known as the Hertz-Signorini-Moreau (HSM) conditions [4]. They consist of a
non-penetration constraint and a non-zero contact force condition. This
results in a set of variational inequalities for the overall problem [4, 5].
One of the earliest approaches for modeling contact for large deformation
bodies involving non-matching meshes in a Lagrangian frame of reference was
the node-to-surface algorithm introduced by Hughes et al. [6] and improved in
some of the later works in [7, 8]. The method is based on the idea of checking
for contact between a node on the slave body and a surface element on the
master body. The contact constraints are then imposed onto this set of slave
nodes (also known as active set) based on the type of solution method i.e.
interior point method, Lagrange multiplier method, penalty method, etc. The
choice of an appropriate slave body and concentrated nodal forces were some of
the limitations of this approach.

Several weak (integral) implementations of the contact constraints came up
later to bypass the drawbacks of the node-to-surface algorithm. The most
significant among them was the surface-to-surface contact approach [8, 9]. The
contact contribution to the weak form is integrated by a set of a priori
defined Gauss quadrature points on each element of the slave surface. This
algorithm also suffers from the master/ slave bias similar to the node-to-
surface method. However, the improved two-pass versions of the above
algorithms, which treat both bodies as slaves in two alternate iterations,
remove this inherent bias [10, 11]. Mortar methods [12, 13, 14] are another
class of approaches for resolving contact in a weak sense as they allow for a
variationally consistent treatment of non-penetration and frictional sliding
conditions. Recently, most of the above methods have been implemented in
isogeometric formulations (IGA) [15, 16] as well for their improved accuracy
and continuity properties. Authors in [17, 18, 19, 20] illustrate the
applicability of IGA for problems involving contact mechanics and impact
dynamics.

In the context of FSI applications, we have the additional traction continuity
condition at the fluid-solid interface. Thus, we need to ensure a smooth
transition between the FSI conditions and the contact conditions over the
solid boundary. The presence of the fluid and the requirement to handle large
solid deformations make the Eulerian frame of reference more suitable for our
applications. Keeping in mind the physical requirements of the problem, we
utilize a fully Eulerian-based FSI framework for the present work, whereby the
space itself is discretized and the evolving solid bodies are transported
through it. A fully Eulerian approach alleviates the issues of dissimilar
kinematics and dynamics between the fluid and solid domains by formulating the
governing conservation laws in a unified framework. The Eulerian approach can
significantly reduce the complexity of contact detection and the satisfaction
of contact constraints. Some of the recent studies have demonstrated the
simplicity of handling solid-solid contact both in isolated [21, 22] and
submerged environments [23, 24] within a fully Eulerian framework.

One of the obvious ways to handle the non-penetration constraint in contact is
to use optimization algorithms for solving the quadratic program. Such an
approach was carried out in [21] where they used the Interior Point method in
computer graphics simulations of colliding bodies. They utilized a finite
difference fully Eulerian solver to carry out the FSI simulations and used the
reference map to evaluate the solid strains. The velocities of the solids were
updated separately to keep track of the individual objects and a collision
detection algorithm [25] was used for broad-phase contact handling. However,
using optimization algorithms is ill-suited and tedious for complex domain
discretizations and diffuse interface surfaces. Hence, in a later work [26] by
the above groups they shifted to an explicit impulse-based approach for
handling contact. An explicit contact force approach is the most common method
for dealing with collisions in fully Eulerian FSI simulations. Such an
explicit force can be represented in a generalized form as:

$\boldsymbol{f}_{c}=\kappa\delta(x),$ (1)

where $\delta(x)$ is the short-range influence function that is dependent on
the gap between the bodies and $\kappa$ is a constant used to scale the force.
Different forms of the influence function have been used in literature. [23]
used a smooth Heaviside function of the relative level set as their influence
function to compute the contact force. Authors in [27, 24] formulated their
influence function as inversely proportional to the gap between the colliding
bodies. Such a form of the forcing function ensures that true contact never
occurs and a fluid layer is always present between the colliding bodies.
However, the contact force in the above model can have very high values and
might cause stability issues in the numerical solver. We propose a contact
force model similar to the above ideas but based on the phase-field model.
This allows us to model the contact force as an additional body force as
opposed to surface traction in sharp interface methods.

Using a fully Eulerian framework, we adopt the phase-field model for the
spatio-temporal evolution of the implicit interface. The phase-field
interpolation and the unified formulation implicitly satisfy the velocity and
stress continuity conditions at the fluid-solid interface. Through stabilized
finite element discretization, we implement a phase-field based contact force
model to handle collisions between solids in FSI applications. This alleviates
the issue of a discontinuous contact force, since now the contact force can
build up gradually as the diffuse boundaries of the colliding bodies overlap
each other. We verify the implemented contact routine with the traditional
Hertz contact model for smooth elastic solids. Such verification studies are
important for ensuring the accuracy of the contact model but are usually
absent in FSI contact works. We further verify our model for elastic-elastic
contact in submerged environments. Finally, we present a novel method for
simulating multibody contact problems for solids with identical physical
properties. The proposed method obviates the need to define multiple phase
functions and strain variables for each solid object, thus saving significant
computational resources otherwise required for solving the extra PDEs.

The present article is structured as follows. In Section 2, we briefly review
the fully Eulerian formulation and present the governing equations deriving
from the conservation laws. We also discuss the accurate evaluation of solid
stresses in a fully Eulerian framework. Section 3 discusses the background of
generalized contact simulations. It builds from the classical contact
constraints and subsequently dives into diffused contact handling for FSI
problems. Multibody contact modeling is introduced in Section 4 along with the
algorithm for the same. Details of the numerical implementation are presented
in Section 5. Section 6 focuses on the verification and convergence studies of
the implemented model for smooth dry contact and FSI contact scenarios. In
Section 7, we present results obtained from two applications to demonstrate
the robustness of the proposed method. Section 8 summarizes the overall
conclusions and provides directions for future work.

##  2 Fully Eulerian Formulation

In this section, we provide the background for the fully Eulerian formulation.
Both the solid and fluid equations are represented in the Eulerian
coordinates. From standard conventions in continuum mechanics, let $\Omega$
denote the domain in the reference configuration $\boldsymbol{X}$ at $t=0$ and
$\Omega^{t}$ denote the deformed state of the domain in spatial coordinates,
$\boldsymbol{x}$. The fluid and the solid make up the two parts of the domain,
such that $\Omega=\Omega_{f}^{t}\cup\Omega_{s}^{t}$ and
$\phi=\Omega_{f}^{t}\cap\Omega_{s}^{t}$. We consider both the solid and the
fluid to be incompressible. We analyze the interactions between a Newtonian
fluid and a hyperelastic solid for all our FSI cases.

###  2.1 Continuum Equations

We can express the unified continuum equations for the physical domains via a
phase indicator $(\phi(\boldsymbol{x},t))$ based interpolation as follows

$\displaystyle\nabla\cdot\boldsymbol{v}$ $\displaystyle=0\ \ $
$\displaystyle\text{on }\Omega,$ (2)
$\displaystyle\rho(\phi)\left(\frac{\partial\boldsymbol{v}}{\partial
t}\bigg{\rvert}_{x}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}\right)$
$\displaystyle=\nabla\cdot(\boldsymbol{\sigma}(\phi))+\boldsymbol{b}(\phi)\ \
$ $\displaystyle\text{on }\Omega,$

where the physical properties are defined by the following interpolation
rules:

$\displaystyle\rho(\phi)$ $\displaystyle=\phi_{f}\rho_{f}+\phi_{s}\rho_{s},$
(3) $\displaystyle\mu(\phi)$ $\displaystyle=\phi_{f}\mu_{f}+\phi_{s}\mu_{s}.$

In Eq. (2), $\boldsymbol{v}$ denotes the velocity and $\boldsymbol{b}$ denotes
the body force at each spatial point $\boldsymbol{x}\in\Omega$. The first
equation ensures mass conservation of the incompressible continuua in the
domain in Eulerian coordinates. The second equation is the momentum balance
equation for the combined fluid-solid domain. We vary the order parameter in
the range $[-1,1]$ in the current work. Thus, we have
$\phi_{f}=\frac{1-\phi}{2}$ and $\phi_{s}=\frac{1+\phi}{2}$ as the
interpolation functions in Eq. (3). We expand more on the phase indicator
function in Section 2.4. Interested readers are referred to our previous works
in [28, 29] to review the process of arriving at the above set of equations in
an Eulerian reference frame.

The semi-discrete variational form for the above differential equation, in
conjunction with the generalized-alpha time integration method [30, 31] can be
presented as follows. Let $\mathcal{S}^{h}$ be the space of trial solutions,
whose values satisfy the Dirichlet boundary conditions and $\mathcal{V}^{h}$
be the space of test functions whose values vanish on the Dirichlet boundary.
Thus, find $[\boldsymbol{v}_{h}^{n+\alpha},p_{h}^{n+1}]\in\mathcal{S}^{h}$
such that $\forall[\boldsymbol{\psi}_{h},q_{h}]\in\mathcal{V}^{h}$,

$\displaystyle\int_{\Omega}\rho(\phi)\left(\partial_{t}\boldsymbol{v}_{\mathrm{h}}^{\mathrm{n}+\alpha_{\mathrm{m}}}+(\boldsymbol{v}_{\mathrm{h}}^{\mathrm{n}+\alpha}\cdot\nabla)\boldsymbol{v}_{\mathrm{h}}^{\mathrm{n}+\alpha}\right)\cdot\boldsymbol{\psi}_{\mathrm{h}}\mathrm{d}\Omega+\int_{\Omega}\boldsymbol{\sigma}_{\mathrm{h}}^{\mathrm{n}+\alpha}:\nabla\boldsymbol{\psi}_{\mathrm{h}}\mathrm{d}\Omega+\int_{\Omega}q_{\mathrm{h}}\left(\nabla\cdot\boldsymbol{v}_{\mathrm{h}}^{\mathrm{n}+\alpha}\right)\mathrm{d}\Omega$
(4)
$\displaystyle+\sum_{\mathrm{e}=1}^{\mathrm{n}_{\mathrm{el}}}\int_{\Omega_{\mathrm{e}}}\frac{\tau_{\mathrm{m}}}{\rho(\phi)}\left(\rho(\phi)(\boldsymbol{v}_{\mathrm{h}}^{\mathrm{n}+\alpha}\cdot\nabla)\boldsymbol{\psi}_{\mathrm{h}}+\nabla
q_{\mathrm{h}}\right)\cdot\boldsymbol{\mathcal{R}}_{\mathrm{m}}(\boldsymbol{v},p)\mathrm{d}\Omega_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{\mathrm{n}_{\mathrm{el}}}\int_{\Omega_{\mathrm{e}}}\nabla\cdot\boldsymbol{\psi}_{\mathrm{h}}\tau_{\mathrm{c}}\rho(\phi)\mathcal{R}_{\mathrm{c}}(\boldsymbol{v})\mathrm{d}\Omega_{\mathrm{e}}$
$\displaystyle=\int_{\Omega}\boldsymbol{b}\left(t^{\mathrm{n}+\alpha}\right)\cdot\boldsymbol{\psi}_{\mathrm{h}}\mathrm{d}\Omega+\int_{\Gamma_{\mathrm{h}}}\boldsymbol{h}\cdot\boldsymbol{\psi}_{\mathrm{h}}\mathrm{d}\Gamma.$

The first line consists of the Galerkin terms for the combined momentum and
continuity equations. The second line contains the Petrov-Galerkin
stabilization terms for the continuum equations.
$\boldsymbol{\mathcal{R}}_{m}$ and $\mathcal{R}_{c}$ are the element-wise
residuals for the momentum and continuity equations respectively. The
stabilization parameters $\tau_{m}$ and $\tau_{c}$ [32, 33] are given by

$\tau_{m}=\left[\left(\frac{2}{\Delta
t}\right)^{2}+\boldsymbol{v}_{h}\cdot\boldsymbol{G}\boldsymbol{v}_{h}+C_{I}\bigg{(}\frac{\mu(\phi)}{\rho(\phi)}\bigg{)}^{2}\boldsymbol{G}:\boldsymbol{G}\right]^{-1/2},\text{
}\tau_{c}=\frac{1}{\mathrm{tr}(\boldsymbol{G})\tau_{m}},$ (5)

where $C_{I}$ is a constant derived from the element-wise inverse estimate and
$\boldsymbol{G}$ is the element contravariant metric tensor. We elaborate on
the strain evolution in the solid domain in the next sub-section.

###  2.2 Solid Strain Evolution

As we do not solve for the solid displacements in an Eulerian setting, we need
an additional equation to evolve the solid strains and close the system of
equations. This obviates the need for a mesh motion equation and helps us
compute the stresses in the structures accurately. We capture the solid
deformation by the evolution of the left Cauchy-Green deformation tensor
($\boldsymbol{B}$). Some of the other approaches used in the literature are
the transport of the reference/ inverse map function ($\boldsymbol{\xi}$) [23,
34] or the deformation gradient tensor ($\boldsymbol{F}$) [35, 36]. Although
$\boldsymbol{\xi}$ has the additional advantage of storing the location
information, it can be inconvenient for actual stress computations as it
involves second-order derivatives. For the choice between $\boldsymbol{F}$ and
$\boldsymbol{B}$, we go with $\boldsymbol{B}$ as it is a symmetric tensor and
saves us the effort of computing four different components in 2-D or nine
different components in 3-D. We emphasize the decomposition of the stress
tensor into its volumetric and deviatoric components to ensure a uniform
definition of pressure in the domain.

The transport equation for $\boldsymbol{B}$ in a solid moving with velocity
$\boldsymbol{v}_{s}$ is given as

$\frac{\partial\boldsymbol{B}}{\partial
t}+(\boldsymbol{v}_{s}\cdot\nabla)\boldsymbol{B}=\nabla\boldsymbol{v}_{s}\boldsymbol{B}+\boldsymbol{B}(\nabla\boldsymbol{v}_{s})^{T}\
\ \ \ \text{on }\Omega^{s}.$ (6)

The velocity of a node in the solid domain can be expressed as
$\boldsymbol{v}_{s}=\left(\frac{1+\phi}{2}\right)\boldsymbol{v}$. Such an
interpolation ensures that the above transport equation is solved only in the
solid domain and smoothly extended into the interface. There is no physical
significance of defining a strain measure in the fluid domain. To avoid any
instabilities that might arise from the fluid side and creep up into the
solid, we impose a Dirichlet condition of $\boldsymbol{B}=\boldsymbol{I}$ in
the fluid domain at a cutoff for the phase-field function of $\phi<-0.95$
[37].

The weak form for the transport of the left Cauchy-Green deformation tensor is
presented below. Let $\mathcal{S}^{h}$ be the space of trial solutions, whose
values satisfy the Dirichlet boundary conditions and $\mathcal{V}^{h}$ be the
space of test functions whose values vanish on the Dirichlet boundary. The
variational form of the left Cauchy-Green tensor equation can be written as:
find $\boldsymbol{B}_{h}^{n+\alpha}\in\mathcal{S}^{h}$ such that
$\forall\boldsymbol{m}_{h}\in\mathcal{V}^{h}$,

$\displaystyle\int_{\Omega}(\boldsymbol{m}_{h}):\bigg{(}\partial_{t}\boldsymbol{B}_{h}^{\mathrm{n+\alpha_{m}}}+\bigg{(}\frac{1+\phi}{2}\bigg{)}(\boldsymbol{v}^{{\mathrm{n+\alpha}}}\cdot\nabla)\boldsymbol{B}_{h}^{\mathrm{n+\alpha}}\bigg{)}\mathrm{d}\Omega-$
(7)
$\displaystyle\int_{\Omega}\bigg{(}\bigg{(}\frac{1+\phi}{2}\bigg{)}\nabla\boldsymbol{v}^{{\mathrm{n+\alpha}}}\boldsymbol{B}_{h}^{\mathrm{n+\alpha}}+\bigg{(}\frac{1+\phi}{2}\bigg{)}\boldsymbol{B}_{h}^{\mathrm{n+\alpha}}(\nabla\boldsymbol{v}^{{\mathrm{n+\alpha}}})^{T}\bigg{)}:(\boldsymbol{m}_{h})\mathrm{d}\Omega=0,$

where $(:)$ represents double dot product. We obtain updated values of
$\boldsymbol{B}$ by solving the above variational equation in the solid domain
which can then be used to calculate the new solid stresses before substituting
in the unified continuum equations in the next time step.

###  2.3 Stress Computation for Incompressible Materials

In this sub-section, we delve deeper into evaluating the stress for
incompressible materials in a fully Eulerian framework. The total stress in
the domain in Eq 2 can be expressed as a weighted average of the fluid and
solid stresses based on the phase-field function:

$\boldsymbol{\sigma}=\phi_{f}\boldsymbol{\sigma}_{f}+\phi_{s}\boldsymbol{\sigma}_{s}.$
(8)

Since we are specifically concerned with incompressible materials in the
present article, it is preferable to compute the pressure field as a separate
independent variable. This calls for a decomposition of the stress fields in
the physical domains into their corresponding volumetric and deviatoric
components. This step will allow us to have a single unified pressure field
for the entire domain and two separate deviatoric stresses for the fluid and
solid respectively as shown:

$\boldsymbol{\sigma}=-p\boldsymbol{I}+\phi_{f}\boldsymbol{\sigma}_{f}^{\prime}+\phi_{s}\boldsymbol{\sigma}_{s}^{\prime},$
(9)

where $p$ is the unified pressure field and $(^{\prime})$ indicates deviatoric
component of the stress tensor, given by
$\boldsymbol{A}^{\prime}=\boldsymbol{A}-\frac{1}{d}tr(\boldsymbol{A})\boldsymbol{I}$
for a tensor $\boldsymbol{A}$, where $d$ is the dimension of the problem. The
deviatoric components in the fluid and solid domains can be evaluated as
follows:

$\boldsymbol{\sigma}_{f}^{\prime}=\mu_{f}(\nabla\boldsymbol{v}_{f}+(\nabla\boldsymbol{v}_{f})^{T}),$
(10)

and

$\boldsymbol{\sigma}_{s}^{\prime}=\mu_{s}(\nabla\boldsymbol{v}_{s}+(\nabla\boldsymbol{v}_{s})^{T})+\mu_{s}^{L}(\boldsymbol{B}-\frac{1}{d}tr(\boldsymbol{B})\boldsymbol{I}),$
(11)

where we have used the incompressible Neo-Hookean model for the hyperelastic
solid. Such a decomposition ensures that the definition of the Lagrange
multiplier for enforcing incompressibility i.e. pressure is uniform across the
domain. During the numerical implementation, we need to compute the Jacobian
of the above stress term with respect to velocity to plug in the continuum
equation. To achieve this, we utilize the transport equation of the left
Cauchy-Green deformation tensor $\boldsymbol{B}$ as elaborated in [29]. We
discuss more on the Jacobian computations in A. This completes the system of
equations for the fluid-solid continuum. We focus on the implicit interface
capturing method in the next sub-section.

###  2.4 Diffuse Interface Representation

We utilize the phase-field model as a mathematical construct to evolve the
interface implicitly during the FSI simulation. The phase-field function/
order parameter $\phi(\boldsymbol{x},t)$ varies sharply within the interface
and has constant values outside of it. The Allen-Cahn equations for phase-
field description can be derived from the process of free energy minimization,
where the Ginzburg-Landau free energy functional is given by:

$\mathcal{E}:\mathcal{H}^{1}(\Omega)\cap\mathcal{L}^{4}(\Omega)\to\mathbb{R}_{\geq
0},\mathcal{E}(\phi(\boldsymbol{x},t))=\int_{\Omega}\left(F(\phi(\boldsymbol{x},t))+\frac{\varepsilon^{2}}{2}|\nabla\phi(\boldsymbol{x},t)|^{2}\right)d\Omega,$
(12)

where $\Omega$ is the bounded physical domain, $H^{1}(\Omega)$ denotes the
space of square-integrable real-valued functions with square-integrable
derivatives on $\Omega$, $L^{4}(\Omega)$ denotes the function space in which
the fourth power of the function is integrable, $\mathbb{R}_{\geq 0}$
represents the set of non-negative real numbers. The first term in the RHS is
the bulk or mixing energy and depends on the local composition of the mixture.
The functional form for the bulk energy $F(\phi)$ has been chosen as the
double-well potential function $F(\phi)=\frac{1}{4}(\phi^{2}-1)^{2}$ in the
present study. Thus, the value of the phase-field function within a pure solid
and fluid domain is given by $1$ and $-1$ respectively. The second term is the
interfacial or gradient energy and depends on the composition of the immediate
environment. The ratio of these two effects controlled by $\varepsilon$
decides the thickness of the diffused interface region. At equilibrium, the
interface thickness is the distance over which $\phi$ varies from $-0.9$ to
$0.9$ which can be estimated as $2\sqrt{2}\tanh^{-1}(0.9)\varepsilon\approx
4\varepsilon$. The final convective form of the Allen-Cahn equation after
addition of the Lagrange multiplier for mass conservation can be derived as
follows:

$\frac{\partial\phi}{\partial
t}+\boldsymbol{v}\cdot\nabla\phi=-\gamma\left(F^{\prime}(\phi)-\varepsilon^{2}\Delta\phi-\beta(t)\sqrt{F(\phi)}\right)\text{
on }\Omega\times[0,T],$ (13)

where $\gamma$ is the mobility parameter (of the order $\mathcal{O}(10^{-3})$
for all the test problems in this work) and $\beta(t)$ is the time-dependent
part of the Lagrange multiplier, given by
$\beta(t)=\frac{\int_{\Omega}F^{\prime}(\phi)d\Omega}{\int_{\Omega}\sqrt{F(\phi)}d\Omega}$.

For the variational form of the Allen-Cahn equation, we define
$\mathcal{S}^{h}$ to be the space of trial solutions, whose values satisfy the
Dirichlet boundary conditions and $\mathcal{V}^{h}$ to be the space of test
functions whose values vanish on the Dirichlet boundary. The final stabilized
form of the Allen-Cahn equation can be stated as: find
$\phi_{h}^{n+\alpha}\in\mathcal{S}^{h}$ such that $\forall
w_{h}\in\mathcal{V}^{h}$,

$\displaystyle\int_{\Omega}(w_{\mathrm{h}}\partial_{\mathrm{t}}\phi_{\mathrm{h}}^{\mathrm{n}+\alpha_{m}}+w_{\mathrm{h}}(\boldsymbol{v}^{\mathrm{n}+\alpha}\cdot\nabla\phi^{\mathrm{n}+\alpha}_{\mathrm{h}}$
$\displaystyle-\gamma(\nabla
w_{\mathrm{h}}\cdot(k\nabla\phi_{\mathrm{h}}^{\mathrm{n}+\alpha})-w_{\mathrm{h}}s\phi_{\mathrm{h}}^{\mathrm{n}+\alpha}+w_{\mathrm{h}}f)))\mathrm{d}\Omega$
(14)
$\displaystyle+\sum_{\mathrm{e}=1}^{\mathrm{n}_{\mathrm{el}}}\int_{\Omega_{e}}\left(\boldsymbol{v}^{\mathrm{n+\alpha}}\cdot\nabla
w_{\mathrm{h}}\right)\tau_{\phi}\mathcal{R}(\phi_{\mathrm{h}})\mathrm{d}\Omega_{e}=0,$

where $k$, $s$ and $f$ are the diffusion coefficient, reaction coefficient and
source terms respectively and $\mathcal{R}\left(\phi_{\mathrm{h}}\right)$ is
the element-wise residual of the Allen-Cahn equation as defined in [38]. The
first line of Eq. (14) contains the Galerkin terms, and the second line
contains the linear stabilization terms with the stabilization parameter,
$\tau_{\phi}$ [32, 33] given by

$\tau_{\phi}=\left[\left(\frac{2}{\Delta
t}\right)+\boldsymbol{v}\cdot\boldsymbol{G}\boldsymbol{v}+9k^{2}\boldsymbol{G}:\boldsymbol{G}+s^{2}\right]^{-1/2}.$
(15)

The boundary condition and initial condition for the Allen-Cahn equation are
as follows:

$\displaystyle\frac{\partial\phi}{\partial
n}\biggr{\rvert}_{\Gamma}=\boldsymbol{n}\cdot\nabla\phi=0\ \ \ \ $
$\displaystyle\text{on }\Gamma\times[0,T],$ (16)
$\displaystyle\phi|_{t=0}=\phi_{0}\ \ \ \ $ $\displaystyle\text{on }\Omega.$
(17)

We define the initial phase function $\phi_{0}$ based on the geometry of the
solid bodies. This completes the set of the governing equations for the
system. We now turn our attention to the central theme of the paper i.e. the
contact modeling in our phase-field based FSI framework.

##  3 Phase-field Contact Modeling

In this section, we discuss the collision modeling approach used in the
present work. We begin by describing the difficulties associated with the
conventional contact problem. Then we focus on the specific case of FSI
contact for diffuse interface bodies. Finally, we present the contact force
model that we use in the current work to resolve contact in FSI applications.

###  3.1 Contact Constraints

In numerical modeling of contact, there are two additional inequalities that
need to be satisfied besides the governing equations. The first inequality
derives from the fact that the deformation of the bodies should be such that
there is no inter-penetration. Hence the gap $g$ between the two bodies should
be non-negative,

$g\geq 0.$ (18)

Additionally, non-adhesive contact requires that there is no tension force in
the contact patch. This condition of compressive-only stresses translates to
the following inequality:

$\sigma_{nn}\leq 0,$ (19)

where the normal stress, $\sigma_{nn}=\boldsymbol{\sigma}:\mathbf{P}_{n}$, and
$\mathbf{P}_{n}=\mathbf{n}\otimes\mathbf{n}$ is the projection tensor.
Condensing these inequality constraints into an equality constraint, we arrive
at the following KKT (Karush-Kuhn-Tucker) form:

$g\geq 0,\text{ }\sigma_{nn}\leq 0,\text{ }\sigma_{nn}g=0.$ (20)

The above constraints are also known as the Hertz-Signiorini-Moreau (HSM)
conditions in the framework of contact problems. The equality condition states
that the part of the boundary in contact has $g=0$ and the part of the
boundary not in contact experiences zero stress, i.e. $\sigma_{nn}=0$. The gap
function for a general situation is given by,
$g=(\boldsymbol{x}_{1}-\boldsymbol{x}_{2})\cdot\boldsymbol{n}$, where
$\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$ are points on the two bodies and
$\boldsymbol{n}$ is the common unit normal. The exact formulation of the gap
depends on the underlying model used for the solid.

###  3.2 Contact Conditions for Submerged Elastic Bodies

Submerged bodies also experience the fluid pressure as an external traction.
To account for this effect of the presence of a thin fluid layer between the
colliding bodies, we modify the HSM conditions accordingly. Thus, ensuring a
continuous transition from the FSI traction to the contact traction, we arrive
at the following conditions:

$\displaystyle g$ $\displaystyle\geq 0,$ (21)
$\displaystyle\sigma^{s}_{nn}-\sigma^{f}_{nn}$ $\displaystyle\leq 0,$ (22)
$\displaystyle g(\sigma^{s}_{nn}-\sigma^{f}_{nn})$ $\displaystyle=0.$ (23)

The last line in the above set of constraints enforces that there is either a
vanishing gap as in solid-solid contact, or there is a vanishing relative
traction as in the case of FSI.

$g$$f_{c}$$g$$f_{c}$

Figure 1: Illustration of contact force vs gap for sharp (left) and diffused
(right) interface of solids.

The definitions for contact and penetration are altered slightly for solid
bodies with diffuse interface representations. In a diffused interface
setting, contact is defined by the overlap of the diffused region of one
object and the boundary of the other. And penetration is said to occur when
the zero iso-contours of the order parameters of the two bodies cross each
other. Thus, the presence of a diffuse interface ameliorates the issue of
contact force discontinuity by allowing for a more gradual contact force (Fig.
1). A non-penetrating contact routine allows the zero iso-contours to be
separated by less than half of the transition width. When the solids are
further apart than this distance, the routine has no effect.

###  3.3 Contact Force Computation

$\Omega_{s}^{1}$$\Omega_{s}^{2}$$\Omega_{f}$$\Gamma_{FSI}^{1}$$\Gamma_{FSI}^{2}$$\boldsymbol{f}_{c}^{2}$$\boldsymbol{f}_{c}^{1}$$e_{1}$$e_{2}$$e_{3}$

Figure 2: Schematic for the contact force between two colliding bodies. Two
solid bodies $\Omega_{s}^{1}$ (black) and $\Omega_{s}^{2}$ (blue) with diffuse
interfaces and immersed in a surrounding fluid ($\Omega_{f}$), approach
contact as shown above. The hatched regions denote the area where the contact
force is applied and the arrows represent the direction of the total contact
force acting on the two bodies.

In the current approach, we define an explicit force similar to [23, 24, 22]
as a consequence of two bodies coming into contact. Any further approach of
the two solids beyond the transition zone is penalized by the force function.
This acts as a substitute of the reaction force that would result from a
collision. This force is a function of the amount of overlap between the two
diffuse solids and the physical properties of the solids. The formulation is
closest to the conventional penalty method for contact handling in terms of
implementation. However, unlike the penalty method, the scope of inter-
penetration between the bodies is very less in the current approach. The
method relies on first identifying the boundaries of the two solids at each
time step and evaluating the minimum gap between them, using a node-node
search. If the gap is below half of the diffuse interface thickness ($g\leq
2\varepsilon$), we call the contact force subroutine to identify the amount of
overlap, and apply the appropriate force. At each time step, we compute the
quantity $d_{12}$ [34], which is defined as:

$\displaystyle d_{12}=\frac{d_{1}-d_{2}}{2},$ (24)

where $d_{1}$ and $d_{2}$ are the minimum distances of the point from the zero
iso-contours of the two bodies. These can be calculated once we have the
boundary nodes on both solids. We note that $d_{12}=0$ denotes a point
equidistant from the two bodies. We can use this quantity to construct a
short-range forcing function and define a common unit normal. The spring-type
influence function for our case is defined as:

$\displaystyle\delta(x)=\left(1-\frac{x}{2\varepsilon}\right),$ (25)

where $4\varepsilon$ is the thickness of the diffuse interface. The force
function in the present work is defined as:

$\displaystyle\boldsymbol{f}_{c}(d_{12})=\kappa\frac{\mu_{s}^{L}}{h}\delta(d_{12})\boldsymbol{n},$
(26)

where $\kappa$ is a constant used to control the magnitude of the force. The
computed contact force is then incorporated in the body force term in the
continuum equation (Eq. 2). In problems involving gravity, the total body
force is given as:

$\boldsymbol{b}(\phi)=\rho(\phi)\boldsymbol{g}+\boldsymbol{f}_{c}.$ (27)

The dependence of the contact force on the shear modulus and the mesh
resolution has been shown in the literature [39]. $\boldsymbol{n}$ is the
common unit normal defined as $\boldsymbol{n}=\pm\frac{\nabla d_{12}}{|\nabla
d_{12}|}$. We also account for the collision of an object with the domain
boundaries (rigid walls), where the distance is computed between the $\phi=0$
boundary of the body to the rigid wall and the unit normal of the fixed wall
is utilized for the force direction.

It is worth mentioning that the method described above is an explicit
interface extraction method as we need to recover the diffused interface and
the $\phi=0$ boundary at every time step to compute the distances and the
amount of overlap. Additionally, when computing the contact forces, we scan
through the domain again to identify the nodes present in the interior of the
solid ($\phi>0$) to apply the forces appropriately. It is also restricted to
two body systems only. We present an improved contact handling algorithm for
multibody systems in the next section that circumvents the aforementioned
issues.

##  4 Mono-field Interface Advancing for Collision Detection

$\Omega_{f}$$\Omega_{s}^{1}$$\Omega_{s}^{2}$$\Omega_{s}^{3}$$\Omega_{s}^{4}$$\Omega_{s}^{5}$$\mathcal{B}_{2}^{n}$$\mathcal{B}_{2}^{n+1}$

Figure 3: Illustration of mono-field interface advancing method for collision
detection. The above schematic shows five different solid bodies ($N=5$)
interacting and colliding in a fluid. The lists at the bottom denote the
sorted rank lists for each object, where the number $i$ stands for body
$\Omega_{s}^{i}$. All the orange-green pairs of objects are passed to the
contact force subroutine as the distance between them ($g_{ij}$) is less than
the specified threshold, whereas the orange-gray pairs are skipped while
traversing the lists. The inset at the top shows the boundary of the second
body $\mathcal{B}_{2}$ at two different time steps translating and deforming
in the fluid domain. The black lines denote the solid at $t=n$
($\mathcal{B}_{2}^{n}$), whereas the red lines denote the body at $t=n+1$
($\mathcal{B}_{2}^{n+1}$). We have highlighted some representative nodes based
on their status before and after the update step. The blue circles represent
the nodes that are present both in $\mathcal{B}_{2}^{n}$ and
$\mathcal{B}_{2}^{n+1}$. The black circles denote the nodes that were present
in $\mathcal{B}_{2}^{n}$ but have been removed in $\mathcal{B}_{2}^{n+1}$.
Finally, the red circles denote the nodes that were absent in
$\mathcal{B}_{2}^{n}$ and have been added newly in $\mathcal{B}_{2}^{n+1}$.

The presented contact force routine can be generalized to multi-body
simulations as well, especially for cases with identical physical properties,
for example, ice floes on the water surface. It is important to distinguish
between the different solid bodies present in the domain to identify the pairs
of bodies that undergo contact at a particular instant. However, the explicit
interface extraction method presented previously is insufficient for problems
involving more than two bodies. Thus, we introduce the mono-field interface
advancing (MFIA) method to keep track of the individual solid interfaces as
they evolve in time, while using one phase-field equation. The traditional
approach for such a problem would be to first assign individual objects with
their own order parameter field and strain/ deformation field and solve these
N equations for the N solid bodies at every time step. However, this can prove
to be computationally expensive with increase in the number of objects. The
proposed procedure is still a pairwise contact model but now involves updating
the boundary node lists for each object instantiated at $t=0$, instead of
reconstructing the interface at every time step. This allows us to get away
with one interface capturing field and one equation for the strain evolution
for the N solid bodies in the physical domain.

Let us consider that we have N objects that make up the solid domain
$\Omega_{s}$ and one fluid phase $\Omega_{f}$ for simplification. We
initialize the phase-field function of each solid body based on its location
in the domain and isolate its boundary nodes $(\mathcal{B}_{i})$ in the mesh
at $t=0$ based on their $\phi$ values such that
$-0.9\leq\phi(\mathcal{B}_{i})\leq 0.9$. The unified phase-field function for
the entire domain is just a scalar sum of these individual functions and off-
setting the sum to constrain between $[-1,1]$. The first stage of the
procedure is the updation of the boundary node lists at each time step. At
every time $t=t$ and for each node list $\mathcal{B}_{i}$, we check for nodes
that no longer belong in the interface zone and remove them from the list. We
then perform a domain scan for the phase indicator values of the nodes and
collect the set of nodes $\Gamma$ that belong to the interface zone at this
step ($-0.9\leq\phi(\Gamma)\leq 0.9$). Using this, we create a new set $S$
such that $S=\Gamma\cap(\cup_{i}\mathcal{B}_{i})$. This step isolates all the
nodes that have recently entered the interface zone in the entire domain. The
next step involves assigning the nodes in the set $S$ to the correct
$\mathcal{B}_{i}$. We do this by computing the distances of the particular
node to the different bodies and finding the minimum among them. This
completes the boundary update step of the procedure and can be visualized in
the inset for Fig. 3.

The next stage deals with creating the rank list for each solid body. For
every solid body with boundary nodes $\mathcal{B}_{i}$, we compute its minimum
gap from other bodies and arrange them in the ascending order of distances
$g_{ij}$, as illustrated in Fig. 3. The rank list for a higher-numbered body
does not contain smaller-numbered bodies as those pairs are already accounted
for in the earlier lists. Finally, we pass the pair of bodies that satisfy the
minimum distance criterion to the contact force sub-routine (orange-green
pairs in Fig. 3). Please note that we have excluded body-wall forces in this
explanation. However, the actual implementation contains the rigid walls in
the rank lists as well. The computed forces from each pair are then included
in the unified momentum balance equation. For example, with reference to Fig.
3, forces between pairs 1-5 and 1-2 are computed separately and included in
the one-field momentum equation. Also note that in the steps involving
distance computations in the above procedure, we trim the boundary nodes to a
tighter threshold such as $-0.1\leq\phi\leq 0.1$. This improves the accuracy
of the distance calculations from the $\phi=0$ boundary. The aforementioned
procedure is summarized in Algorithm 1.

Algorithm 1 Mono-field interface advancing (MFIA) for collision detection

1: Initialize phase functions for N bodies and create boundary lists
$\mathcal{B}_{i}$ for every solid body

2: Loop over time steps $\mathrm{n}=1,2,\cdots$

3: Loop over each boundary list $B_{i}$

4: Remove nodes that do not lie in the interface zone: $-0.9\leq\phi\leq 0.9$

5: Loop over all nodes and find the ones that are in the interface zone at
$t=t_{n}$

6: Extract the set of nodes that have recently entered the interface zone

7: Assign these nodes to the correct boundary list $\mathcal{B}_{i}$

8: Loop over the solid bodies

9: Create the rank list based on the proximities from other bodies

10: Obtain the pairs of potentially colliding bodies based on the minimum gap
criterion

11: Pass these pairs to the contact force sub-routine (Eq. 26)

The novel MFIA approach allows us to simulate multiple bodies of identical
physical parameters immersed in a surrounding fluid. It eliminates the need
for solving N number of level set/ phase field equations and solid strain
equations for each of the N bodies. This also prevents creating new lists at
each time step by tweaking the existing lists by a small amount. The above
approach can be considered similar to the fast marching method used in level
set problems [40]. The fast marching method involves solving the convection
equation for the level set in a narrow band around the existing front in an
upwind fashion (or normal to the front) and maintaining the values for the
far-away points. However, in the current approach, we solve the unified phase-
field equation in the full domain to maintain the thermodynamic consistency
and only carry out the above procedure to distinguish the different bodies.
The assumptions involved are that the mesh size and time step size are
sufficiently small to avoid excessive displacements, which are reasonable for
any contact simulation. It is also restricted to convex bodies. Again, this
assumption is reasonable as any arbitrary object can be considered as
comprising several smaller convex bodies.

######  Remark 1.

Another related approach is the interface reconstruction technique in the cut-
cell method or in X-FEM which is carried out to maintain the sharp interface
nature of the method [41, 42]. We would like to highlight, however, that in
the present approach, we are not concerned with the location of the exact
sharp interface. We are only interested in separating the updated diffuse
regions of each body after we solve the Allen-Cahn equation. Thus, the above
procedure does not pertain to any form of interface or geometry-preserving
ideas [29] but can be used in conjunction with one.

######  Remark 2.

The proposed method avoids the use of any material/ Lagrangian points as done
in particle level set [43] or additional indicator functions as done for front
tracking methods [44] and relies entirely on the Eulerian mesh nodes.

This completes the description of our Eulerian phase-field based contact
formulation. We will now present some implementation details of the
variational finite element solver before turning our attention to verification
and applications.

##  5 Implementation Details

In this section, we present the details of the numerical implementation of the
proposed framework. We first initialize the velocities, pressure, order
parameters and deformation tensors in the domain and apply the boundary
conditions depending on the problem at hand. We initialize the boundary lists
$\mathcal{B}_{i}$ for each body $i$ at $t=0$. The boundary nodes on each wall
of the domain are stored in separate arrays in addition to the above lists and
the wall normals are computed and saved. The solver then enters the time loop
where the nonlinear set of equations are solved in a partitioned manner [28].
Inside the time loop, we first check for potentially colliding bodies and
impose the appropriate contact forces before entering into the non-linear
Newton iterations. Implicit discretizations are carried out for the underlying
equations for improved efficiency and robustness of the solver. The
incremental quantities are computed by employing Newton-Raphson type
iterations and the generalized-$\alpha$ time integration is used to update the
solver. The Jacobians are calculated and the matrices are assembled for each
of the equations separately and then solved sequentially.

The updated velocities from the solution of the unified continuum equations
are used for evolving the solid strains via the $\boldsymbol{B}$ transport
equation and updating the order parameter values through the Allen-Cahn
equation. Solid stresses are computed from the new deformation tensor
$(\boldsymbol{B})$ values. The updated order parameter values are further
utilized in the interpolation of density $\rho(\phi)$, viscosity $\mu(\phi)$
and stress $\boldsymbol{\sigma}(\phi)$. The updated phase functions also help
in updating the solid boundaries and estimating the gap between the solid
bodies during the broad-phase collision detection step. The pairwise contact
forces obtained for the various elastic-elastic and elastic-rigid collisions
are computed based on the amount of overlap of the diffuse interfaces. We
evaluate the total contact force acting on a given body as a vector sum of all
the collision forces acting on it. These additional forces that contribute to
the body force along with the other updated physical quantities are fed to the
unified continuum equation for the next iteration.

##  6 Verification and Convergence

We provide a systematic verification of the above contact force implementation
for a dry contact scenario and an FSI contact scenario in this section. We
utilize the analytical Hertzian contact solutions for the dry contact or
solid-only contact verification and the problem of two-body collision in a
Taylor-Green vortex field for the FSI contact verification. We also present
the convergence of the implemented contact routine with respect to the mesh
resolution for a fixed interface resolution ($\frac{\varepsilon}{h}$) for the
Hertz verification case.

###  6.1 Dry Contact Verification

The analytical solutions available for smooth Hertzian contact between two
linear elastic bodies are utilized for the assessment of our phase-field model
in dry contact scenarios.

####  6.1.1 Hertz Analytical Solution

Analytical solutions exist for most of the common Hertzian contact problems.
They exploit the idea of superposition for linear elastic problems and solve
the governing equations in different directions independently to compute the
stresses. The strains and displacements are evaluated consequently. Our focus
in the current work is on 2D contact verification, and thus we consider the
cylindrical Hertzian contact here. This problem involves the normal contact of
two right circular cylinders with their axes parallel to each other. Let the
cylinders have radii $R_{1}$ and $R_{2}$, Young’s moduli $E_{1}$ and $E_{2}$
and Poisson’s ratios $\nu_{1}$ and $\nu_{2}$ respectively. We assume the
length of the cylinders to be unity for convenience. Their contact results in
an elliptic stress distribution of width $2a$ across the contact zone. The
half-width $a$ for this scenario is given by:

$a=\left(\frac{4PR_{eq}}{\pi E_{eq}}\right)^{\frac{1}{2}},$ (28)

where $P$ is the total contact force which is equal and opposite for both
bodies, $R_{eq}$ is the equivalent radius of curvature,
$R_{eq}=\frac{R_{1}R_{2}}{R_{1}+R_{2}}$ and $E_{eq}$ is the equivalent Young’s
modulus,
$\frac{1}{E_{eq}}=\frac{(1-\nu_{1}^{2})}{E_{1}}+\frac{(1-\nu_{2}^{2})}{E_{2}}$.
The maximum contact pressure is given by,

$p_{0}=\frac{2P}{\pi a}.$ (29)

The elliptic traction profile is given by,

$\frac{p(r)}{p_{0}}=\left(1-(r/a)^{2}\right)^{\frac{1}{2}}\text{ in }0\leq
r\leq a.$ (30)

The above relations apply to the contact of a cylinder with a plane surface as
well, where $R=\infty$ for the plane surface.

####  6.1.2 Deformable Cylinder on Elastic Plane Surface

$\Omega_{s}^{1}$$\Omega_{s}^{2}$$R$$v_{x}=-u_{0}$

Figure 4: Schematic of a simple Hertzian contact problem with a deformable
cylinder contacting an elastic plane. The plane wall on the right approaches
the fixed cylinder on the left with a constant velocity boundary condition at
the right boundary.

(a)

(b)

(c)

(d)

Figure 5: Hertz contact problem on an elastic plane: (a) mesh for the contact
problem between a deformable cylinder on an elastic plane. The mesh is refined
within the expected contact region. (b) Pressure contour for the contact
problem, (c) normalized traction profile (for $\varepsilon=0.03$) comparison
with analytical Hertz solution and (d) convergence of the relative L2 error in
the traction profile with respect to $h$. The numerical traction profile
plotted in (c) has been smoothed out by carrying out a moving average of the
obtained data. The interface resolution is similar for all the cases,
$\frac{\varepsilon}{h}\approx 10$.

In this test case, we consider a deformable cylinder with $R=1$ that is
compressed by an elastic plane. We model the second phase as air to emulate
dry solid-solid contact. We use the incompressible Neo-Hookean model to
represent both the solids and assign a Young’s modulus of 1. Other physical
properties are: $\rho_{f}=10^{-3}$, $\rho_{s}=1$, $\mu_{f}=10^{-6}$. All the
physical parameters are non-dimensional in nature. The center of the cylinder
is located at (0,0). The initial condition for the order parameter is defined
as:

$\phi(x,y,0)=1+\mathrm{tanh}\left(\frac{R-\sqrt{(x-0)^{2}+(y-0)^{2}}}{\sqrt{2}\varepsilon}\right)+\mathrm{tanh}\left(\frac{x-1.1}{\sqrt{2}\varepsilon}\right).$
(31)

We initialize the elastic plane in the region $1.1\leq x\leq 2$ so that it
slowly approaches the cylinder and presses on it. The traction profiles are
extracted when the total contact force ($P$) on the cylinder reaches a
prescribed value ($\approx 0.005$ in our case). The cylinder is fixed on the
left boundary. A constant velocity of $R/100$ is applied at the right boundary
to emulate the slow approach of the plane toward the cylinder. During post-
processing, we compute the traction profile by integrating the contact body
force $\boldsymbol{f}_{c}$ normal to the interface. The total contact force is
given by the volume integral,
$P=\int_{\Omega_{cyl}}\boldsymbol{f}_{c}\cdot\boldsymbol{n}\mathrm{d}\Omega$
for the cylinder. Numerically, the traction computation is done as follows. We
first define a new variable that stores the quantity
$f=\boldsymbol{f}_{c}\cdot\boldsymbol{n}$. We define a local rectangular grid
around the contact patch to approximate the diffused contact region. We
interpolate the values of the newly computed quantity over the grid and carry
out quadrature integration over $x$-points for each $y$ location as
$p(y)=\int_{x_{1}}^{x_{2}}f(x,y)dx$. This gives us the traction profile $p(y)$
along the contact patch.

Figure 5(c) shows the comparison of the obtained traction profile with the
analytical Hertz profile. Overall, the numerical traction profile follows the
elliptic trend of the analytical solution. However, some oscillations can be
observed in the solution which can be minimized by using better reconstruction
approaches while computing the contact pressure. Some of the other possible
causes for the deviation might be (i) the use of a diffused interface profile
which results in a volumetric contact force (area force in 2D) instead of a
surface force (line force in 2D); (ii) an incompressible hyperelastic material
law compared to a compressible linear elastic material in Hertz theory; and
(iii) the finite size of the computational domain as compared to the semi-
infinite assumption of Hertz. We have plotted the convergence of the relative
$L_{2}$ error in the traction profile with different mesh resolutions in the
contact zone $h$ in Fig. 5(d). We define the relative $L_{2}$ error in the
obtained traction profile as:

$||e||_{2}=\frac{||p(y)-p_{Hertz}||_{2}}{||p_{Hertz}||_{2}},$ (32)

where $p_{Hertz}$ is the traction profile obtained via Hertz theory using the
numerically computed total contact force $P$. We can observe that as we reduce
the mesh size, the error in the traction profile decreases monotonically. The
convergence rate is found to be between first and second order ($\approx
1.8$).

###  6.2 Elastic-Elastic Contact in Fluid Environment

To ensure the accuracy of the implemented model for generalized FSI contact
cases, we present a simple verification of two hyperelastic bodies colliding
in the Taylor-Green vortex field. We compare the trajectories of the centroids
of the two bodies with other works in the literature and present discussions
on the energy transfer mechanisms and volume conservation for this case.
Finally, we also illustrate the cost-effectiveness of the proposed MFIA
approach over a simplistic interface extraction method for collision detection
and contact handling using this test case.

####  6.2.1 Two Body Collision in Taylor-Green Vortex Field

(a)

(b)

(c)

(d)

$\omega_{z}:$

Figure 6: Two elastic body collision problem in a Taylor-Green voretx field.
Snapshots of the vorticity contours for the two body collision problem at $t=$
(a) $0$, (b) $0.3$, (c) $0.6$ and (d) $0.9$. The black solid lines denote the
$\phi=0$ iso-contours of the two bodies.

(a)

(b)

Figure 7: Two elastic body collision problem in a Taylor-Green voretx field.
Time evolution of the centroids of the top and bottom bodies: (a) comparison
for different mesh resolutions ($h$) and (b) comparison with previous works by
Bhosale et al. (2021) [45] and Jain et al. (2019) [46] for the $h=0.005$ case.

To simulate this problem, we initialize two incompressible hyperelastic bodies
in a Taylor-Green vortex field and allow them to collide and separate
eventually. This problem has been carried out previously in [46, 45] to
simulate soft solid contact. The non-dimensional geometric and physical
parameters for the problem are as follows. We define a unit square over
$[0,1]\times[0,1]$ as the computational domain for this simulation. Two
incompressible solids that are initially circular with radii $R_{1}=R_{2}=1/6$
are placed at $(0.5,0.7)$ and $(0.5,0.3)$. The phase-field function for the
two body system is initialized as follows:

$\displaystyle\phi(x,y,0)=1$
$\displaystyle+\mathrm{tanh}\left(\frac{R_{1}-\sqrt{(x-0.5)^{2}+(y-0.7)^{2}}}{\sqrt{2}\varepsilon}\right)$
(33)
$\displaystyle+\mathrm{tanh}\left(\frac{R_{2}-\sqrt{(x-0.5)^{2}+(y-0.3)^{2}}}{\sqrt{2}\varepsilon}\right).$

The initial Taylor-Green vortex field given by the streamfunction
$\psi=\psi_{0}\mathrm{sin}(k_{x}x)\mathrm{sin}(k_{y}y)$ is imposed using the
following velocity distribution in the domain:

$\displaystyle v_{x}(x,y)$
$\displaystyle=\psi_{0}k_{y}\mathrm{sin}(k_{x}x)\mathrm{cos}(k_{y}y)$ (34)
$\displaystyle v_{y}(x,y)$
$\displaystyle=-\psi_{0}k_{x}\mathrm{cos}(k_{x}x)\mathrm{sin}(k_{y}y),$ (35)

where $\psi_{0}=\frac{1}{2\pi}$, $k_{x}=2\pi$ and $k_{y}=2\pi$ for the present
case. Other physical parameters have been chosen the same as those in previous
works ([45]) with fluid viscosity $\mu_{f}=\frac{0.01}{2\pi}$, shear modulus
$\mu_{s}^{L}=2$ and density $\rho_{f}=\rho_{s}=1$. The interface resolution is
considered to be $\frac{\varepsilon}{h}=1$ and the time step size is 0.001 for
all simulations of this test case. We impose the periodic boundary conditions
on all sides for the momentum equation and homogeneous Neumann conditions for
the phase-field equation. Figure 6 shows the temporal evolution of the
centroids of the two solids as they approach collision and separate
eventually.

####  6.2.2 Assessment of Volume Conservation and Energy Transfer

(a)

(b)

Figure 8: Comparison of temporal evolution of (a) relative volume error for
the top body and (b) kinetic and strain energies of the system for different
mesh resolutions. The dashed black line in the left figure shows the ideal
case of constant volume.

We measure the error in volume conservation of the solids during the collision
and separation processes for the present case. For an incompressible solid,
this quantity gives us an estimate of the deviation of
$J=\mathrm{det}(\boldsymbol{F})$ from unity. We define the relative volume
error for the top body using the following definition:

$e_{vol}=\frac{\int_{\Omega_{s,1}^{t}}\mathrm{d}\Omega-\int_{\Omega_{s,1}^{0}}\mathrm{d}\Omega}{\int_{\Omega_{s,1}^{0}}\mathrm{d}\Omega}.$
(36)

We plot this quantity for different mesh resolutions in Fig. 8(a). It can be
concluded that the magnitude of the relative volume error is within $5\%$ for
all the cases and within $2\%$ for the finer resolution cases ($h\leq 0.01$).
The deviation from the ideal scenario of constant volume can be attributed to
a lower interface resolution ($\frac{\varepsilon}{h}=1$) of the diffused
interface body. This issue can be alleviated by the application of adaptive
mesh refinement techniques as shown in [28].

We also track the evolution of energy for this test case to illustrate the
transfer of kinetic energy to solid strain energy and vice versa (Fig. 8(b)).
Kinetic energy of the system is defined as:

$KE=\int_{\Omega}\frac{1}{2}\rho||\boldsymbol{v}||^{2}\mathrm{d}\Omega$ (37)

and strain energy for the incompressible neo-Hooekan solids is given by

$SE=\int_{\Omega_{s}}\frac{1}{2}\mu_{s}^{L}(tr(\boldsymbol{B})-d)\mathrm{d}\Omega$
(38)

where $d$ is the dimension of the problem i.e. $2$ in our case. We can observe
that the kinetic energy is maximum at the initial time step due to the imposed
velocity conditions and dissipates gradually because of the viscous fluid.
However, the low-frequency oscillations during this dissipation process are
caused by the conversion of the kinetic energy into the solid strain energy
and vice versa. The kinetic energy reaches its minimum when the solids are
closest to each other. Concurrently, the solid strain energy peaks since the
solids have achieved their maximum deformation. Eventually, the total energy
decays owing to the viscosity of the fluid.

####  6.2.3 Efficiency of MFIA

(a)

(b)

Figure 9: Comparison of (a) the non-dimensional energies and (b) the total
elapsed time for the two contact handling approaches for the two body
collision problem. The time axis (x-axis) has been normalized with the total
simulation time.

In this section, we contrast the MFIA with the interface extraction process
for detecting collision and resolving contact. We utilize the two-body
collision problem here as it is not feasible to use the interface extraction
approach for multibody contact problems. Figure 9(a) illustrates the
comparison of the energy transfer curves between the two algorithms to verify
the results. We plot the total elapsed time for this test case using both
approaches in Fig. 9(b) to illustrate the cost effectiveness of the MFIA. We
observe a $10\%$ reduction in the total elapsed time for the MFIA as compared
to the interface extraction approach. This can be attributed to two important
factors. First, the interface extraction approach discards the old interfaces
at every timestep and constructs new ones by comparing the radius vectors and
outward normals. Second, when imposing the contact force, we do not have to
scan through the entire domain again to identify eligible solid nodes that
satisfy the minimum gap threshold. Instead, we re-use the saved interface
arrays and scan through them to identify the eligible nodes. Hence, the
difference in the elapsed time is expected to rise as we increase the domain
size and/or reduce the volume occupied by the solid bodies in the
computational domain.

##  7 Applications

In this section, we present two applications to illustrate the robustness of
the proposed phase-field based contact model. The first problem of a two-body
collision system highlights an actuation mechanism for elastic bodies in fluid
flow. The second test case demonstrates the capability of the solver to handle
multibody contact problems involving complex hydrodynamic interactions.

###  7.1 Two Body Collision with Actuation

Our approach also admits a simple method for solid actuation. This is
especially advantageous when simulating clamped solids or prescribing a
particular motion to the solid with applications ranging from active soft
matter [47] to bio-locomotion studies [48, 49, 45].

####  7.1.1 Actuation of Bodies in Fluid Flow

We utilize the constitutive model of the solid and the imposed kinematic
conditions to achieve our objective. To actuate a particular solid, a smaller
anchored region is defined inside the larger solid object. The strain tensor
(i.e. the left Cauchy-Green deformation tensor $\boldsymbol{B}$ in the present
case) is then set to be equal to the identity tensor $\boldsymbol{I}$ within
the anchor zone. This restrains any deformation of the solid inside the anchor
during the simulation. The next step is to specify the desired kinematic
conditions in terms of imposed velocity in the anchored region. The rest of
the body deforms in response to the imposed motion in the anchor. During the
numerical implementation, we identify the nodes in the anchor as fixed nodes
and remove them from the linear solve process. Another approach to achieve
solid actuation can be to specify the desired deformation within the anchored
zone as demonstrated in [34]. However, it is cumbersome to use this approach
for cases where we do not know the exact deformation of the body or for cases
where we would ideally want to impose a particular kinematic constraint.

$\Omega_{s}^{1}$$\Omega_{s}^{2}$$\Omega_{f}$(fixed)(moving)

Figure 10: Schematic for the two body collision problem. The anchored regions
have been shaded for both the bodies. The anchor on the left body is fixed
whereas a harmonic x-velocity is provided to the anchor on the right body.

(a)

(b)

(c)

(d)

Figure 11: Two body collision problem: Snapshots of the pressure contours at
non-dimensional times $t=$ (a) $0$, (b) $0.36$, (c) $0.5$ and (d) $0.72$. The
physical time has been non-dimensionalized with the time period of the imposed
motion. The left ball is anchored at $(-1.1,0)$ and the right ball is given a
prescribed velocity via an anchored region initially centred at $(1.1,0)$. The
black solid lines denote the $\phi=0$ iso-contours of the two bodies.

We present a test case to demonstrate the effectiveness of the proposed
actuation approach. To illustrate the robustness of the unified approach, two
incompressible hyperelastic solid bodies are considered that are actuated
differently and undergo contact. This problem has also been carried out by
[23] in the context of an imposed reference/inverse map condition within a
level set based finite difference discretization. The domain is a square of
side length 4 on $[-2,2]\times[-2,2]$ and the radius of the circular discs are
$R_{1}=R_{2}=0.7$ initially. The physical parameters for the problem are
densities, $\rho_{f}=\rho_{s}=1$, fluid viscosity, $\mu_{f}=0.1$, shear
modulus, $\mu_{s}^{L}=10$ and solid viscosity, $\mu_{s}=0.007$. Both bodies
are provided with an anchored region of radius 0.25 at their center which
controls their motion. The left body is fixed at its initial location while
the right body is given a harmonic motion of the following form:

$\displaystyle v_{x}(t)=-\frac{\pi}{25}\sin\left(\frac{2\pi t}{25}\right).$
(39)

We simulate this test case on a $256\times 256$ grid with a time step of
$0.005$. We normalize the physical time with the time period of the above
imposed motion, i.e. $T=25$. The initial condition for the order parameter in
this problem is considered as follows:

$\displaystyle\phi(x,y,0)=1$
$\displaystyle+\mathrm{tanh}\left(\frac{R_{1}-\sqrt{(x+1.1)^{2}+(y-0)^{2}}}{\sqrt{2}\varepsilon}\right)$
(40)
$\displaystyle+\mathrm{tanh}\left(\frac{R_{2}-\sqrt{(x-1.1)^{2}+(y-0)^{2}}}{\sqrt{2}\varepsilon}\right).$

The interface resolution is assumed to be $\frac{\varepsilon}{h}=2$.

Figure 11 illustrates the pressure contours obtained from the above
implementation for two colliding bodies immersed in a fluid. The results match
qualitatively with those presented in [23]. The snapshots confirm the accurate
motion of the two bodies via the active solid actuation approach. We can
observe that the pressure builds up between the two bodies as they approach
each other. The right body reaches its extremum at a non-dimensional time of
$t=0.5$ and separates from the left body gradually due to the application of
the contact forces. We are able to avoid any numerical sticking between the
objects via the implemented contact routine. A negative pressure region is
created as it moves away and returns to its initial configuration.

###  7.2 Free-falling of Elastic Bodies in Fluid

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12: Free-falling of elastic bodies: Snapshots of the pressure contours
for the problem at $t=$ (a) $0$, (b) $1$, (c) $2$, (d) $3$, (c) $4$ and (f)
$5$. The six bodies start at the top of the domain and fall freely under the
action of gravity, undergoing complex hydrodynamic interactions and collision.
The black solid lines denote the $\phi=0$ iso-contours of the bodies.

In this section, we present a test case for demonstrating the robustness of
the implemented multibody contact routine. We choose a sedimentation problem
of six elastic bodies in an open tank [50]. The problem setup is similar to
studies of suspension and sedimentation of particles carried out in [51, 52,
53] with applications in fluidized beds, sediment transport in rivers and
blood flow in arteries. The elastic solids are initialized at the top of the
tank and allowed to fall freely under the action of gravity. The problem
involves complex hydrodynamic interactions, elastic-elastic collisions between
the solid bodies and elastic-rigid collisions between the bodies and the rigid
walls. The domain is a rectangle on $[-3,3]\times[-7,1]$ and the radius of the
circular discs are 0.5 initially. The physical properties for the problem are
as follows: $\rho_{f}=1$, $\rho_{s}=5$, $\mu_{f}=0.001$ and
$\mu_{s}^{L}=10^{4}$. A gravitational acceleration of $2$ units is applied in
the downward direction. All physical parameters are non-dimensional in nature.
The initial condition for the phase-field function is defined as:

$\phi(x,y,0)=5+\sum_{i=1}^{6}\mathrm{tanh}\left(\frac{R_{i}-\sqrt{(x-x_{C}^{i})^{2}+(y-y_{C}^{i})^{2}}}{\sqrt{2}\varepsilon}\right),$
(41)

where $(x_{C}^{i},y_{C}^{i})$ is the coordinate of the centroid and $R_{i}$ is
the radius of the $i^{th}$ body. The interface resolution is considered to be
$\frac{\varepsilon}{h}=2$ and a time step size of 0.001 is chosen. The no-slip
boundary conditions for velocity are imposed on the left, right and bottom
walls with the do-nothing boundary condition at the top to simulate an open
tank. Homogeneous Neumann conditions are applied on all the boundaries for the
Allen-Cahn equation.

Figure 12 illustrates the pressure contours for the above test case at six
different time instants. The local pressure increases as the bodies approach
the rigid walls or collide against each other. The applied contact forces
ensure non-penetration and allow the bodies to separate after colliding. After
a long time, all the bodies settle at the bottom of the tank.

##  8 Conclusions

In the current work, we presented an efficient way of modeling multibody
contact problems in a phase-field-based fully Eulerian framework. Our finite-
element Eulerian framework helps us to model large deformations and topology
changes such as contact mechanics in a comparatively simple manner on a fixed
mesh. The phase-field method allows for a gradually increasing contact force
implementation via its diffused interface representation. We verified our
contact force model with the traditional Hertz contact problem for smooth, dry
contact between a cylinder and a plane and obtained the expected elliptical
contact force profile. The spatial convergence rate was around $1.8$. We also
assessed our model for a more general FSI contact scenario of two hyperelastic
bodies colliding in a Taylor-Green vortex field against available literature.
We showcased the volume conservation properties and energy transfer mechanisms
for this problem via our Eulerian phase-field model. The relative volume error
was shown to be within $5\%$ for all the cases simulated with this problem. We
presented an approach to model flexible actuated solids in submerged
environments and illustrated its utility in a two-body collision problem. The
proposed MFIA method avoids the use of redundant interface capturing and
strain evolution equations for multiple solids with identical physical
properties. It also precludes the necessity of extracting the solid boundaries
at every time step in case of diffused interface descriptions such as phase-
field. These factors helped in reducing the computational time by around
$10\%$ compared to the explicit interface extraction method. Using our
proposed phase-field Eulerian framework, we demonstrated two test cases namely
contact between two submerged bodies with actuation and the free-falling of
multiple elastic bodies in an open tank.

In the current work, we have restricted ourselves to 2D problems for
verification and demonstration purposes. The next step involves extending the
present algorithm to 3D and integrating it with a multi order parameter solver
to simulate ship-ice-water-air systems. This will allow us to capture the rich
multiphysics phenomena involving complex hydrodynamic and solid-solid
interactions observed in such an application.

##  Appendix A Jacobian Computation for the Solid Stress Terms

In order to obtain the linearized form of the solid strain evolution, we
utilize the generalized alpha time integration scheme. As mentioned in Section
2.3, we decompose the solid stress into its volumetric and deviatoric
components to have a uniform description for pressure. Hence, we will focus on
deriving the Jacobian for the deviatoric component of a 2D stress state in
this appendix. We expand the deviatoric part of the solid stresses in terms of
the left Cauchy-Green tensor, utilizing the incompressible neo-Hookean
material model:

$\boldsymbol{\sigma}_{s}^{\prime}=\mu^{L}_{s}(\boldsymbol{B}^{\mathrm{n}+\alpha}-\frac{1}{2}tr(\boldsymbol{B}^{\mathrm{n}+\alpha})\boldsymbol{I}).$
(1)

We need to evaluate $\boldsymbol{B}^{\mathrm{n}+\alpha}$ to substitute in the
above equation to calculate the stresses. For this, we make use of the
generalized-$\alpha$ time integration and the evolution equation for left
Cauchy-Green tensor. From the generalized-$\alpha$ time integration, we have

$\displaystyle\boldsymbol{B}^{\mathrm{n}+\alpha}$
$\displaystyle=\boldsymbol{B}^{\mathrm{n}}+\alpha(\boldsymbol{B}^{\mathrm{n}+1}-\boldsymbol{B}^{\mathrm{n}}),$
(2a) $\displaystyle=\boldsymbol{B}^{\mathrm{n}}+\alpha\Delta
t\left(\partial_{t}\boldsymbol{B}^{\mathrm{n}}+\frac{\varsigma}{\alpha_{m}}(\partial_{t}\boldsymbol{B}^{\mathrm{n+\alpha_{m}}}-\partial_{t}\boldsymbol{B}^{\mathrm{n}})\right),$
(2b) $\displaystyle=\boldsymbol{B}^{\mathrm{n}}+\alpha\Delta
t\left(\left(1-\frac{\varsigma}{\alpha_{m}}\right)\partial_{t}\boldsymbol{B}^{\mathrm{n}}+\frac{\varsigma}{\alpha_{m}}\partial_{t}\boldsymbol{B}^{\mathrm{n+\alpha_{m}}}\right),$
(2c)

where $\alpha$, $\alpha_{m}$ and $\zeta$ are the generalized-$\alpha$
parameters which are dependent on the user-defined spectral radius
$\rho_{\infty}$:

$\displaystyle\alpha=\frac{1}{1+\rho_{\infty}},\text{
}\alpha_{m}=\frac{1}{2}\left(\frac{3-\rho_{\infty}}{1+\rho_{\infty}}\right),\text{
}\zeta=\frac{1}{2}-\alpha+\alpha_{m}.$ (3)

Using the transport equation for the left Cauchy-Green tensor, we have

$\partial_{t}\boldsymbol{B}^{\mathrm{n+\alpha_{m}}}=\nabla\boldsymbol{v}^{\mathrm{n+\alpha}}\boldsymbol{B}^{\mathrm{n+\alpha}}+\boldsymbol{B}^{\mathrm{n+\alpha}}(\nabla\boldsymbol{v}^{\mathrm{n+\alpha}})^{T}-(\boldsymbol{v}^{\mathrm{n+\alpha}}\cdot\nabla)\boldsymbol{B}^{\mathrm{n+\alpha}}.$
(4)

For simplification, we assume $\varsigma=\alpha_{m}$ i.e. $\rho_{\infty}=1$.
Substituting Equation (4) in Equation (2c) with the above simplification, we
obtain

$\boldsymbol{B}^{\mathrm{n}+\alpha}=\boldsymbol{B}^{\mathrm{n}}+\alpha\Delta
t\bigg{(}\nabla\boldsymbol{v}^{\mathrm{n+\alpha}}\boldsymbol{B}^{\mathrm{n+\alpha}}+\boldsymbol{B}^{\mathrm{n+\alpha}}(\nabla\boldsymbol{v}^{\mathrm{n+\alpha}})^{T}-(\boldsymbol{v}^{\mathrm{n+\alpha}}\cdot\nabla)\boldsymbol{B}^{\mathrm{n+\alpha}}\bigg{)}.$
(5)

The Jacobian or the tangent stiffness matrix can thus be evaluated as:

$\boldsymbol{K}=\frac{\partial\boldsymbol{\sigma}_{s}^{\prime}}{\partial\boldsymbol{v}^{\mathrm{n}+\alpha}}=\mu_{s}^{L}\left(\frac{\partial\boldsymbol{B}^{\mathrm{n}+\alpha}}{\partial\boldsymbol{v}^{\mathrm{n}+\alpha}}-\frac{1}{2}\frac{\partial(tr(\boldsymbol{B}^{\mathrm{n}+\alpha})\boldsymbol{I})}{\partial\boldsymbol{v}^{\mathrm{n}+\alpha}}\right).$
(6)

The terms for the above equation can be evaluated as shown in [54].

##  Appendix B Oscillating Disk in Taylor-Green Vortex Field

(a)

(b)

(c)

(d)

$\omega_{z}:$

Figure 1: Oscillating disk in a Taylor-Green voretx field: Snapshots of the
vorticity contours for the oscillating disk problem at $t=$ (a) $0$, (b)
$0.25$, (c) $0.5$ and (d) $1$. The black solid lines denote the $\phi=0$ iso-
contours of the two bodies. Figure 2: Comparison of temporal variation of
kinetic energy KE and strain energy SE with previous works [55, 46, 45].

We present our numerical results for the FSI case of an oscillating disk
placed in a Taylor-Green vortex field and compare it with available literature
to verify our energy computations. A unit square box over $[0,1]\times[0,1]$
with periodic boundaries is defined as the computational domain for this
problem. An incompressible hyperelastic solid is placed at the center of the
domain ($0.5,0.5$) at $t=0$. The stress-free undeformed configuration of the
solid is a circle of radius $R=0.2$. The initial condition for the phase-field
function is considered as follows:

$\phi(x,y,0)=1+\mathrm{tanh}\left(\frac{R-\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}}{\sqrt{2}\varepsilon}\right).$
(7)

The interface resolution is chosen $\frac{\varepsilon}{h}=1$. A Taylor-Green
vortex field, similar to Section 6.2, is initialized with the following
streamfunction distribution:

$\psi(x,y)=\psi_{0}\mathrm{sin}(k_{x}x)\mathrm{sin}(k_{y}y),$ (8)

where $\psi_{0}=0.05$ and $k_{x}=k_{y}=2\pi$. Other physical parameters for
the simulation are fluid viscosity $\mu_{f}=10^{-3}$, shear modulus
$\mu_{s}^{L}=1$ and density $\rho_{f}=\rho_{s}=1$. A structured grid of
$100\times 100$ and a time step size of $0.001$ is used for this problem.

Figure 1 illustrates the vorticity contours for this problem at different time
steps. The dynamics of the solid body resemble that of a damped oscillator.
The solid is deformed initially by the imposed vorticity and retracts due to
its elastic response. This generates oscillations in the system that dissipate
eventually due to the viscosity of the fluid. We track the temporal evolution
of the kinetic energy of the system and strain energy of the solid as defined
in Eqs. 37 and 38. Figure 2 shows the comparison of energy plots obtained from
our approach with previous works in the literature. This verifies our energy
calculations and highlights that our results are consistent with a variety of
other methods.

## References

  * [1] W. Lu, H.-M. Heyn, R. Lubbad, S. Løset, A large scale simulation of floe-ice fractures and validation against full-scale scenario, International Journal of Naval Architecture and Ocean Engineering 10 (3) (2018) 393–402.
  * [2] F. Li, M. Kõrgesaar, P. Kujala, F. Goerlandt, Finite element based meta-modeling of ship-ice interaction at shoulder and midship areas for ship performance simulation, Marine Structures 71 (2020) 102736.
  * [3] Y.-H. Shi, D.-Q. Yang, W.-W. Wu, Numerical analysis method of ship-ice collision-induced vibration of the polar transport vessel based on the full coupling of ship-ice-water-air, Journal of Ocean Engineering and Science.
  * [4] N. Kikuchi, J. T. Oden, Contact problems in elasticity: a study of variational inequalities and finite element methods, SIAM, 1988.
  * [5] P. Wriggers, T. A. Laursen, Computational contact mechanics, Vol. 2, Springer, 2006\.
  * [6] T. J. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai, A finite element method for a class of contact-impact problems, Computer methods in applied mechanics and engineering 8 (3) (1976) 249–276.
  * [7] P. Wriggers, T. V. Van, E. Stein, Finite element formulation of large deformation impact-contact problems with friction, Computers & Structures 37 (3) (1990) 319–331.
  * [8] R. L. Taylor, P. Papadopoulos, On a patch test for contact problems in two dimensions, Computational methods in nonlinear mechanics 690 (1991) 702.
  * [9] G. Zavarise, P. Wriggers, A segment-to-segment contact strategy, Mathematical and Computer Modelling 28 (4-8) (1998) 497–515.
  * [10] N. El-Abbasi, K.-J. Bathe, Stability and patch test performance of contact discretizations and a new solution algorithm, Computers & Structures 79 (16) (2001) 1473–1486.
  * [11] G. Zavarise, L. De Lorenzis, The contact patch test for linear contact pressure distributions in 2d frictionless contact, Trends in Computational Contact Mechanics (2011) 79–100.
  * [12] F. B. Belgacem, P. Hild, P. Laborde, The mortar finite element method for contact problems, Mathematical and Computer Modelling 28 (4-8) (1998) 263–271.
  * [13] T. McDevitt, T. Laursen, A mortar-finite element formulation for frictional contact problems, International Journal for Numerical Methods in Engineering 48 (10) (2000) 1525–1547.
  * [14] P. Farah, A. Popp, W. A. Wall, Segment-based vs. element-based integration for mortar methods in computational contact mechanics, Computational Mechanics 55 (2015) 209–228.
  * [15] T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement, Computer methods in applied mechanics and engineering 194 (39-41) (2005) 4135–4195.
  * [16] J. A. Cottrell, T. J. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [17] L. De Lorenzis, P. Wriggers, T. J. Hughes, Isogeometric contact: a review, GAMM-Mitteilungen 37 (1) (2014) 85–123.
  * [18] F. Fahrendorf, L. De Lorenzis, The isogeometric collocated contact surface approach, Computational Mechanics 70 (4) (2022) 785–802.
  * [19] I. Temizer, P. Wriggers, T. Hughes, Contact treatment in isogeometric analysis with nurbs, Computer Methods in Applied Mechanics and Engineering 200 (9-12) (2011) 1100–1112.
  * [20] P. Otto, L. De Lorenzis, J. F. Unger, Explicit dynamics in impact simulation using a nurbs contact interface, International Journal for Numerical Methods in Engineering 121 (6) (2020) 1248–1267.
  * [21] D. I. Levin, J. Litven, G. L. Jones, S. Sueda, D. K. Pai, Eulerian solid simulation with contact, ACM Transactions on Graphics (TOG) 30 (4) (2011) 1–10.
  * [22] F. Lorez, M. Pundir, D. S. Kammer, Eulerian framework for contact between solids represented as phase fields, Computer Methods in Applied Mechanics and Engineering 418 (2024) 116497.
  * [23] B. Valkov, C. H. Rycroft, K. Kamrin, Eulerian method for multiphase interactions of soft solid bodies in fluids, Journal of Applied Mechanics 82 (4) (2015) 041011.
  * [24] S. Frei, Eulerian finite element methods for interface problems and fluid-structure interactions, Ph.D. thesis (2016).
  * [25] S. Le Grand, Broad-phase collision detection with CUDA, GPU gems 3 (2007) 697–722.
  * [26] Y. Teng, D. I. Levin, T. Kim, Eulerian solid-fluid coupling, ACM Transactions on Graphics (TOG) 35 (6) (2016) 1–8.
  * [27] E. Burman, M. A. Fernández, S. Frei, A Nitsche-based formulation for fluid-structure interactions with contact, ESAIM: Mathematical Modelling and Numerical Analysis 54 (2) (2020) 531–564.
  * [28] B. Rath, X. Mao, R. K. Jaiman, An interface preserving and residual-based adaptivity for phase-field modeling of fully Eulerian fluid-structure interaction, Journal of Computational Physics 488 (2023) 112188.
  * [29] X. Mao, R. Jaiman, An interface and geometry preserving phase-field method for fully Eulerian fluid-structure interaction, Journal of Computational Physics 476 (2023) 111903.
  * [30] J. Chung, G. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method.
  * [31] K. E. Jansen, C. H. Whiting, G. M. Hulbert, A generalized-$\alpha$ method for integrating the filtered Navier–Stokes equations with a stabilized finite element method, Computer methods in applied mechanics and engineering 190 (3-4) (2000) 305–319.
  * [32] F. Shakib, T. J. Hughes, Z. Johan, A new finite element formulation for computational fluid dynamics: X. the compressible Euler and Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 89 (1-3) (1991) 141–219.
  * [33] C. Johnson, Numerical solution of partial differential equations by the finite element method, Courier Corporation, 2012.
  * [34] C. H. Rycroft, C.-H. Wu, Y. Yu, K. Kamrin, Reference map technique for incompressible fluid–structure interaction, Journal of Fluid Mechanics 898 (2020) A9.
  * [35] C. Liu, N. J. Walkington, An Eulerian description of fluids containing visco-elastic particles, Archive for rational mechanics and analysis 159 (3) (2001) 229–252.
  * [36] J. A. Trangenstein, P. Colella, A higher-order godunov method for modeling finite deformation in elastic-plastic solids, Communications on Pure and Applied mathematics 44 (1) (1991) 41–100.
  * [37] K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi, Y. Matsumoto, A full eulerian finite difference approach for solving fluid–structure coupling problems, Journal of Computational Physics 230 (3) (2011) 596–627.
  * [38] V. Joshi, R. K. Jaiman, C. Ollivier-Gooch, A variational flexible multibody formulation for partitioned fluid–structure interaction: Application to bat-inspired drones and unmanned air-vehicles, Computers & Mathematics with Applications 80 (12) (2020) 2707–2737.
  * [39] A. Džiugys, B. Peters, An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers, Granular matter 3 (2001) 231–266.
  * [40] J. A. Sethian, A fast marching level set method for monotonically advancing fronts., proceedings of the National Academy of Sciences 93 (4) (1996) 1591–1595.
  * [41] J. Parvizian, A. Düster, E. Rank, Finite cell method, Computational Mechanics 41 (1) (2007) 121–133.
  * [42] W. A. Wall, P. Gamnitzer, A. Gerstenberger, Fluid–structure interaction approaches on fixed grids based on two different domain decomposition ideas, International Journal of Computational Fluid Dynamics 22 (6) (2008) 411–427.
  * [43] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A hybrid particle level set method for improved interface capturing, Journal of Computational physics 183 (1) (2002) 83–116.
  * [44] S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of computational physics 100 (1) (1992) 25–37.
  * [45] Y. Bhosale, T. Parthasarathy, M. Gazzola, A remeshed vortex method for mixed rigid/soft body fluid–structure interaction, Journal of Computational Physics 444 (2021) 110577.
  * [46] S. S. Jain, K. Kamrin, A. Mani, A conservative and non-dissipative eulerian formulation for the simulation of soft solids in fluids, Journal of Computational Physics 399 (2019) 108922.
  * [47] A. M. Menzel, Tuned, driven, and active soft matter, Physics reports 554 (2015) 1–45.
  * [48] D. Trivedi, C. D. Rahn, W. M. Kier, I. D. Walker, Soft robotics: Biological inspiration, state of the art, and future research, Applied bionics and biomechanics 5 (3) (2008) 99–117.
  * [49] B. E. Griffith, N. A. Patankar, Immersed methods for fluid–structure interaction, Annual review of fluid mechanics 52 (2020) 421–448.
  * [50] B. Liu, D. Tan, A nitsche stabilized finite element method for embedded interfaces: Application to fluid-structure interaction and rigid-body contact, Journal of Computational Physics 413 (2020) 109461.
  * [51] R. Glowinski, T.-W. Pan, T. I. Hesla, D. D. Joseph, J. Periaux, A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow, Journal of computational physics 169 (2) (2001) 363–426.
  * [52] A. D. Rakotonirina, A. Wachs, Grains3d, a flexible dem approach for particles of arbitrary convex shape-part ii: Parallel implementation and scalable performance, Powder technology 324 (2018) 18–35.
  * [53] M. Reder, D. Schneider, F. Wang, S. Daubner, B. Nestler, Phase-field formulation of a fictitious domain method for particulate flows interacting with complex and evolving geometries, International Journal for Numerical Methods in Fluids 93 (8) (2021) 2486–2507.
  * [54] X. Mao, V. Joshi, R. Jaiman, A variational interface-preserving and conservative phase-field method for the surface tension effect in two-phase flows, Journal of Computational Physics 433 (2021) 110166.
  * [55] H. Zhao, J. B. Freund, R. D. Moser, A fixed-mesh method for incompressible flow–structure systems with finite solid deformations, Journal of Computational Physics 227 (6) (2008) 3114–3140.