3.5cm3.5cm3cm3cm

# Implicit like time discretization for the one-phase Hele-Shaw problem with
surface tension

Ido Lavi , Nicolas Meunier , Olivier Pantz Departament de Física de la Matèria
Condensada, Universitat de Barcelona, Avinguda Diagonal 647, 08028 Barcelona,
SpainLaMME, UMR CNRS 8071, Université Évry Val d’Essonne, France,
<EMAIL_ADDRESS>Jean Alexandre Dieudonné, UMR
CNRS 7351, Université de Nice<EMAIL_ADDRESS>
###### Abstract

In this work, we propose and compare three numerical methods to handle the
one-phase Hele-Shaw problem with surface tension in dimension two by using
three variational approaches in the spirit of the seminal works [17, 13].

##  1 Introduction

Consider the classical experiment: a droplet of viscous fluid is trapped
between two narrowly spaced horizontal glass plates (the Hele-Shaw cell see
Fig. 1). Because of this special geometry, the motion of the viscous fluid is
strongly overdamped; it is assumed to be governed by Darcy’s law. The effect
of surface tension at the interface is such that the droplet is at rest if and
only if its cross-section is circular.

Figure 1: Viscous droplet confined between two parallel plates.

In equations, the continuous time moving-boundary-value problem associated
with the confined droplet is the following. Let $\Omega_{0}=\Omega(t=0)$ be a
given smooth (convex and bounded) domain of $\mathbb{R}^{2}$. For all $t>0$,
the velocity $\mathbf{u}(t,\mathbf{x})$, the pressure $p(t,\mathbf{x})$ and
the domain describing the droplet $\Omega(t)$ together with its boundary
$\Gamma(t):=\partial\Omega(t)$ satisfy

$\mathbf{u}+\nabla p=0,$

with $p$ solution of

(1.1) $\begin{cases}\Delta p=0&\text{in }\Omega(t)\\\ p=\sigma\kappa&\text{on
}\Gamma(t)\end{cases}$

together with the kinematic condition, stating that the sharp interface is
transported by the fluid

(1.2)
$V_{\text{n}}=V_{\Gamma(t)}\cdot\mathbf{n}=\mathbf{u}\cdot\mathbf{n}\quad\text{on
}\Gamma(t),$

where $\kappa$ is the (mean) curvature (positive for a circle) of the evolving
free-boundary $\Gamma(t)$, $V_{n}$ and $\mathbf{n}$ are the normal velocity
field and the outer unit normal field over $\Gamma(t)$, and $\sigma>0$ the
surface tension is a given constant.

Equations (1.1) – (1.2) describe a situation where a partial differential
equation is solved for an unknown function, the fluid velocity
$\mathbf{u}(\mathbf{x},t)$ and the pressure $p(\mathbf{x},t)$, but at the same
time the exact extent of the computational domain is also unknown. This is a
free boundary problem. Consequently, one of the computational tasks will be to
treat the domain shape throughout the simulation.

It is well-known that the perimeter $\mathcal{P}(\Omega)$ is a Liapunov
functional for the evolution problem (1.1) – (1.2), see [2] e.g., that is

$\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}\mathcal{P}(\Omega(t))=-\frac{1}{\sigma}\int_{\Omega(t)}|\nabla
p|^{2}\operatorname{d}\\!x.$

Furthermore, in [17, 13], it was proved that the evolution problem (1.1) –
(1.2) can be understood as the gradient flow of the perimeter. It can be
approximated by the so-called JKO scheme, which, starting from an initial data
$\Omega_{0}$, is given by a sequence $k\in\mathbf{N}$ of variational problems
of the form

$\displaystyle\Omega^{\tau}_{(k+1)}\textrm{ minimizes }$
$\displaystyle\qquad\qquad\qquad
W_{2}^{2}(\chi_{\Omega},\chi_{\Omega^{\tau}_{(k)}})+2\tau\mathcal{P}(\Omega)$
(1.3) $\displaystyle\qquad\qquad\qquad\mbox{ among all the sets of finite
perimeter with }|\Omega|=|\Omega_{0}|,$

where $k$ is the time step, $\tau$ is its size,
$\chi_{\Omega}\in\mathrm{BV}(\Omega;\\{0,1\\})$ is the characteristic function
of the set $\Omega$, that is

$\chi_{\Omega}(x)=\begin{cases}1&\textrm{ in }\Omega\,,\\\ 0&\textrm{
else}\,,\end{cases}$

$|\Omega|$ is the area of $\Omega$ and $W_{2}$ is the usual Wasserstein
distance, defined by

(1.4) $\displaystyle
W^{2}_{2}\left(\chi_{\Omega}-\chi_{\Omega^{\prime}}\right)$ $\displaystyle=$
$\displaystyle\inf_{\pi\in\Pi(\mu,\nu)}\iint_{\Omega\times\Omega}|x-y|^{2}\operatorname{d}\\!\pi(x,y)$
$\displaystyle=$
$\displaystyle\inf\int_{\Omega}|x-T(x)|^{2}\operatorname{d}\\!x$

where $\Pi(\mu,\nu)$ denotes the set of all probability measures
$\pi\in\mathbb{P}(\Omega\times\Omega)$ with marginals $\mu$ and $\nu$ and
where in the second line the infimum runs over all transport maps, i.e.,
measure preserving $T:\mathbb{R}^{2}\to\mathbb{R}^{2}$ such that
$T_{\sharp\chi_{\Omega}}=\chi_{\Omega^{\prime}}$, that is

$\int_{\mathbb{R}^{2}}T(x)\chi_{\Omega}(x)\operatorname{d}\\!x=\int_{\mathbb{R}^{2}}\chi_{\Omega^{\prime}}(y)\operatorname{d}\\!y\,.$

Under the assumption that the perimeter $\mathcal{P}(\Omega_{k})$ converges
toward the perimeter $\mathcal{P}(\Omega)$, by using the JKO scheme defined by
(1.14) in [13], it was possible to prove the existence of weak solution to
(1.1) – (1.2), see [17, 13]. Therefore, the natural numerical method would be
to use this scheme. But it is very nonlinear. In this work, we present
numerical discretziation schemes by using sequences of variational problems in
the spirit of (1.3) but which are less nonlinear and we compare their results.
From a practical point of view, there are several difficulties to solve:

  1. 1.

Evaluation of the Wasserstein distance.

  2. 2.

Which discretization in space?

  3. 3.

How to impose the incompressibility constraint?

Throughout this work, we will consider smooth domains and we will use for the
Wasserstein distance the formulation that is based on the minimization of a
diffeomorphism between the two forms whose distance we want to evaluate.

Let us give a brief overview of the content of this article. In Section 2 we
present the different time discrete problems associated with (2.1) – (2.3). In
Section 3, we recall some basic facts that will be useful. In Sections 4 and
5, we study the explicit scheme and (2.10) – (2.11). In Section 6 we study
(2.12) – (2.13) and in Section 7 we study (2.15). Finally we end with
numerical simulations in Section 8.

##  2 Different time-discretisations of (2.1) – (2.3)

In free boundary problems, representations of interface and bulk fields are
coupled through: (i) interface kinematics: the transport of the Lagrangian or
Eulerian interface description by the Eulerian velocity field, or (ii)
interface dynamics: problems with boundary (or jump) conditions associated
with the sharp interface. The solution of the interface kinematic problem has
seen major progress in the past twenty years (see [18, 14]). However, the
situation for interface dynamics, and surface tension in particular, is more
complex and a wide range of methods or their combinations have been proposed,
see, e.g., [12].

Typically, Hele-Shaw flow problems with a sharp moving interface are solved
numerically (if not analytically) by taking advantage of elegant meshless
techniques such as conformal mapping and the vortex-sheet method (see [11] for
a review).

However, to simulate the type of problem that is involved in some recent
developments in mathematical biology, in which the fluid flow ($\mathbf{u}$)
is coupled to the bulk dynamics of some concentration field ($c$), see e.g.
[19, 1, 5, 6, 7, 10, 16], it is crucial to determine, at each time step, the
deformed geometry on which $c$ and $\mathbf{u}$ are defined. In principle, one
could avoid explicitly tracing the interface by using a phase-field (or level-
set) method, but this entails formulating a new model which would only
approximate our equations of motion at a computationally-expensive limit [9].
In our attempt to maximize precision at efficient computation costs, we chose
to build a dynamic-interface simulation based on the finite element method
(FEM), which is known as one of the most powerful numerical techniques for
solving PDEs on arbitrary domains. In a further work [15], we will apply the
three methods described here to the problems described in [19, 1, 5, 6, 7, 10,
16].

Rewriting the equation (1.1) in a slightly different form

(2.1) $\begin{cases}\mathbf{u}+\nabla p=0&\text{in }\Omega(t)\,,\\\
\nabla\cdot\mathbf{u}=0&\text{in }\Omega(t)\,,\\\ p=\sigma\kappa&\text{on
}\Gamma(t)\,,\end{cases}$

allows us to derive a mixed variational formulation of the PDE model in which
the pressure boundary condition is included as a natural boundary condition.

Let us first detail the role played by the normal component of the velocity on
the boundary. Let $\mathbf{X}(\cdot,t):I\to\Gamma(t)$ be a parameterization of
$\Gamma(t)$. The kinematic condition (1.2) writes

(2.2) $\begin{cases}\mathbf{X}(\cdot,t)&=\Gamma(t)\,,\\\
\frac{\operatorname{d}\\!\mathbf{X}(s,t)}{\operatorname{d}\\!t}&=\left(\mathbf{u}\left(\mathbf{X}(s,t)\right)\cdot\mathbf{n}\left(\mathbf{X}(s,t)\right)\right)\mathbf{n}\left(X(s,t)\right)+f(s,t)\mathbf{\boldsymbol{\tau}}\left(X(s,t)\right)\,,\end{cases}$

where $s$ is the parameterization variable, $I$ is the parameterization
interval, $\mathbf{n}$ is the outwards pointing normal vector of $\Gamma$,
$\mathbf{\boldsymbol{\tau}}$ is the tangent vector, and $f(s,t)$ is any smooth
function. According to (2.2), the time-varying set of points $\Gamma(t)$ only
depends on the normal component of the velocity ($\mathbf{u}\cdot\mathbf{n}$).
Indeed, any movement along the tangent $\mathbf{\boldsymbol{\tau}}$ simply
serves to re-parameterize $\Gamma$. In other words, the shape of the deformed
interface is determined by the normal velocity alone. Therefore, the
tangential flow component on the boundary
($\mathbf{u}\cdot\mathbf{\boldsymbol{\tau}}$) is completely irrelevant for the
continuous problem. In fact, we can modify the equation of motion of
$\Gamma(t)$ so that it is

(2.3)
$\frac{\operatorname{d}\\!\mathbf{X}(s,t)}{\operatorname{d}\\!t}=\mathbf{u}\left(X(s,t)\right).$

Much of the computational accuracy and stability depends on the discretization
of the term responsible for the surface tension force in (2.1). For most
numerical approaches the evaluation of the curvature $\kappa\mathbf{n}$ is
very difficult as it contains second derivatives.

###  2.1 Explicit time-discretisation of (2.1) – (2.3)

Knowing the domain configuration at $t=t_{n}$, and seeking the solution at
$t_{n+1}$ of (2.1), an explicit treatment of the surface tension term implies
that the curvature vector $\mathbf{H}:=\kappa\mathbf{n}$ is obtained using the
known domain configuration $X(\cdot,t_{n})$. For simplicity, we only consider
one time-step $\Delta t>0$ of the time-discrete problem. This allows us to
drop the time index notation and we omit the superscript $i$ for brevity.

More precisely, let $\Omega$ and $\Gamma$ be the domain and its interface at
time $t_{n}=n\Delta t$, $n\in\mathbf{N}$. The explicit time-discrete version
of (1.2) – (2.1) is: find $\mathbf{u}\in\left(L^{2}(\Omega)\right)^{2}$ and
$p\in H^{1}(\Omega)$ solutions of

(2.4) $\begin{cases}\mathbf{u}+\nabla p=0&\text{in }\Omega\,,\\\
\nabla\cdot\mathbf{u}=0&\text{in }\Omega\,,\\\
p\mathbf{n}=\sigma\mathbf{H}&\text{on }\Gamma\,,\end{cases}$

where $\mathbf{H}$ and $\mathbf{n}$ are the curvature vector and the unit
outwards normal vector of $\Gamma$.

According to the kinematic condition, (1.2), the new domain is defined by

(2.5) $\Omega_{\varphi}=(\textrm{Id}+\Delta
t\,\mathbf{u})(\Omega):=\varphi(\Omega)\,,$

and

(2.6) $\mathbf{X}_{\varphi}(s):=\varphi\left(\mathbf{X}(s)\right)\,,$

where $\mathbf{u}$ has been extended to $\overline{\Omega}$. Note here that
$s=\left(\mathbf{X}\right)^{-1}(\mathbf{x})$ for all $\mathbf{x}\in\Gamma$,
hence (2.6) re-writes as

(2.7)
$\mathbf{X}_{\varphi}\left(\left(\mathbf{X}\right)^{-1}(\mathbf{x})\right)=\mathbf{x}+\Delta
t\,\mathbf{u}(\mathbf{x})\,,\quad\textrm{ for all }\mathbf{x}\in\Gamma\,.$

The use of a finite time-step for the new position of the interface introduces
a time-discretization error (that is a difference between
$\mathbf{X}_{\varphi}$ given by (2.6) and the exact position where the
interface ”should” be for the continuous time case (1.2)) which is directly
related to the values of $\Delta t$ and of the velocity $\mathbf{u}$.

In (2.4), the unknowns ($\mathbf{u}$, $p$) are explicit, meaning that their
computation is done from the geometry of the previous time, that is $\Omega$.

The solution to (2.4) is the solution of the minimization problem

$\inf_{\mathbf{u}\in V(\Omega)}J_{\rm explicit}(\mathbf{u}),$

where

(2.8) $J_{\rm
explicit}(\mathbf{u})=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\operatorname{d}\\!x+\frac{\sigma}{\Delta
t}\mathcal{P}(\Omega)+\sigma\langle\mathcal{P}^{\prime}(\Omega),\mathbf{u}\rangle\,,$

and $V(\Omega)$ is the functional space

(2.9)
$V(\Omega)=\left\\{\mathbf{u}\in\left(L^{2}\left(\Omega\right)\right)^{2}\textrm{
such that }\mathrm{div}\,\,\mathbf{u}=0\right\\}\,.$

###  2.2 Other formulations

The main result of the present work is the construction of numerical solutions
in the spirit of the implicit time discretization (1.3) proposed in [17, 13].
To do so, we consider three different directions: a fully boundary
minimization problem, a curl penalization and a fully non-linear treatment of
the incompressibility constraint.

####  2.2.1 Boundary variation problem

The boundary minimization problem is:

(2.10) $\inf_{\psi\in BV(\Gamma)}J_{\rm mod}(\psi),$

with

(2.11) $J_{\rm
mod}(\psi):=\|\psi\cdot\mathbf{n}\|_{H^{-1/2}(\Gamma)}+V\left((\textrm{Id}+\Delta
t\,\psi)(x)\right),$

where $V(\cdot)$ is the total variation and $BV$ is the space of functions
with bounded variations.

####  2.2.2 A problem with curl penalization

The minimization problem with a curl penalization is:

(2.12) $\inf_{\mathbf{u}\in V(\Omega)}J_{\alpha}(\Omega,\mathbf{u})\,,$

with $\alpha>0$ and

(2.13)
$J_{\alpha}(\Omega,\mathbf{u})=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\,\operatorname{d}\\!x+\alpha\int_{\Omega}|\nabla\wedge\mathbf{u}|^{2}\,\operatorname{d}\\!x+\frac{\sigma}{\Delta
t}\mathcal{P}\left((\textrm{Id}+\Delta t\,\mathbf{u})(\Omega)\right)\,.$

######  Remark 2.1.

As we will see below, one can not take $\alpha=0$. Indeed the minimization
problem

(2.14) $\inf_{\mathbf{u}\in V(\Omega)}J_{0}(\Omega,\mathbf{u})\,,$

is ill-posed.

####  2.2.3 A fully nonlinear problem

The fully nonlinear minimization problem is

(2.15) $\inf_{\mathbf{u}\in W(\Omega)}J(\Omega,\mathbf{u})\,,$

where $W(\Omega)$ is the functional space

(2.16)
$W(\Omega)=\left\\{\mathbf{u}\in\left(H^{1}\left(\Omega\right)\right)^{2}\textrm{
such that }\det\left(\textrm{Id}+\Delta t\,\nabla\mathbf{u}\right)=1\text{ in
}\Omega\right\\}\,.$

##  3 Some material

In this part, we recall some very classical facts on differential geometry, on
eulerian derivatives and on functions with bounded variations that will be
used later on.

###  3.1 Computation of the first variation of the deformed permiter

In order to define the eulerian derivative of the perimeter we embed problem
(2.1) into a family of perturbed problems which are defined on perturbations
of a $\mathcal{C}^{2,1}$ reference domain $\Omega$ constructed by perturbing
the identity.

Let $U$ be a convex bounded domain of class $\mathcal{C}^{2,1}$ such that
$\bar{\Omega}\subset U$ and let

(3.1)
$\mathcal{S}=\\{\mathbf{u}\in\mathcal{C}^{2,1}(\bar{U},\mathbb{R}^{2})\,:\,\mathbf{u}=0\textrm{
on }\partial U\textrm{ and }\mathbf{u}\cdot\mathbf{n}=0\textrm{ on }\Gamma\\}$

be the space of feasible deformation fields endowed with the natural norm in
$C^{2}(\bar{U},\mathbb{R}^{2})$. For a fixed field $\mathbf{u}\in\mathcal{S}$
and for all $\Delta t>0$ define the mapping from $\bar{U}$ to $\mathbb{R}^{2}$
by $\varphi=\textrm{Id}+\Delta t\mathbf{u}$. For $\Delta t$ sufficiently small
$\varphi$ defines a family of $\mathcal{C}^{2}$-diffeomorphisms of $U$ onto
itself. For such $\Delta t$ one sets

(3.2) $\Omega_{\varphi}=\varphi(\Omega)=(\textrm{Id}+\Delta
t\,\mathbf{u})(\Omega)\subset\mathbb{R}^{2}\,,$

and

$\Gamma_{\varphi}=\varphi(\Gamma),$

hence $\Omega_{0}=\Omega$, $\Gamma_{0}=\Gamma$ and
$\partial\varphi(\Omega):=\varphi(\partial\Omega)$.

We define the matrix $\nabla\varphi$ by
$\left(\nabla\varphi\right)_{ij}=\partial_{j}\varphi_{i}$ for $1\leq i,j\leq
2$. This matrix represents the differential of $\varphi$ in the sense that
$\varphi(\mathbf{x}+h)-\varphi(\mathbf{x})=\nabla\varphi(\mathbf{x})h+o(\|h\|)$
for $\mathbf{x}$ and $\mathbf{x}+h$ in $\Omega$.

####  3.1.1 On the curvature vectors of $\partial\Omega$ and of
$\varphi(\partial\Omega)$

Let $\mathbf{\boldsymbol{\tau}}$ and $\mathbf{\boldsymbol{\tau}}_{\varphi}$
denote the tangent vectors along $\partial\Omega$ and
$\partial\varphi(\Omega)$ respectively, and let $\mathbf{T}$ represent the
tangent vector on $\partial\varphi(\Omega)$ ”pulled back” to the reference
domain $\partial\Omega$. In other words, the input of $\mathbf{T}$ is any
$\mathbf{x}$ on $\partial\Omega$ and the output is
$\mathbf{\boldsymbol{\tau}}_{\varphi}(\mathbf{x}_{\varphi})$ where
$\mathbf{x}_{\varphi}:=\mathbf{x}+\Delta t\,\mathbf{u}(\mathbf{x})$ and

(3.3)
$\mathbf{T}:=\mathbf{\boldsymbol{\tau}}_{\varphi}\circ\varphi=\frac{\nabla\varphi\,\mathbf{\boldsymbol{\tau}}}{|\nabla\varphi\,\mathbf{\boldsymbol{\tau}}|}=\frac{(\textrm{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}}{|(\textrm{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|}\,,$

where I represents the $2\times 2$ identity matrix and we use the notation
$\nabla\varphi\,\mathbf{\boldsymbol{\tau}}=(\mathbf{\boldsymbol{\tau}}\cdot\nabla)\varphi$.

Let $\mathbf{X}(s)$ for $s\in I$ be a parameterization of $\partial\Omega$ and
$\mathbf{X}_{\varphi}(s)$ be a parametrization of the interface
$\partial\varphi(\Omega)$:

(3.4) $\mathbf{X}_{\varphi}(s)=\varphi(\mathbf{X}(s))=\mathbf{X}(s)+\Delta
t\,\mathbf{u}\left(\mathbf{X}(s)\right)\,.$

Using differential geometry tools [20], the curvature vectors
$\mathbf{H}:=\kappa\mathbf{n}$ on $\partial\Omega$ and
$\mathbf{H}_{\varphi}:=\kappa_{\varphi}\mathbf{n}_{\varphi}$ on
$\partial\varphi(\Omega)$ are defined by

$\mathbf{H}=-|\mathbf{X}^{{}^{\prime}}(s)|^{-1}\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\mathbf{\boldsymbol{\tau}}(\mathbf{X}(s))=-\nabla\mathbf{\boldsymbol{\tau}}\,\mathbf{\boldsymbol{\tau}}\,,$
$\mathbf{H}_{\varphi}=-|\mathbf{X}_{\varphi}^{{}^{\prime}}(s)|^{-1}\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\mathbf{\boldsymbol{\tau}}_{\varphi}(\mathbf{X}_{\varphi}(s))=-\nabla\mathbf{\boldsymbol{\tau}}_{\varphi}\,\mathbf{\boldsymbol{\tau}}_{\varphi}\,.$

Define the vector $\widetilde{\mathbf{H}_{\varphi}}$ by

(3.5) $\displaystyle\widetilde{\mathbf{H}_{\varphi}}$
$\displaystyle=-|\mathbf{X}^{{}^{\prime}}(s)|^{-1}\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\left(\mathbf{T}(\mathbf{X}(s))\right)=-\nabla\mathbf{T}\,\mathbf{\boldsymbol{\tau}}\,.$

Let $\operatorname{Cof}\nabla\varphi$ be the cofactor matrix associated with
the gradient matrix $\nabla\varphi$. Since

$\operatorname{Cof}\nabla\varphi^{T}\nabla\varphi=\det(\nabla\varphi)\textrm{Id}\,,$

we see that

$\widetilde{\mathbf{H}_{\varphi}}=|\operatorname{Cof}\nabla\varphi|\,\mathbf{H}_{\varphi}\circ\varphi\,.$

####  3.1.2 First variation of the perimeter
$\mathcal{P}\left(\varphi(\Omega)\right)$

The Eulerian derivative of the permiter $\mathcal{P}$ at $\Omega$ in the
direction $\mathbf{u}$ is defined as

(3.6)
$\langle\mathcal{P}^{\prime}\left(\Omega\right);\mathbf{u}\rangle:=\lim_{\Delta
t\to
0}\frac{\mathcal{P}\left(\varphi(\Omega)\right)-P\left(\Omega\right)}{\Delta
t}\,.$

In order to evaluate the first variation of $P(\Omega_{\varphi})$ we transform
$P(\Omega_{\varphi})$ to an integral over the reference domain $\Omega$. The
general strategy we use in this work is adapted from the one used for shape
derivative calculations. It consists in transferring the problem on the
original boundary before writing the variational problem.

######  Lemma 3.1.

The first variation of the perimeter $\mathcal{P}\left(\varphi(\Omega)\right)$
is given by

(3.7)
$\langle\mathcal{P}^{\prime}\left(\varphi(\Omega)\right);\mathbf{v}\rangle=\Delta
t\int_{\partial\varphi(\Omega)}\mathbf{H}_{\varphi}\cdot(\mathbf{v}\circ\varphi^{-1})\operatorname{d}\\!s_{\varphi}=\Delta
t\int_{\partial\Omega}\widetilde{\mathbf{H}_{\varphi}}\cdot\mathbf{v}\operatorname{d}\\!s\,.$

###### Proof.

Assume that the boundary $\partial\Omega$ is parametrized by the arc length
$s$, we have

(3.8)
$\mathcal{P}\left(\varphi(\Omega)\right)=\int_{\partial\Omega}|\nabla\varphi\,\mathbf{\boldsymbol{\tau}}|\operatorname{d}\\!s=\int_{\partial\Omega}|(\textrm{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|\operatorname{d}\\!s\,.$

Let us perform an asymptotic expansion of order one of
$\mathcal{P}\left(\varphi(\Omega)\right)$. We first find

$\displaystyle|(\textrm{I}+\Delta
t\,\nabla(\mathbf{u}+\delta\mathbf{u}))\mathbf{\boldsymbol{\tau}}|^{2}$
$\displaystyle=$ $\displaystyle|(\textrm{I}+\Delta
t\,\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|^{2}+2\Delta
t\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot(\mathbf{\boldsymbol{\tau}}+\Delta
t\,\nabla\mathbf{u}\,\mathbf{\boldsymbol{\tau}})+O(||\nabla\delta\mathbf{u}||^{2})\,.$

Let $\operatorname{d}\\!S=|(\textrm{I}+\Delta
t\,\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|$, using the Taylor expansion
of $\sqrt{\operatorname{d}\\!S^{2}+x}$, we obtain

$|(\textrm{I}+\Delta
t\,\nabla(\mathbf{u}+\delta\mathbf{u}))\mathbf{\boldsymbol{\tau}}|=|(\textrm{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|+\Delta
t\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}+O(||\nabla\delta\mathbf{u}||^{2})\,,$

where we recall the definition of $\mathbf{T}$:

$\mathbf{T}=\frac{\nabla\varphi\,\mathbf{\boldsymbol{\tau}}}{|\nabla\varphi\,\mathbf{\boldsymbol{\tau}}|}=\frac{(\textrm{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}}{|(\textrm{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|}\,,$

and thus

$\mathcal{P}((\textrm{Id}+\Delta t\,\mathbf{u}+\Delta
t\delta\mathbf{u})(\Omega))=\mathcal{P}(\varphi(\Omega))+\Delta
t\int_{\partial\Omega}\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}\operatorname{d}\\!s+O(||\nabla\delta\mathbf{u}||^{2})\,.$

Consequently, for $\mathbf{v}:\partial\Omega\to\mathbb{R}^{2}$ we have

(3.9)
$\langle\mathcal{P}^{\prime}\left(\varphi(\Omega)\right);\mathbf{v}\rangle=\Delta
t\int_{\partial\Omega}\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}\operatorname{d}\\!s\,.$

Recalling that $s$ is the arc length coordinate of $\partial\Omega$, it
follows that
$\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}=\frac{\operatorname{d}\\!\mathbf{v}}{\operatorname{d}\\!s}$.
Then, integrating (3.1.2) by parts we obtain

$\langle\mathcal{P}^{\prime}\left(\varphi(\Omega)\right);\mathbf{v}\rangle=-\Delta
t\int_{\partial\Omega}\mathbf{v}\cdot\frac{\operatorname{d}\\!\mathbf{T}}{\operatorname{d}\\!s}\operatorname{d}\\!s\,.$

Finally, since
$\widetilde{\mathbf{H}_{\varphi}}=-\frac{\operatorname{d}\\!\mathbf{T}}{\operatorname{d}\\!s}$
we deduce (3.7). ∎

###  3.2 On functions with bounded variation

We recall the definition of the total variation for a function $f\in
L^{1}(\Gamma)$:

$V(f):=\int_{\Gamma}|\textrm{D}f|\operatorname{d}\\!s=\sup\left\\{\int_{\Gamma}f\,g^{\prime}\operatorname{d}\\!s\,;\,g\in\mathcal{C}^{1}(\Gamma,\mathbb{R}^{2})\,,\,|g(s)|\leq
1\textrm{ for all }s\in\Gamma\right\\},$

where $f^{\prime}$ is the derivative of $f$. The space of functions with
bounded variation is then defined as

$BV(\Gamma)=\left\\{f\in
L^{1}(\Gamma)\,;\,\int_{\Gamma}|\textrm{D}f|\operatorname{d}\\!s<\infty\right\\},$

and is equipped with the norm

$\|f\|_{BV(\Gamma)}=\|f\|_{L^{1}(\Gamma)}+\int_{\Gamma}|\textrm{D}f|\operatorname{d}\\!s.$

######  Remark 3.2.

If $f\in W^{1,1}(\Gamma)$, this coincide with the usual $W^{1,1}(\Gamma)$
norm. However $W^{1,1}(\Gamma)$ is a proper subset of $BV(\Gamma)$, since the
derivative of $BV$ functions are in general measures not functions in
$L^{1}(\Gamma)$.

We recall the following classical results for BV functions.

######  Theorem 3.3.

Let $(\mathbf{u}_{k})_{k\in\mathbf{N}}$ be a bounded sequence in $BV(\Gamma)$.
Then there exists a subsequence which converges strongly in $L^{1}(\Gamma)$.

A crucial lemma in the study of BV functions is the following approximation
result:

######  Lemma 3.4 ((Approximation by $\mathcal{C}^{\infty}$ functions)).

Let $\mathbf{u}\in BV(\Gamma)$. There exists a sequence
$\mathbf{u}_{k}\in\mathcal{C}^{\infty}(\Gamma)$ such that
$\mathbf{u}_{k}\to\mathbf{u}$ in $L^{1}(\Gamma)$ and

$\int_{\Gamma}|\mathbf{u}_{k}^{\prime}|\operatorname{d}\\!s\to\int_{\Gamma}|\textrm{D}\mathbf{u}_{k}|\operatorname{d}\\!s.$

######  Remark 3.5.

Note that we cannot expect to have

$\int_{\Gamma}|\mathbf{u}_{k}^{\prime}-\mathbf{u}^{\prime}|\operatorname{d}\\!s\to
0,$

in general since that would imply that $\mathbf{u}\in W^{1,1}(\Gamma)$.

##  4 Explicit time discretization of (1.2) – (2.1)

From the time-discrete equations (2.4) – (2.7), we derive the variational
formulation by standard techniques.

Assuming that $\Omega$ is smooth enough, say $C^{2}$ e.g., the varitional
formulation associated with (2.4) is:

(4.1) $\begin{cases}\textrm{find
}\mathbf{u}\in\left(L^{2}(\Omega)\right)^{2}\textrm{ and }p\in
H^{1}(\Omega)\textrm{ that solves}\\\
\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\int_{\Omega}\nabla
p\cdot\mathbf{v}\,\operatorname{d}\\!x+\sigma\int_{\Gamma}\mathbf{H}^{i}\cdot\mathbf{v}\,\operatorname{d}\\!s=0\,,\\\
\int_{\Omega}(\nabla\cdot\mathbf{u})q\,\operatorname{d}\\!x=0\,,\\\ \textrm{
for any smooth test functions }\mathbf{v}:\Omega\to\mathbb{R}^{2}\textrm{ s.t.
}\mathrm{div}\,\mathbf{v}=0\\\ \textrm{ and }q:\Omega\to\mathbb{R}\textrm{
s.t. }\mathbf{v}_{\Gamma}=0\,.\end{cases}$

Since

$\frac{1}{\Delta t}\mathcal{P}\left((\textrm{Id}+\Delta
t\,\mathbf{u})(\Omega)\right)=\frac{1}{\Delta
t}\left(\mathcal{P}\left(\Omega\right)+\Delta
t\langle\mathcal{P}^{\prime}\left(\Omega\right);\mathbf{u}\rangle\right),$

We start with an existence result for the minimization of (2.8). Recall that
the space $V(\Omega)$ is defined by (2.9).

######  Lemma 4.1.

There exists a unique solution of the minimization problem

$\inf_{\mathbf{u}\in V(\Omega)}J_{\rm explicit}(\mathbf{u}).$

###### Proof.

The proof of this result is classical. The term
$\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\operatorname{d}\\!x$ is coercive on
$\left(L^{2}(\Omega)\right)^{2}$. The set $\\{\mathrm{div}\,\,\mathbf{u}=0\\}$
is closed in $\left(L^{2}(\Omega)\right)^{2}$ weak. Finally, since
$\nabla\kappa$ belongs to $\left(L^{2}(\Omega)\right)^{2}$, the linear part of
$J_{\rm explicit}$ is continuous. ∎

######  Proposition 4.2.

$(\mathbf{u},p)\in\left(L^{2}(\Omega)\right)^{2}\times H^{1}(\Omega)$ is
solution of (2.4) iff

$\mathbf{u}=\operatorname*{arg\,min}_{\mathbf{v}\in
V(\Omega)}J(\mathbf{v})\,.$

###### Proof.

Let $f\in\left(L^{2}(\Omega)\right)^{2}$. Consider first the problem with
homogeneous Dirichlet condition

(4.2) $\begin{cases}\mathbf{u}+\nabla p=f&\text{in }\Omega\,,\\\
\nabla\cdot\mathbf{u}=0&\text{in }\Omega\,,\\\ p=0&\text{on
}\Gamma\,.\end{cases}$

Let

(4.3)
$I(\mathbf{u})=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\operatorname{d}\\!x-\int_{\Omega}f\mathbf{u}\operatorname{d}\\!x\,.$

######  Lemma 4.3.

$(\mathbf{u},p)\in\left(L^{2}(\Omega)\right)^{2}\times H^{1}_{0}(\Omega)$ is
solution of (4.2) iff

$\mathbf{u}=\operatorname*{arg\,min}_{\mathbf{v}\in
V(\Omega)}I(\mathbf{v})\,.$

###### Proof.

Let $(\mathbf{u},p)\in\left(L^{2}(\Omega)\right)^{2}\times H^{1}_{0}(\Omega)$
a solution of (4.2). For all
$\delta\mathbf{u}\in\left(L^{2}(\Omega)\right)^{2}$ such that
$\mathrm{div}\,\,\delta\mathbf{u}=0$, we have

$\displaystyle\langle I^{\prime}(\mathbf{u}),\delta\mathbf{u}\rangle$
$\displaystyle=$
$\displaystyle\int_{\Omega}\mathbf{u}\cdot\delta\mathbf{u}\operatorname{d}\\!x-\int_{\Omega}f\delta\mathbf{u}\operatorname{d}\\!x$
$\displaystyle=$ $\displaystyle-\int_{\Omega}\nabla
p\cdot\delta\mathbf{u}\operatorname{d}\\!x$ $\displaystyle=$
$\displaystyle\int_{\Omega}p\,\mathrm{div}\,\left(\delta\mathbf{u}\right)\operatorname{d}\\!x$
$\displaystyle=$ $\displaystyle 0\,,$

hence $\mathbf{u}=\operatorname*{arg\,min}_{\mathbf{v}\in
V(\Omega)}I(\mathbf{v})$.

Conversely, assume that $\mathbf{u}=\operatorname*{arg\,min}_{\mathbf{v}\in
V(\Omega)}I(\mathbf{v})$. Let

$L(\mathbf{v}):=\langle I^{\prime}(\mathbf{u}),\mathbf{v}\rangle\,,$

then $L(\mathbf{v})=0$ for all $\mathbf{v}\in\left(L^{2}(\Omega)\right)^{2}$
such that $\mathrm{div}\,\,\mathbf{v}=0$. We can apply De Rham’s theorem (see
theorem A.1) to
$L:\left(L^{2}(\Omega)\right)^{2}\to\left(L^{2}(\Omega)\right)^{2}$ and we get
that there exists $\tilde{p}\in H_{0}^{1}(\Omega)$ such that

$L(\mathbf{v})=\int_{\Omega}\nabla\tilde{p}\cdot\mathbf{v}\,\operatorname{d}\\!x\,.$

Integrating by parts it comes

$L(\mathbf{v})=-\int_{\Omega}\tilde{p}\,\mathrm{div}\,\mathbf{v}\,\operatorname{d}\\!x\,.$

Setting $p=-\tilde{p}$, we deduce that
$(\mathbf{u},p)\in\left(L^{2}(\Omega)\right)^{2}\times H^{1}_{0}(\Omega)$ is a
solution of (4.2). ∎

We now apply Lemma 4.3 to the particular case where $f=-\sigma\nabla\kappa$.
Indeed $(\mathbf{u},p)$ is a solution of (4.2) iff
$(\mathbf{u},p+\sigma\kappa)$ is a solution of (2.4).

∎

######  Lemma 4.4.

The variational formulation (4.1) is the Euler equation for the minimization
problem:

$\inf_{\mathbf{u}\in V(\Omega)}J_{\rm explicit}(\mathbf{u})\,,$

where $\Omega$ is a given smooth subset of $\mathbb{R}^{2}$ ($C^{2}$ e.g.) and
the functional $J_{\rm explicit}$ is defined by (2.8) over the functional
space $V(\Omega)$ defined by (2.9).

###### Proof.

Define the lagrangian $\mathcal{L}$ on
$\left(L^{2}\left(\Omega\right)\right)^{2}\times H^{1}\left(\Omega\right)$ by

$\mathcal{L}(\mathbf{u},p)=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\,\operatorname{d}\\!x+\int_{\Omega}\mathbf{u}\cdot\nabla
p\,\operatorname{d}\\!x+\sigma\int_{\Omega}\nabla\kappa\cdot\mathbf{u}\,\operatorname{d}\\!x\,.$

The functional $\mathcal{L}$ is Gâteaux differentiable with respect to
$(\mathbf{u},p)$. Let

$\mathbf{u}=\operatorname*{arg\,min}_{\mathbf{v}\in V(\Omega)}J_{\rm
explicit}(\mathbf{v})\,.$

Using the arguments given in the proof of Proposition 4.2, we deduce that for
all $\mathbf{v}\in V(\Omega)$ and $q\in H^{1}_{0}(\Omega)$, we have

$\langle\mathcal{L}^{\prime}(\mathbf{u},p);(\mathbf{v},q)\rangle=\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\int_{\Omega}\mathbf{v}\cdot\nabla
p\operatorname{d}\\!x+\sigma\int_{\Omega}\nabla\kappa\cdot\mathbf{v}\,\operatorname{d}\\!s+\int_{\Omega}\mathbf{u}\cdot\nabla
q\operatorname{d}\\!x=0\,.$

Integrating by parts we obtain

$\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\int_{\Omega}\mathbf{v}\cdot\nabla
p\,\operatorname{d}\\!x+\sigma\int_{\partial\Omega}\mathbf{H}\cdot\mathbf{v}\,\operatorname{d}\\!s-\int_{\Omega}\left(\mathrm{div}\,\mathbf{u}\right)q\,\operatorname{d}\\!x=0\,,$

that is (4.1). ∎

######  Remark 4.5.

Recalling that the vector curvature $\mathbf{H}:=\kappa\mathbf{n}$ of $\Gamma$
is defined by

$\mathbf{H}=-|\mathbf{X}^{{}^{\prime}}(s)|^{-1}\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\mathbf{\boldsymbol{\tau}}(\mathbf{X}(s))=-\nabla\mathbf{\boldsymbol{\tau}}\,\mathbf{\boldsymbol{\tau}}\,,$

and integrating by parts on the boundary term in (4.1), we get

(4.4)
$\int_{\Gamma^{i}}\mathbf{H}^{i}\cdot\mathbf{v}=\int_{\Gamma^{i}}\mathbf{\boldsymbol{\tau}}\cdot\nabla\left(\mathbf{\boldsymbol{\tau}}\mathbf{v}\right)\,.$

When $s$ is the arc length coordinate of $\Gamma$, it comes
$\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}=\frac{\operatorname{d}\\!\mathbf{v}}{\operatorname{d}\\!s}$.
Hence, the explicit scheme is very simple to implement but like all explicit
schemes it requires a very small time step, see Section 8.

##  5 Study of (2.10) – (2.11)

In this part we study the minimization problem (2.10) – (2.11). We start with
its origin and the existence of a solution and then we study the variational
formulation and its implementation.

###  5.1 Origin of (2.10) – (2.11)

The incompressibility constraint

$\nabla\cdot\mathbf{u}=0\textrm{ in }\Omega,$

can be rewritten as

$\int_{\Gamma}\mathbf{u}\cdot\mathbf{n}\operatorname{d}\\!s=0.$

Moreover, consider functions such that

$\int_{\Gamma}\psi\cdot\mathbf{n}\operatorname{d}\\!s=0,$

we can define a norm $\|\cdot\|_{H^{-1/2}(\Gamma)}$ on the space
$H^{-1/2}(\Gamma)$ by

$\|\mathbf{\psi}\cdot\mathbf{n}\|_{H^{-1/2}(\Gamma)}:=\inf_{\tiny\begin{array}[]{l}\nabla\cdot\mathbf{u}=0\\\
\mathbf{u}\cdot\mathbf{n}=\psi\cdot\mathbf{n}\end{array}}\|\mathbf{u}\|_{L^{2}(\Omega)}.$

Consider the minimization problem (2.10) – (2.11) that we recall

(5.1) $\inf_{\psi\in BV(\Gamma)}J_{\rm mod}(\psi),$

with

$J_{\rm
mod}(\psi):=\|\psi\cdot\mathbf{n}\|_{H^{-1/2}(\Gamma)}+V\left((\textrm{Id}+\Delta
t\,\psi)(x)\right),$

where $V(\cdot)$ is the total variation and $BV(\Gamma)$ is the space of
functions with bounded variations.

######  Lemma 5.1.

There exists a solution to (2.10).

###### Proof.

Consider a minimizing sequence $(\psi_{k})_{k\in\mathbf{N}}$ of the
minimization problem (2.10). We first note that $(\psi_{k})_{k\in\mathbf{N}}$
is bounded in $BV$, hence, up to a subsequence, $(\psi_{k})_{k}$ converges to
a limit $\psi\in BV$, and we have for almsot every point of the boundary

$\lim_{k\to\infty}(\psi_{k}\cdot\mathbf{n}_{|\Gamma})=\psi\cdot\mathbf{n}_{|\Gamma},$

hence the result. ∎

Furthermore, we can perform an analysis on the tangential component of the
solution of (5.1). In this case, it is obvious that the normal component of
$\mathbf{u}$ will be small and assuming that the solution behaves correctly,
the tangential component of $\mathbf{u}$ should also be small.

For simplicity let us write
$\mathbf{u}_{\mathbf{\boldsymbol{\tau}}}:=\mathbf{u}\cdot\mathbf{\boldsymbol{\tau}}$,
$\mathbf{u}_{\mathbf{n}}:=\mathbf{u}\cdot\mathbf{n}$ and
$\dot{f}:=\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}f$.

######  Lemma 5.2.

Assume that $\Omega$ is smooth enough, $\mathcal{C}^{2}$ for example. Let
$\mathbf{u}\in BV(\Gamma)$ be such that

$J_{\rm mod}(\mathbf{u})=\inf_{\psi\in BV(\Gamma)}J_{\rm mod}(\psi).$

Assume that $\mathbf{u}_{\mathbf{\boldsymbol{\tau}}}$ and $\mathbf{u}_{n}$ are
small. Then, almost everywhere on $\Gamma$, we have

$\kappa=0\qquad\textrm{ or
}\qquad\mathbf{u}_{\mathbf{\boldsymbol{\tau}}}=\frac{1}{\sigma\kappa}\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\mathbf{u}_{\mathbf{n}}.$

###### Proof.

Define $P$ by

$P(\psi)=\int_{\Gamma}|\mathbf{\boldsymbol{\tau}}+\dot{\psi}|\operatorname{d}\\!s,$

and consider the minimization problem

$\inf_{\psi\in BV(\Gamma)}P(\psi).$

For all $\phi=\phi_{\tau}\tau$, we have

$<P^{\prime}(\psi),\phi>=\int_{\Gamma}\frac{\mathbf{\boldsymbol{\tau}}+\dot{\psi}}{|\mathbf{\boldsymbol{\tau}}+\dot{\psi}|}\cdot\dot{\phi}\operatorname{d}\\!s.$

Integating by parts, we get

$<P^{\prime}(\psi),\phi>=-\int_{\Gamma}\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\left(\frac{\mathbf{\boldsymbol{\tau}}+\dot{\psi}}{|\mathbf{\boldsymbol{\tau}}+\dot{\psi}|}\right)\cdot\phi\operatorname{d}\\!s=0$

Hence, almost evrywhere on $\Gamma$, we deduce that

$\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\left(\frac{\mathbf{\boldsymbol{\tau}}+\dot{\psi}}{|\mathbf{\boldsymbol{\tau}}+\dot{\psi}|}\right)\cdot\mathbf{\boldsymbol{\tau}}=0.$

Therefore, there exists a funtion $\beta:\Gamma\rightarrow\mathbb{R}$ such
that

$\left(\frac{\mathbf{\boldsymbol{\tau}}+\dot{\psi}}{|\mathbf{\boldsymbol{\tau}}+\dot{\psi}|}\right)=\beta\mathbf{\boldsymbol{\tau}},$

thus,

$\dot{\psi}=(\beta|\mathbf{\boldsymbol{\tau}}+\dot{\psi}|-1)\mathbf{\boldsymbol{\tau}}.$

Consequently $\dot{\psi}$ is colinear to $\mathbf{\boldsymbol{\tau}}$, that is
$\dot{\psi}\cdot\mathbf{n}=0$.

Since

$\psi=\mathbf{u}_{\mathbf{\boldsymbol{\tau}}}\mathbf{\boldsymbol{\tau}}+\mathbf{u}_{\mathbf{n}}\mathbf{n},$

we see that

$\dot{\psi}=\dot{\mathbf{u}}_{\mathbf{\boldsymbol{\tau}}}\mathbf{\boldsymbol{\tau}}-\mathbf{u}_{\mathbf{\boldsymbol{\tau}}}\kappa\mathbf{n}+\dot{\mathbf{u}}_{\mathbf{n}}\mathbf{n}+\kappa\mathbf{u}_{\mathbf{n}}\mathbf{\boldsymbol{\tau}}.$

Recalling that the normal composant of $\dot{\psi}$ is zero, we obtain

$-\kappa\mathbf{u}_{\mathbf{\boldsymbol{\tau}}}+\dot{\mathbf{u}}_{\mathbf{n}}=0.$

∎

######  Remark 5.3.

It is not surprising to have no information on the tangential component of the
deformation on the regions where the curvature of $\Gamma$ is zero. Indeed, in
these regions, displacement along the tangent simply corresponds to a
reparameterization whose cost function is independent.

######  Remark 5.4.

Using (3.8), considering that $s$ is the arc-length and denoting
$\dot{\mathbf{u}}:=\frac{\operatorname{d}\\!}{\operatorname{d}\\!s}\mathbf{u}$,
we see that $J_{0}(\Omega,\mathbf{u})$ defined by (2.13), with $\alpha=0$,
rewrites as

$J(\Omega,\mathbf{u})=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\,\operatorname{d}\\!x+\frac{\sigma}{\Delta
t}\int_{\Gamma}|\mathbf{\boldsymbol{\tau}}+\dot{\mathbf{u}}|\,\operatorname{d}\\!s.$

Consider a minimizing sequence $(\mathbf{u}_{k})_{k\in\mathbf{N}}$ of the
minimization problem (2.14), that we recall

(5.2) $\min_{\tiny\begin{array}[]{l}\mathbf{u}\in L^{2}(\Omega)\\\
\nabla\cdot\mathbf{u}=0\end{array}}J(\Omega,\mathbf{u}).$

We first note that $(\mathbf{u}_{k})_{k\in\mathbf{N}}$ is bounded in
$H(\operatorname{div})$. Hence, up to a subsequence, the sequence
$(\mathbf{u}_{k})_{k}$ converges to a limit $\mathbf{u}\in
H(\operatorname{div})$, and we have for almsot every point of the boundary

$\lim_{k\to\infty}(\mathbf{u}_{k}\cdot\mathbf{n}_{|\Gamma})=\mathbf{u}\cdot\mathbf{n}_{|\Gamma},$

where $\mathbf{n}$ is the normal to the boundary.

However, without any control on the curl of $\mathbf{u}$, there is no reason
that the sequence $(\mathbf{u}_{k}\cdot\mathbf{\boldsymbol{\tau}}_{|\Gamma})$
converges towards $\mathbf{u}\cdot\mathbf{\boldsymbol{\tau}}_{|\Gamma}$.
Indeed as we will in Section 8, even with a frequent remeshing, we obtain a
shift between the tangential displacement of the internal nodes and the
boundary nodes. For this reasons, it is difficult to use the original
formulation numerically and we propose a second method.

As a conclusion, the perimeter minimization problem on zero divergence fields
(5.2) is ill-posed, but correctly reformulated by (5.1) on the domain
boundary. For (5.1) and for small displacements, the tangential component
depends on the derivative of the normal component. We therefore expect a loss
of regularity of the tangential component with respect to the normal component
outside the parts of the boundary where $\kappa$ is zero. We cannot expect any
regularity (other than bounded variation) of the tangential component in the
areas where $\kappa=0$. Finally, in case of discretization of the initial
problem (5.2), one should see boundary layers appearing on the tangential
component of $\mathbf{u}$ whose thickness depends on the mesh used. We will
see some situations where it is not the case in Section 8.

###  5.2 Variational formulation associated with (2.10) – (2.11)

Define the functional space $W(\Omega)$ by

$\displaystyle W(\Omega)$ $\displaystyle=$
$\displaystyle\Big{\\{}\left(\mathbf{U},\mathbf{u}\right)\in\left(L^{2}(\Omega)\right)^{2}\times\left(H^{-\frac{1}{2}}\left(\Gamma\right)\cap
BV\left(\Gamma\right)\right)$ $\displaystyle\qquad\textrm{ s.t. for all }p\in
H^{1}(\Omega)\,\int_{\Gamma}p\mathbf{u}\cdot\mathbf{n}\,\operatorname{d}\\!s=\int_{\Omega}\mathbf{U}\cdot\nabla
p\,\operatorname{d}\\!x\Big{\\}},$

and the functional $J^{*}_{\rm mod}$ by

(5.3) $J^{*}_{\rm
mod}\left(\mathbf{U},\mathbf{u}\right):=\frac{1}{2}\|\mathbf{U}\|^{2}_{L^{2}(\Omega)}+\int_{\Gamma}\left|\mathbf{\boldsymbol{\tau}}+\Delta
t\,\dot{\mathbf{u}}\right|\,\operatorname{d}\\!s.$

######  Proposition 5.5.

The minimization problem (2.10) – (2.11) is equivalent to the following one:

(5.4) $\inf_{\left(\mathbf{U},\mathbf{u}\right)\in W(\Omega)}J^{*}_{\rm
mod}\left(\mathbf{U},\mathbf{u}\right).$

###### Proof.

The proof relies on the following equivalence:

$\int_{\Gamma}\mathbf{u}\cdot\mathbf{n}\,\operatorname{d}\\!s=0,$

is equivalent to

$\begin{cases}\textrm{there exists
}\mathbf{U}\in\left(L^{2}(\Omega)\right)^{2}\textrm{ such that }\\\
\mathrm{div}\,\mathbf{U}=0\textrm{ in }\Omega\textrm{ and
}\mathbf{U}\cdot\mathbf{n}=\mathbf{u}\cdot\mathbf{n}\textrm{ almost everywhere
on }\Gamma.\end{cases}$

∎

Define the variational formulation:

(5.5) $\begin{cases}\textrm{find
}\left(\mathbf{U},\mathbf{u}\right)\in\left(L^{2}(\Omega)\right)^{2}\times\left(H^{-\frac{1}{2}}\left(\Gamma\right)\cap
BV\left(\Gamma\right)\right)\\\ \textrm{and }p\in H^{1}(\Omega)\textrm{ that
solve}\\\
\int_{\Omega}\mathbf{U}\cdot\mathbf{V}\,\operatorname{d}\\!x+\int_{\Omega}\nabla
p\cdot\mathbf{V}\,\operatorname{d}\\!x+\sigma\int_{\Gamma}\widetilde{\mathbf{H}}\cdot\mathbf{v}\,\operatorname{d}\\!s=0\,,\\\
\int_{\Gamma}q\,\mathbf{u}\cdot\mathbf{n}\,\operatorname{d}\\!s=0\,,\\\
\textrm{ for any smooth test functions
}\mathbf{V}:\Omega\to\mathbb{R}^{2}\,,\mathbf{v}\in BV(\Gamma^{i})\\\ \textrm{
and }q:\Omega^{i}\to\mathbb{R}\textrm{ s.t.
}\mathbf{v}_{\Gamma}=0\,.\end{cases}$

######  Proposition 5.6.

The variational formulation (5.5) is the Euler equation for the minimization
problem (2.10) – (2.11).

###### Proof.

Define the lagrangian $\mathcal{L}_{\rm mod}$ by

$\mathcal{L}_{\rm
mod}(\mathbf{u},\mathbf{U},p)=\frac{1}{2}\int_{\Omega}|\mathbf{U}|^{2}\,\operatorname{d}\\!x+\int_{\Gamma}\left(\mathbf{u}\cdot\mathbf{n}\right)p\,\operatorname{d}\\!s-\int_{\Omega}\mathbf{U}\cdot\nabla
p\,\operatorname{d}\\!x+\sigma\mathcal{P}((\textrm{Id}+\Delta
t\,\mathbf{u})(\Omega))\,.$

Let $F(\mathbf{u})$ be the perimeter functional in the deformed configuration,
i.e.,

(5.6) $F(\mathbf{u})=\mathcal{P}((\textrm{Id}+\Delta
t\,\mathbf{u})(\Omega))=\int_{\Gamma}\left|(\mathbf{I}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}\right|\,\operatorname{d}\\!s\,.$

Then,

$F^{\prime}(\mathbf{u})(\mathbf{v})=\mathcal{P}^{\prime}\left((\textrm{Id}+\Delta
t\,\mathbf{u})(\Omega);\mathbf{v}\right)\,,$

and using Lemma 3.1, we compute the first order expansion of F at
$\mathbf{u}$.

######  Lemma 5.7.

It holds that

(5.7)
$\int_{\Gamma}\widetilde{\mathbf{H}}\cdot\mathbf{v}\,\operatorname{d}\\!s=\frac{1}{\Delta
t}F^{\prime}(\mathbf{u})(\mathbf{v})=\int_{\Gamma^{i}}\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}\,\operatorname{d}\\!s.$

The result then follows. ∎

######  Remark 5.8.

Note that the expression on the right-hand side of (5.7) absolves us from the
difficult task of computing $\widetilde{\mathbf{H}}$ directly. That being
said, our problem is still nonlinear through the dependence of $\mathbf{T}$ on
$\mathbf{u}$, see (3.3).

###  5.3 A Newton algorithm to solve (5.5)

The main difficulty in solving (5.5) is to find a method to handle the
nonlinear boundary term, (5.7).

We solve the problem (5.5) by a Newton method by seeking, for $\mathbf{u}$
given, a correction $\delta\mathbf{u}$ such that $(\delta\mathbf{u},p)$ is
solution of the linearized system

(5.8) $\begin{aligned}
\int_{\Omega}\left(\mathbf{u}+\delta\mathbf{u}\right)\cdot\mathbf{v}\,\operatorname{d}\\!x-\int_{\Omega}p(\nabla\cdot\mathbf{v})\,\operatorname{d}\\!x&-\int_{\Omega}\left(\nabla\cdot\left(\mathbf{u}+\delta\mathbf{u}\right)\right)q\,\operatorname{d}\\!x\\\
+&\frac{\sigma}{\Delta
t}F^{\prime}(\mathbf{u})(\mathbf{v})+\frac{\sigma}{\Delta
t}F^{\prime\prime}(\mathbf{u})\left(\delta\mathbf{u},\mathbf{v}\right)=0\,.\end{aligned},$

######  Remark 5.9.

The last term in Eq. (5.8) is a bilinear form. It is obtained by finding the
second order asymptotic expansion of $F$ at $\mathbf{u}$.

More precisely, we proceed iteratively. At each time step, it consists of
computing a sequence $(u^{k})_{k}$ where $\mathbf{u}^{0}=0$ and
$\mathbf{u}^{k+1}\in\left(H^{1}(\Omega)\right)^{2}$, $p\in L^{2}(\Omega)$ are
solutions to the following variational problem

(5.9)
$\displaystyle\int_{\Omega}\mathbf{u}^{k+1}\cdot\mathbf{v}\,\operatorname{d}\\!x$
$\displaystyle-\int_{\Omega}p(\nabla\cdot\mathbf{v})\,\operatorname{d}\\!x-\int_{\Omega}(\nabla\cdot\mathbf{u}^{k+1})q\,\operatorname{d}\\!x$
$\displaystyle+$ $\displaystyle\frac{\sigma}{\Delta
t}F^{\prime}(\mathbf{u}^{k})(\mathbf{v})+\frac{\sigma}{\Delta
t}F^{\prime\prime}(\mathbf{u}^{k})\left(\mathbf{u}^{k+1}-\mathbf{u}^{k},\mathbf{v}\right)=0$

for any arbitrary smooth test functions $\mathbf{v}$, $q$.

Assuming that this method converges, i.e.,
$\lim_{k\to\infty}|\mathbf{u}^{k+1}-\mathbf{u}^{k}|\to 0$, we shall denote by
$\mathbf{u}$ the limit of $(u^{k})_{k}$ for $k\to\infty$. The remaining
challenge is to compute $F^{\prime\prime}$.

#### Second order expansion of the deformed perimeter

To apply our Newton like method we must first perform an asymptotic expansion
of order two of the perimeter functional in the deformed configuration.

Let us compute the second order Taylor expansion of $F$ defined by (5.6).

Let $\operatorname{d}\\!s$, $\mathbf{\boldsymbol{\tau}}$ and $\mathbf{n}$
denote respectively the unit length, the tangent and the outward pointing unit
normal vectors in the current configuration $\Omega$, and
$\operatorname{d}\\!S$, $\mathbf{T}$ and $\mathbf{N}$ are respectively the
unit length, the tangent and the outward pointing unit normal vectors in the
deformed configuration $(\textrm{Id}+\Delta t\mathbf{u})(\Omega)$, ”pulled
back” in the coordinate system of the current configuration.

######  Lemma 5.10.

It holds that

(5.10) $F^{\prime\prime}(\mathbf{u})(\delta\mathbf{u},\mathbf{v})=\Delta
t^{2}\int_{\Gamma}\frac{\left(\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{N}\right)\left(\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{N}\right)}{\operatorname{d}\\!S}\operatorname{d}\\!s\,.$

###### Proof.

We see that

(5.11) $\displaystyle F(\mathbf{u}+\delta\mathbf{u})=$ $\displaystyle
F(\mathbf{u})+\Delta
t\int_{\Gamma}\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}\operatorname{d}\\!s$
$\displaystyle+\frac{\Delta
t^{2}}{2}\int_{\Gamma}\frac{\left(\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{N}\right)^{2}}{\operatorname{d}\\!S}\operatorname{d}\\!s+O(||\nabla\delta\mathbf{u}||^{3})\,.$

∎

#### First linearized problem

Substituting (5.7) and (5.11) into our Newton like method, (5.9), gives

(5.12)
$\displaystyle\int_{\Omega}\mathbf{u}^{k+1}\cdot\mathbf{v}\,\operatorname{d}\\!x-\int_{\Omega}p(\nabla\cdot\mathbf{v})\,\operatorname{d}\\!x-\int_{\Omega}(\nabla\cdot\mathbf{u}^{k+1})q\,\operatorname{d}\\!x$
$\displaystyle+\sigma\int_{\Gamma}\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}^{k}\operatorname{d}\\!s+\sigma\Delta
t\int_{\Gamma}\frac{(\nabla\delta\mathbf{u}^{k+1}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{N}^{k})(\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{N}^{k})}{\operatorname{d}\\!S^{k}}\operatorname{d}\\!s=0$

where $\delta\mathbf{u}^{k+1}:=(\mathbf{u}^{k+1}-\mathbf{u}^{k})$ and

$\operatorname{d}\\!S^{k}=|(\textrm{I}+\Delta
t\,\nabla\mathbf{u}^{k})\mathbf{\boldsymbol{\tau}}|,\quad\mathbf{T}^{k}=\frac{(\textrm{I}+\Delta
t\,\nabla\mathbf{u}^{k})\mathbf{\boldsymbol{\tau}}}{|(\textrm{I}+\Delta
t\,\nabla\mathbf{u}^{k})\mathbf{\boldsymbol{\tau}}|},\quad\mathbf{N}^{k}=\mathbf{T}^{k^{\perp}}\,.$

######  Remark 5.11.

We emphasize the fact that (5.12) does not contain second-order spatial
derivatives, meaning that the curvature vector has been transformed into a
term involving only the first spatial derivatives (hence less regularity is
needed).

This system in the variables ($\mathbf{u}^{k+1}$, $p$) does not necessarily
admit a solution (let alone unique) due to the possible lack of coercivity of
the bilinear form. Moreover, even if solutions to (5.12) exist, the
convergence of this iterative Newton method, (5.9) with (5.12), is not
granted. The straighforward remedy is to replace the bilinear term in (5.12)
by a coercive one, which can be done in various ways. The classical scheme is
to compute the eigenvectors and eigenvalues of the matrix and to construct the
modified term using only the eigenspaces associated with the positive
eignevalues.

For simplicity and robustness of our algorithm, we prefer to predetermine a
closed formula for the modified matrix which is coercive.

#### Modified problem

The main obstacle in our current method is that the bilinear form in (5.12) is
not positive definite. Therefore, we want to define a modified problem of
(5.9) – (5.12) with a positive-definite matrix.

######  Lemma 5.12.

It holds that

(5.13) $F(\mathbf{u}+\delta\mathbf{u})\leq F(\mathbf{u})+\Delta
t\int_{\Gamma}\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}\operatorname{d}\\!s+\frac{\Delta
t^{2}}{2}\int_{\Gamma}\frac{(\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}})\cdot(\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}})}{\operatorname{d}\\!S}\operatorname{d}\\!s\,.$

###### Proof.

From the inequality $\sqrt{1+x}\leq 1+x/2$ it follows that

(5.14) $\displaystyle|(\textrm{I}+\Delta
t\,\nabla(\mathbf{u}+\delta\mathbf{u}))\mathbf{\boldsymbol{\tau}}|$
$\displaystyle\quad=\Big{(}|(\textrm{I}+\Delta
t\,\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|^{2}+2\Delta
t\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot(\mathbf{\boldsymbol{\tau}}+\Delta
t\,\nabla\mathbf{u}\,\mathbf{\boldsymbol{\tau}})+\Delta
t^{2}|\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}|^{2}\Big{)}^{1/2}$
$\displaystyle\quad\leq\,|(\textrm{I}+\Delta
t\,\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}}|+\Delta
t\,\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}+\frac{\Delta
t^{2}|\nabla\delta\mathbf{u}\,\mathbf{\boldsymbol{\tau}}|^{2}}{2\operatorname{d}\\!S}\,,$

where $\operatorname{d}\\!S$ and $\mathbf{T}$ are defined by (5.12). ∎

We stress that here, unlike (5.11), the bilinear form on the right-hand side
of (5.13) is positive definite. We choose to adopt this term in our modified
Newton like method.

To recapitulate, in each time step we omit the index $i$ and compute a
sequence $(u^{k})_{k}$, where $\mathbf{u}^{0}=0$ and $\mathbf{u}^{k+1}\in
H^{1}(\Omega)^{2}$, $p\in L^{2}(\Omega)$ are solutions of the following
variational problem

(5.15)
$\displaystyle\int_{\Omega}\mathbf{u}^{k+1}\cdot\mathbf{v}\,\operatorname{d}\\!x$
$\displaystyle-\int_{\Omega}p(\nabla\cdot\mathbf{v})\,\operatorname{d}\\!x-\int_{\Omega}(\nabla\cdot\mathbf{u}^{k+1})q\,\operatorname{d}\\!x$
$\displaystyle+\sigma\int_{\Gamma}\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}}\cdot\mathbf{T}^{k}\operatorname{d}\\!s+\sigma\Delta
t\int_{\Gamma}\frac{(\nabla\delta\mathbf{u}^{k+1}\,\mathbf{\boldsymbol{\tau}})\cdot(\nabla\mathbf{v}\,\mathbf{\boldsymbol{\tau}})}{\operatorname{d}\\!S^{k}}\operatorname{d}\\!s=0$

for all arbitrary smooth test functions $\mathbf{v}:\Omega\to\mathbb{R}^{2}$,
$q:\Omega\to\mathbb{R}$, and where

$\operatorname{d}\\!S^{k}=|(\textrm{I}+\Delta
t\,\nabla\mathbf{u}^{k})\mathbf{\boldsymbol{\tau}}|,\quad\mathbf{T}^{k}=\frac{(\textrm{I}+\Delta
t\,\nabla\mathbf{u}^{k})\mathbf{\boldsymbol{\tau}}}{|(\textrm{I}+\Delta
t\,\nabla\mathbf{u}^{k})\mathbf{\boldsymbol{\tau}}|}$

In our algorithm, the method is applied recursively until the stopping
criteria based on the computation of the global residual is satisfied. We set
the Newton tolerance to $10^{-5}$ in our computations, i.e.,

$\int_{\Omega}\left|\delta\mathbf{u}^{k+1}\right|^{2}=\int_{\Omega}\left|\mathbf{u}^{k+1}-\mathbf{u}^{k}\right|^{2}<10^{-5}$

When this condition is satisfied, we take ($\mathbf{u}^{k+1}$, $p$) as the
approximate solution for ($\mathbf{u}^{i+1}$, $p^{i+1}$) in the time-discrete
PDE problem. All that is left is to propagate the domain via Eq. (2.5).

##  6 Variational formulation associated with (2.12) – (2.13)

Let $\alpha>0$, define

(6.1)
$J_{\alpha}(\Omega,\mathbf{u})=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\,\operatorname{d}\\!x+\frac{\sigma}{\Delta
t}\mathcal{P}(\varphi(\Omega))+\frac{\alpha}{2}\int_{\Omega}|\nabla\wedge\mathbf{u}|^{2}\,\operatorname{d}\\!x.$

######  Proposition 6.1.

Let $\Omega$ be a $C^{2}$ open set of $\mathbb{R}^{2}$. Any minimizer
$\mathbf{u}$ of $J_{\alpha}(\Omega,\cdot)$ with respect to $\mathbf{u}$ over
the set of zero divergence fields satisfies

(6.2)
$\left\\{\begin{array}[]{ll}\mathbf{u}-\alpha\nabla^{\perp}(\nabla\wedge\mathbf{u})+\nabla
p=0&\text{ in }\Omega,\\\ \nabla\cdot\mathbf{u}=0&\text{ in }\Omega,\\\
p\mathbf{n}=\sigma\widetilde{\mathbf{H}_{\varphi}}+\alpha(\nabla\wedge\mathbf{u})\mathbf{\boldsymbol{\tau}}&\text{
on }\partial\Omega.\end{array}\right.$

###### Proof.

For any test function $v$, we have

$\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\alpha\int_{\Omega}(\nabla\wedge\mathbf{u})(\nabla\wedge\mathbf{v})\,\operatorname{d}\\!x+\sigma\int_{\partial\Omega}\widetilde{\mathbf{H}_{\varphi}}\cdot\mathbf{v}\,\operatorname{d}\\!s=-\int_{\Omega}\nabla
p\cdot\mathbf{v}\,\operatorname{d}\\!x+\int_{\partial\Omega}p(\mathbf{v}\cdot\mathbf{n})\,\operatorname{d}\\!s.$

By integration by parts, we see that

$\int_{\Omega}(\nabla\wedge\mathbf{u})(\nabla\wedge\mathbf{v})\,\operatorname{d}\\!x=-\int_{\Omega}\nabla^{\perp}(\nabla\wedge\mathbf{u})\cdot\mathbf{v}\,\operatorname{d}\\!x+\int_{\partial\Omega}(\nabla\wedge\mathbf{u})\mathbf{n}\wedge\mathbf{v}\,\operatorname{d}\\!s.$

∎

Define the variational formulation:

(6.3) $\begin{cases}\textrm{find }\mathbf{u}\in L^{2}(\Omega)\textrm{and }p\in
H^{1}(\Omega)\textrm{ that solve}\\\
\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\alpha\int_{\Omega}(\nabla\wedge\mathbf{u})\cdot(\nabla\wedge\mathbf{v})\,\operatorname{d}\\!x+\int_{\Omega}\nabla
p\cdot\mathbf{v}\,\operatorname{d}\\!x+\sigma\int_{\Gamma}\widetilde{\mathbf{H}}\cdot\mathbf{v}\,\operatorname{d}\\!s=0\,,\\\
\int_{\Gamma}q\,\mathbf{u}\cdot\mathbf{n}\,\operatorname{d}\\!s=0\,,\\\
\textrm{ for any smooth test functions
}\mathbf{v}:\Omega\to\mathbb{R}^{2}\textrm{ and
}q:\Omega^{i}\to\mathbb{R}\,.\end{cases}$

######  Proposition 6.2.

The variational formulation (6.3) is the Euler equation for the minimization
problem of (2.12) – (2.13).

###### Proof.

The functional $J_{\alpha}$ is Gâteaux differentiable with respect to
$\mathbf{u}$. Denoting by $J_{\alpha}^{\prime}$ its Gâteaux derivative and
using lemma 3.1, we obtain for all
$\mathbf{v}\in\left(H^{1}(\Omega)\right)^{2}$,

(6.4) $\displaystyle\langle
J_{\alpha}^{\prime}(\Omega,\mathbf{u});\mathbf{v}\rangle$ $\displaystyle=$
$\displaystyle\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\sigma\int_{\partial\varphi(\Omega)}\mathbf{H}_{\varphi}\cdot(\mathbf{v}\circ\varphi^{-1})\,\operatorname{d}\\!s_{\varphi}$
$\displaystyle\qquad+\alpha\int_{\Omega}(\nabla\wedge\mathbf{u})\cdot(\nabla\wedge\mathbf{v})\,\operatorname{d}\\!x$
$\displaystyle=$
$\displaystyle\int_{\Omega}\mathbf{u}\cdot\mathbf{v}\,\operatorname{d}\\!x+\sigma\int_{\partial\Omega}\widetilde{\mathbf{H}_{\varphi}}\cdot\mathbf{v}\,\operatorname{d}\\!s$
$\displaystyle\qquad+\alpha\int_{\Omega}(\nabla\wedge\mathbf{u})\cdot(\nabla\wedge\mathbf{v})\,\operatorname{d}\\!x\,.$

Define the lagrangian $\mathcal{L}_{\alpha}^{*}$ by

$\mathcal{L}^{*}(\mathbf{u},p)=\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\,\operatorname{d}\\!x+\frac{\alpha}{2}\int_{\Omega}|\nabla\wedge\mathbf{u}|^{2}\,\operatorname{d}\\!x-\int_{\Omega}\left(\mathrm{div}\,\mathbf{u}\right)p\,\operatorname{d}\\!x+\sigma\int_{\partial\Omega}\widetilde{\mathbf{H}_{\varphi}}\cdot\mathbf{u}\,\operatorname{d}\\!s\,.$

We proceed as in the proof of Proposition 4.2 by applying de Rahm’s Theorem
A.1 to the linear form
$\langle\mathcal{L}_{\alpha}^{*^{\prime}}(\mathbf{u},p);\cdot\rangle$, to
deduce that if $\mathbf{u}$ is a minimizer of $J(\Omega,\cdot)$ over the set
of divergence-free fields, then there exists $p$ such that for all Lipschitz
test function $\mathbf{v}$ and $q$,

(6.5)
$\langle\mathcal{L}_{\alpha}^{*^{\prime}}(\mathbf{u},p);(\mathbf{v},q)\rangle=\int_{\Omega}p\nabla\cdot\mathbf{v}\,\operatorname{d}\\!x\,.$

∎

######  Remark 6.3.

Note that a priori, it does not seem obvious that the obtained scheme is
consistent, due to the presence of the term $\alpha\nabla\wedge\mathbf{u}$. By
taking the rotational of the first equation of the system (6.2) and the dot
product with the tangent $\tau$ at the boundary for the last equation of
(6.2), we get

$\left\\{\begin{array}[]{ll}\nabla\wedge\mathbf{u}-\Delta(\alpha\nabla\wedge\mathbf{u})=0&\text{
in }\Omega\\\
\alpha(\nabla\wedge\mathbf{u})=-\sigma\widetilde{\mathbf{H}_{\varphi}}\cdot\mathbf{\boldsymbol{\tau}}&\text{
on }\partial\Omega\end{array}\right.$

The boundary term is of order $\Delta t$. Indeed, if $\mathbf{u}$ is smooth,
$C^{2}$ e.g., we have

$\widetilde{\mathbf{H}_{\varphi}}=\nabla((\textrm{Id}+\Delta
t\nabla\mathbf{u})\mathbf{\boldsymbol{\tau}})\cdot\mathbf{\boldsymbol{\tau}}=\Delta
t\partial^{2}_{s}\mathbf{u}\cdot\mathbf{\boldsymbol{\tau}}.$

In particular we deduce

$\alpha\|\nabla\wedge\mathbf{u}\|_{H^{1}(\Omega)}\leq C(\Omega)\,\Delta
t\,\sigma\|\partial_{s}^{2}\mathbf{u}\cdot\mathbf{\boldsymbol{\tau}}\|_{H^{1/2}(\partial\Omega)}.$

The stabilization term therefore introduces an error of order $\Delta t$. Note
that the introduced consistency error is independent of $\alpha$.

##  7 Variational formulation associated with (2.15) – (2.16)

In order to alleviate the problem of the existence of the time-discretized
scheme proposed in the previous section, another option consists in replacing
the linearized incompressibility constraint $\nabla\cdot\mathbf{u}=0$ by the
nonlinear constraint $\det(\nabla\varphi)=1$.

######  Proposition 7.1.

Let $\mathbf{u}$ be a minimizer of $J$ under the constraint

$\det(\nabla\varphi)=1\text{ in }\Omega,$

with $\varphi=\textrm{Id}+\Delta t\mathbf{u}$. Then, there exists a Lagrange
multiplier $p$ such that

$\left\\{\begin{array}[]{rcl}\mathbf{u}+(\operatorname{Cof}\nabla\varphi)\nabla
p=0&\text{ in }&\Omega\\\ \det(\nabla\varphi)=1&\text{ in }&\Omega\\\
p\mathbf{n}=\mathbf{H}_{\varphi}\circ\varphi&\text{ on
}&\partial\Omega.\end{array}\right.$

###### Proof.

In a classical way, to take into account the constraint of nonlinear
incompressibility, one introduces the Lagrangian

(7.1) $\mathcal{L}(\Omega;\mathbf{u},p)=J(\Omega;\mathbf{u})-\frac{1}{\Delta
t}\int_{\Omega}(\det(\nabla\varphi)-1)p\,\operatorname{d}\\!x.$

Define $V=\mathbf{v}\circ\varphi^{-1}$ and $P=p\circ\varphi^{-1}$. We first
see that

$\nabla_{\mathbf{y}}V=\nabla_{\mathbf{x}}\mathbf{v}(\nabla_{\mathbf{x}}\varphi)^{-1}\,,$

with the notation $\mathbf{y}=\varphi(\mathbf{x})$.

Moreover, we compute

$\displaystyle\int_{\Omega}\det(\nabla\varphi+\Delta
t\nabla\mathbf{v})p\,\operatorname{d}\\!x$ $\displaystyle=$
$\displaystyle\int_{\varphi(\Omega)}\det(\nabla_{\mathbf{x}}\varphi)^{-1}\det(\nabla_{\mathbf{x}}\varphi+\Delta
t\nabla_{\mathbf{x}}\mathbf{v})P\,\operatorname{d}\\!x$ $\displaystyle=$
$\displaystyle\int_{\varphi(\Omega)}\det(\textrm{Id}+\Delta
t\nabla_{\mathbf{y}}V)P\,\operatorname{d}\\!x$ $\displaystyle=$
$\displaystyle\int_{\varphi(\Omega)}(1+\Delta t\nabla_{\mathbf{y}}\cdot
V)P\,\operatorname{d}\\!x+o(\mathbf{v}).$

Therefore, the Gâteaux derivation of the last term of (7.1) is

$-\frac{1}{\Delta t}\int_{\varphi(\Omega)}(\nabla_{\mathbf{y}}\cdot
V)P\,\operatorname{d}\\!x,$

and integrating by parts we get

$\frac{1}{\Delta t}\int_{\varphi(\Omega)}(\nabla_{\mathbf{y}}\cdot
V)P\,\operatorname{d}\\!x=-\int_{\varphi(\Omega)}V\cdot\nabla_{\mathbf{y}}P\,\operatorname{d}\\!x+\int_{\partial\varphi(\Omega)}V\cdot\mathbf{n}_{\varphi}P\operatorname{d}\\!s.$

Finally, recalling that $P\circ\varphi=p$, it follows that

$\nabla_{\mathbf{y}}P\nabla\varphi=\nabla_{\mathbf{x}}p,$

hence

$\nabla\varphi^{T}\nabla_{\mathbf{y}}P=\nabla_{\mathbf{x}}p\,$

that is

$\nabla_{\mathbf{y}}P=\nabla\varphi^{-T}\nabla_{\mathbf{x}}p.$

We thus deduce that the Gâteaux derivation of the last term of (7.1) is

$\displaystyle-\frac{1}{\Delta
t}\int_{\varphi(\Omega)}(\nabla_{\mathbf{y}}\cdot V)P\,\operatorname{d}\\!x$
$\displaystyle=$
$\displaystyle\int_{\Omega}\mathbf{v}\cdot\nabla\varphi^{-T}\nabla_{\mathbf{x}}p\det(\nabla\varphi)\,\operatorname{d}\\!x-\int_{\partial\varphi(\Omega)}V\cdot\mathbf{n}_{\varphi}P\,\operatorname{d}\\!s$
$\displaystyle=$
$\displaystyle\int_{\Omega}\mathbf{v}\cdot(\operatorname{Cof}\nabla\varphi)\nabla_{x}p\,\operatorname{d}\\!x-\int_{\partial\varphi(\Omega)}V\cdot\mathbf{n}_{\varphi}P\,\operatorname{d}\\!s\,.$

If $\mathbf{u}$ is a minimizer of $J$ with respect to $\mathbf{u}$ under the
constraint $\det(\nabla\varphi)=1$, then $\mathbf{u}$ is a critical point of
the Lagrangian. Using compactly supported test functions in $\Omega$, we
deduce that

$\mathbf{u}+(\operatorname{Cof}\nabla\varphi)\nabla p=0\quad\text{ in
}\Omega\,,$

and by using any test functions in $\Omega$ we obtain the boundary condition

$p\mathbf{n}_{\varphi}=\mathbf{H}_{\varphi}\quad\text{ on
}\partial\varphi(\Omega)\,,$

that is

$p\mathbf{n}_{\varphi}\circ\varphi=\mathbf{H}_{\varphi}\circ\varphi\quad\text{
on }\partial\Omega\,.$

∎

##  8 Verification: simulation-theory comparisons

###  8.1 Space discretization

Each subdomain is covered by a regular triangulation $\mathcal{T}_{h}$, with
maximum mesh size $h$, and such that it is globally a conforming triangulation
of $\Omega$, i.e., $\mathcal{T}_{h}$ contains a piecewise affine approximation
$\Gamma_{h}$ of the interface $\Gamma$.

In a classical manner, we approximate each component of the velocity in each
element $K\in\mathcal{T}_{h}$ by a polynomial of degree one enriched with a
“bubble” function (a polynomial of degree 3 defined as the product of the
barycentric coordinates in $K$ and vanishing on the faces of $K$) and the
pressure in each element by a polynomial of degree one. Both approximations
are continuous across the element faces except for the pressure at the
interface $\Gamma_{h}$. Hence, we consider the following discretizations of
the spaces $X:=\left(H^{1}(\Omega)\right)^{2}$ and $M:=L^{2}(\Omega)$:

$X_{h}=\\{v_{h}\in\mathcal{C}^{0}(\bar{\Omega})^{2}\,;\ \forall
K\in\mathcal{T}_{h}\,,\ v_{h|K}\in\left(\mathbb{P}_{1}+b_{K}\right)^{2}\\}\cap
X\,,$

and

$M_{h}=\\{q_{h}\in\mathcal{C}^{0}(\bar{\Omega})\,;\ \forall
K\in\mathcal{T}_{h}\,,\ q_{h|K}\in\mathbb{P}_{1}\\}\cap M\,,$

The numerical experiments presented in the following sections were conducted
using FreeFem++.

In this section, we validate our FEM simulation by comparing numerical
experiments with theoretical predictions.

To simulate the passive droplet problem (using our code in FreeFem++) one must
specify the tension $\sigma$, the numerical time step $\Delta t$ and the
initial finite element domain (meaning the triangulation mesh
$\mathcal{T}_{h}$). Let us first comment on the choice of the time step size.
The convergence of the Newton algorithm is indeed quite sensitive to $\Delta
t$. Values which are too large may lead to starting iterations in our
recursive method which are far from the expected solutions. There exist
several strategies for improving the choice of $\Delta t$ but we leave this
topic beyond the scope of the present study. In each numerical experiment, we
choose a value of $\Delta t$ that is 2–4 orders of magnitude smaller than the
physical timescale. Building the finite element domain is done by inputting an
explicit counter-clockwise paremeterization of the closed interface. Here, we
use the polar parameterization
$\left\\{x(\theta),y(\theta)\right\\}=\left\\{R^{0}(\theta)\cos(\theta),R^{0}(\theta)\sin(\theta)\right\\}$
for $\theta\in(0,2\pi)$. We define $R^{0}(\theta)=1+\delta R^{0}(\theta)$,
where we write $\delta R^{0}(\theta)$ in terms of Fourier modes

(8.1) $\delta R^{0}(\theta)=\sum_{m}\left(\delta
R_{\text{c}m}^{0}\cos(m\theta)+\delta R_{\text{s}m}^{0}\sin(m\theta)\right)$

Then, by specifying the number (or density) of vertices along the
parameterized boundary, FreeFem++ automatically generates the internal
triangulation mesh (see Fig. 2a,b). Note that one could also define an
adaptive (non-uniform) mesh, as in Fig. 2c, which is designed to have a finer
definition of vertices along the boundary. The motivation behind such a mesh
is to improve the resolution of the shape itself while not ’hyper-meshing’ the
bulk and thereby drastically increasing computation time (of order $h^{-2}$).

Figure 2: Building different finite element domains for the same interface.
The parameterization $R^{0}(\theta)=1+\delta R^{0}(\theta)$ (red, Eq. (8.1))
is used to define the affine approximation of the initial interface,
$\Gamma_{h}^{0}$ (black). In turn, FreeFem++ uses $\Gamma_{h}^{0}$ to generate
the internal tiangulation mesh $\mathcal{T}_{h}^{0}$ (gray). The density of
the mesh is controlled by specifying the number of vertices on the interface.
This number is increased from (a) to (b). The adaptive mesh (c) is constructed
with two starting interfaces: one (dense) interface on the boundary and
another (coarse) interface in the bulk. These define two triangulation meshes
(a ”disk” and a ”ring”) which are then joined to form a single connected mesh.

In each time step, we follow the algorithm outlined at the end of Section 5.3
(Eq. (5.15)). The simulation data is saved at some fixed interval of time
iterations (of order 10-100, depending on $\Delta t$ and the duration of the
simulation). In each imported frame, we are able to reconstruct the finite-
element domain and the interpolation functions for $\mathbf{u}$ and $p$. In
addition, we use the interface vertices to construct a polar piece-wise
interpolation function of the boundary, $R^{\text{sim}}(\theta)$.

Our objective is to square the simulation results with known characteristics
of the passive droplet, namely:

  * •

Conservation of the droplet area, $\dot{A}=0$.

  * •

External force balance, $\mathbf{u}_{\text{cm}}=0$.

  * •

Morphological relaxation dynamics of linear shape perturbations.

######  Remark 8.1.

There are several ways of computing $\mathbf{u}_{\text{cm}}$. We define

$\mathbf{u}_{\text{cm}}=A^{-1}\int_{\partial\Omega}\mathbf{x}(\mathbf{u}\cdot\mathbf{n})\operatorname{d}\\!l.$

######  Remark 8.2.

To find the numerical growth rate of each normal mode, we first decompose
$R^{\text{sim}}(\theta)$ into Fourier components,

(8.2) $\displaystyle\delta
R^{\text{sim}}_{\text{c}m}=\frac{1}{\pi}\int_{0}^{2\pi}\left(R^{\text{sim}}(\theta)-1\right)\cos(m\theta)d\theta$
$\displaystyle\delta
R^{\text{sim}}_{\text{s}m}=\frac{1}{\pi}\int_{0}^{2\pi}\left(R^{\text{sim}}(\theta)-1\right)\sin(m\theta)d\theta$

The growth rate of each such component is then obtained by fitting $\delta
R_{\text{c}m}^{\text{sim}}(t)$ and $\delta R_{\text{s}m}^{\text{sim}}(t)$ to
an exponential function $\varepsilon e^{s_{m}^{\text{N}}t}$. The idea is to
compare the fitted $s_{m}^{\text{N}}$ with the classical cubic dispersion
relation, $s_{m}=-\sigma m(m^{2}-1)$.

Figure 3: Decaying shape perturbation (mode $m=2$): We set $\delta
R_{\text{c}2}^{0}=.05$ in Eq. (8.1), fixed $\sigma=0.5$, and chose $\Delta
t=0.0005\tau_{2}$, where $\tau_{2}=|s_{2}|^{-1}=1/3$. At the top we present
simulation snapshots demonstrating the decay of the shape perturbation. The
density plot in the bulk represents the pressure $p$ (from low in white to
high in blue). (a) Time series for the domain area $A(t)$. (b) Time series for
the absolute center of mass velocity $|\mathbf{u}_{\text{cm}}(t)|$. (c) The
predicted linear behaviour of the perturbed mode, $\log|\delta
R_{\text{c}2}^{\text{lsa}}(t)|=\log|\delta R_{\text{c}2}^{0}|-3t$ (continuous
orange line) vs. the simulation time series, $\log|\delta
R_{\text{c}2}^{\text{sim}}(t)|$ (black, computed via Eq. (8.2)). We find that
the fitted growth rate, $s_{2}^{\text{N}}\simeq-2.995$, compares well with the
classical linear-stability growth rate, $s_{2}=-3$ (giving a deviation of
$(s_{2}^{\text{N}}-s_{2})/s_{2}\sim 2\times 10^{-3}$). Figure 4: Decaying
shape perturbation (mode $m=3$): We set $\delta R_{\text{c}3}^{0}=.05$ in Eq.
(8.1), fixed $\sigma=0.5$, and chose $\Delta t=0.0005\tau_{3}$, where
$\tau_{3}=|s_{3}|^{-1}=1/12$. At the top we present simulation snapshots
demonstrating the decay of the shape perturbation. The density plot in the
bulk represents the pressure $p$ (from low in white to high in blue). (a) Time
series for the domain area $A(t)$. (b) Time series for the absolute center of
mass velocity $|\mathbf{u}_{\text{cm}}(t)|$. (c) The predicted linear
behaviour of the perturbed mode, $\log|\delta
R_{\text{c}3}^{\text{lsa}}(t)|=\log|\delta R_{\text{c}3}^{0}|-12t$ (continuous
green line) vs. the simulation time series, $\log|\delta
R_{\text{c}3}^{\text{sim}}(t)|$ (black, computed via Eq. (8.2)). We find that
the fitted growth rate, $s_{3}^{\text{N}}\simeq-11.996$, compares well with
the classical linear-stability growth rate, $s_{3}=-12$ (giving a deviation of
$(s_{2}^{\text{N}}-s_{2})/s_{2}\sim 4\times 10^{-4}$). Figure 5: Decaying
shape perturbations (mixed modes $m=2$–$5$): We set $\delta
R_{\text{c}2}=.03$, $\delta R_{\text{s}3}^{0}=-.03$, $\delta
R_{\text{s}4}=.03$, and $\delta R_{\text{c}5}=-.03$ in Eq. (8.1), fixed
$\sigma=0.5$, and chose $\Delta t=0.005\tau_{5}$, where
$\tau_{5}=|s_{5}|^{-1}=1/120$. At the top we present simulation snapshots
demonstrating the sequential decay of the perturbed shape modes. The density
plot in the bulk represents the pressure $p$ (from low in white to high in
blue). (a) Time series for the droplet area $A(t)$. (b) Time series for the
absolute center of mass velocity $|\mathbf{u}_{\text{cm}}(t)|$. (c) The
predicted linear behaviour of each perturbed mode, $\log|\delta
R_{\text{c,s}m}^{\text{lsa}}(t)|=\log|R_{\text{c,s}m}^{0}|-s_{m}t$ (continuous
colored lines) vs. the simulation time series, $\log|\delta
R_{\text{c,s}m}^{\text{sim}}(t)|$ (darker colors, computed via Eq. (8.2)). We
find that the fitted growth rates, ($s_{2}^{\text{N}}\simeq-2.99$,
$s_{3}^{\text{N}}\simeq-11.91$, $s_{4}^{\text{N}}\simeq-29.83$,
$s_{5}^{\text{N}}\simeq-59.50$) are all in good quantitative agreement with
the classical linear-stability growth rates ($s_{2}=-3$, $s_{3}=-12$,
$s_{4}=-30$, $s_{5}=-60$).

Results of three simulations are represented in Figs. 3 – 5. In Figs. 3 and 4
we introduced an initial shape perturbation strictly in one Fourier mode
($m=2$ and $m=3$, respectively), whereas in Fig. 5 we introduced a
superposition of small perturbations in $m=2$–$5$. For the time step and mesh
density chosen, we find numerical deviations in $A(t)$ and
$|\mathbf{u}_{\text{cm}}(t)|$ as low as order $10^{-7}$. Moreover, the fitted
numerical growth rates of the perturbed Fourier modes are also in good
quantitative agreement with the classical dispersion relation $s_{m}$ (see
details in figure captions). We stress that small deviations in the fitted
growth rates may also arise from nonlinear effects which have been neglected
in the calculation of $s_{m}$. As expected, we found through further
experimentation with the numerics that precision is gained by decreasing
$\Delta t$, increasing the overall mesh density and/or decreasing the initial
perturbation amplitudes.

Figure 6: Convergence time for the different schemes

##  Appendix A De Rham’s Theorem

In the setting of partial differential equations, the De Rham’s theorem solves
an over-determined system of linear partial differential equations of order
one.

We recall de Rham’s theorem in the case of homogeneous flows (or currents) of
dimension one on a Euclidean space:

######  Theorem A.1 (de Rahm).

Let $L:L^{2}(\Omega)^{2}\to\mathbb{R}$ be a continuous linear form such that
$L(\mathbf{v})=0$ for all $\mathbf{v}\in L^{2}(\Omega)^{2}$ with
$\nabla\cdot\mathbf{v}=0$. Then, there exists $q\in H_{0}^{1}(\Omega)$ such
that for all $\mathbf{v}\in L^{2}(\Omega)^{2}$

$L(\mathbf{v})=\int_{\Omega}\nabla
q\cdot\mathbf{v}\operatorname{d}\\!\mathbf{x}\,.$

##  Appendix B Linear stability analysis of the Hele-Shaw model

The Hele-Shaw model writes

(B.1) $-\nabla\cdot u=\Delta p=0\qquad\textrm{ in }\Omega(t)\,,$

where $\Omega(t)$ is the domain occupied by the fluid and $u=-\nabla p$ is the
fluid velocity.

The dynamic boundary condition (or normal force balance) is the given by the
Young-Laplace pressure drop

(B.2) $p=\sigma\kappa\qquad\textrm{ on }\partial\Omega(t)\,,$

where $p$ is the fluid pressure (minus a constant), $\sigma$ is the surface
tension and $\kappa$ the local curvature. The free-boundary evolves with the
kinematic condition, which states that the normal velocity of the sharp
interface equals the normal velocity of the fluid $V_{n}=u\cdot\mathbf{n}$,
recalling that $u=-\nabla p$, in terms of the pressure this gives

(B.3) $V_{n}=-\nabla p\cdot\mathbf{n}\qquad\textrm{ on }\partial\Omega(t)\,,$

where $\mathbf{n}$ is the unit vector pointing outward.

We start with the stationary solution of a disk
$\Omega_{0}=\\{(r,\theta)|r\leq R_{0}\\}$. The curvature across the boundary
is $\kappa=1/R_{0}$ and the solution of the pressure is simply
$p_{0}=\sigma/R_{0}$ in $\Omega_{0}$.

We perturb the edge of the domain so that it is defined in terms of the polar
angle

$R(\theta,t)=R_{0}+\varepsilon R_{1}(\theta,t)\,.$

We want to analyze how the perturbation $R_{1}(\theta,t)$ evolves in time; in
particular we consider perturbations of the form
$R_{1}(\theta,t)=R_{m}(t)\cos(m\theta)$, where $R_{m}(t)$ is the amplitude and
$m$ is the wave number. If, for a given $m$, $R_{m}(t)$ grows in time, that
particular wave number is unstable; if $R_{m}(t)$ decreases, that wave number
is stable.

Since $\varepsilon$ is small, we will neglect all the terms which are
proportional to $\varepsilon^{n}$, $n>1$. Assume that the pressure can be
expanded as follows:

$p(r,\theta,t)=p_{0}+\varepsilon p_{1}(r,\theta,t)\,,$

with $p_{1}(r,\theta,t)=p_{m}(r,t)\cos(m\theta)$.

Since $p_{0}$ satisfies Laplace equation, $p_{1}$ will satisfy the Laplace
equation $\Delta p_{1}=0$. In polar coordinates, this writes

(B.4)
$\left(\partial_{r}^{2}+r^{-1}\partial_{r}-r^{-2}m^{2}\right)p_{m}(r,t)=0\,,$

which is solved by

$p_{m}(r,t)=A_{m}(t)r^{m}+B_{m}(t)r^{-m}\,.$

Discarding singularities at $r=0$ we set $B_{m}(t)=0$.

Concerning the curvature a direct computation gives that

$\kappa=\frac{1}{R_{0}}+\varepsilon\kappa_{1}\,,$

with

(B.5)
$\kappa_{1}=-\frac{R_{1}+\partial_{\theta}^{2}R_{1}}{R_{0}^{2}}=\frac{m^{2}-1}{R_{0}^{2}}R_{m}(t)\cos(m\theta)\,.$

From the dynamic boundary condition (B.2) it follows that

$(p_{0}+\varepsilon p_{1})_{|R_{0}(1+\varepsilon
R_{1})}=\sigma\kappa_{|R_{0}(1+\varepsilon
R_{1})}=\sigma\left(\frac{1}{R_{0}}+\varepsilon\frac{m^{2}-1}{R_{0}^{2}}R_{m}(t)\cos(m\theta)\right)\,,$

hence

(B.6) $p_{1|R_{0}(1+\varepsilon
R_{1})}=\sigma\frac{m^{2}-1}{R_{0}^{2}}R_{m}(t)\cos(m\theta)\,.$

On the other hand, neglecting all the terms which are proportional to
$\varepsilon^{n}$, $n>1$, we obtain

$p_{1|R_{0}(1+\varepsilon
R_{1})}=p_{m}(r=R_{0},t)\cos(m\theta)=A_{m}(t)R_{0}^{m}\,,$

hence

$A_{m}(t)=\frac{\sigma(m^{2}-1)}{R_{0}^{m+2}}R_{m}(t)\cos(m\theta)\,.$

We can now analyze the time evolution of the interface. From the kinematic
condition (B.3) it follows that
$\varepsilon\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}R_{1}=-\varepsilon\partial_{r}p_{1}$
on the boundary. Hence,

$\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}R_{1}=-\partial_{r}\left(A(t)r^{m}\right)_{|r=R_{0}}=-\sigma
R_{0}^{-3}m(m^{2}-1)R_{m}(t)\cos(m\theta)\,.$

Consequently, one has

$\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}R_{m}=-\partial_{r}\left(A(t)r^{m}\right)_{|r=R_{0}}=-\sigma
R_{0}^{-3}m(m^{2}-1)R_{m}(t)\,.$

The cubic dispersion relation $\omega=-R^{-3}\sigma m(m^{2}-1)$ shows that the
modes $m=0$ (expansion $m_{0}$ of the circular droplet) and $m=1$
(infinitesimal translation of the circular droplet) are marginally stable,
alluding to mass conservation and translational symmetry. On the other hand,
all $m\geq 2$ modes (morphological deformations) are stabilized by the surface
tension $\sigma$.

Indeed, the solution of the last equation is then

$R_{m}(t)=R_{m}(t=0)e^{\omega t}\,.$

The quantity $\omega$ is often called growth rate, since it determines the
”growth” of the perturbation $R_{m}$. If $\omega>0$, $R_{m}$ grows, and the
perturbation is unstable; otherwise, it is stable. Stability obviously depends
on the wave number of the perturbation, with the perturbations characterized
by large $m$ being stable (due to the term proportional to $-m^{3}$). This is
what we expect due to known stabilizing effect of surface tension.

This result expresses the competition between the destabilizing effect of
viscosity contrast (destabilizing), and the surface tension (stabilizing).

## References

  * [1] T. Alazard, M. Magliocca, N. Meunier, Mathematical study of a coupled incompressible Darcy’s free boundary problem with surface tension, submitted.
  * [2] T. Alazard, N. Meunier and D. Smets, Lyapounov functions, identities and the Cauchy problem for the Hele-Shaw equation, Commun. Math. Phys., volume 377, 1421-1459 (2020)
  * [3] G. Allaire, Conception optimale des structures, Springer.
  * [4] L. Ambrosio, N. Fusco and D. Pallara, Functions of bounded variation and free discontinuity problems Oxford mathematical monographs, 2000.
  * [5] L. Berlyand, J. Fuhrmann, V. Rybalko, Bifurcation of traveling waves in a Keller-Segel type free boundary model of cell motility, Comm. Math. Sciences (2018).
  * [6] L. Berlyand, V. Rybalko, Emergence of traveling waves and their stability in a free boundary model of cell motility. Transactions of the American Mathematical Society, 376, 1799–1844 (2023).
  * [7] C. Blanch-Mercader and J., Casademunt, Spontaneous motility of actin lamellar fragments, Phys. Rev. Lett. (2013), 110, 078102,
  * [8] J.U. Brackbill, D.B. Kothe and C. Zemach, A continuum method for modeling surface tension. Journal of Computational Physics; 100:335–354, 1992.
  * [9] A. Cucchi, A. Mellet, N. Meunier, A Cahn-Hilliard model for cell motility, SIAM J. Math. Anal., 52(4), 3843–3880, 2020.
  * [10] A. Cucchi, A. Mellet, N. Meunier, Self polarization and traveling wave in a model for cell crawling migration. Discrete and dynamical systems, 42 (2022), 2381–2407.
  * [11] M. C. Dallaston, Mathematical models of bubble evolution in a hele-shaw cell. PhD thesis, Queensland University of Technology, 2013.
  * [12] O. Gallinato, M. Ohta, C. Poignard and T. Suzuki, Free boundary problem for cell protrusion formations: theoretical and numerical aspects. J. Math Biol.,75(2):263-307, 2017.
  * [13] L. Giacomelli and F. Otto, Variatonal formulation for the lubrication approximation of the Hele-Shaw flow. Calculus of Variations and Partial Differential Equations, 13(3), 377–403, 2001.
  * [14] C. Hirt and B. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, Journal of Computational Physics, 39(1):201–225, 1981.
  * [15] I. Lavi, N. Meunier, O. Pantz, Implicit time discretization for coupled one-phase Hele-Shaw problem with surface tension, in redaction.
  * [16] I. Lavi, N. Meunier, R. Voituriez, J., Casademunt, Motility and morphodynamics of confined cells, Phys. Rev. E. , 110, 078102 (2020).
  * [17] F. Otto, Dynamics of labyrinthine pattern formation in magnetic fluids: a mean-field theory. Arch. Rational Mech. Anal. 141 (1998)
  * [18] C. Peskin, The immersed boundary method. Acta Numer. 11, 479–517, 2002.
  * [19] F. Ziebert and I. Aronson, Computational approaches to substrate-based cell motility. npj Computational Materials, 2016.
  * [20] S. Walker, The Shapes of Things: A Practical Guide to Differential Geometry and the Shape Derivative 1st, vol. 28, Advances in Design and Control, SIAM, 2015.