# Energy estimate of a Discrete Duality Finite Volume scheme for a phase-field model with surfactants

Margherita Castellano, Ludovic Goudenège, Flore Nabet

► To cite this version:

Margherita Castellano, Ludovic Goudenège, Flore Nabet. Energy estimate of a Discrete Duality Finite Volume scheme for a phase-field model with surfactants. DD29 - 29th International Conference on Domain Decomposition Methods, Jun 2025, Milan, Italy. hal-05473939

HAL Id: hal-05473939

https://hal.science/hal-05473939v1

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Energy estimate of a Discrete Duality Finite
Volume scheme for a phase-field model with
surfactants

Margherita Castellano and Ludovic Goudenège [0000-0002-6449-5888] and Flore Nabet [0000-0001-7828-251X]

**Abstract** In this work we establish a discrete energy estimate for a phase-field model describing a ternary mixture of two immiscible fluids (air and water) and surfactants, firstly introduced by [3]. We mainly rely on the phase-field formulation presented in [5], where the authors include a stabilization term to the model analyzed in [4]. The stabilization term ensures that the free energy remains bounded from below, which is essential for deriving the discrete energy estimate. This will later guarantee the existence of an approximate solution and the convergence of the numerical scheme (which will be covered in a forthcoming paper). We propose a spatial discretization using the *Discrete Duality Finite Volume* (DDFV) method, which is particularly well-suited for the simulation of non-linear diffusive problems. This method can handle very general (possibly non-conforming or anisotropic) meshes, while preserving key structural properties of the continuous problem at the discrete level.

# 1 Definition of the model

The phase-separation of a binary system in presence of surfactants is modeled through two continuous order-parameters $\phi$ and $c$. The first one, $\phi$, is the phase-field variable describing the local density of the two fluids, taking values close to $\phi = 1$ in water and $\phi = -1$ in air. The second one, $c$, represents the local concentration of surfactants, taking values between 0 in absence of surfactants (i.e. the case in which the system is only locally occupied either by water or air) and 1 when the

<footnote>Margherita Castellano, Flore Nabet
CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120, Palaiseau, France
e-mail: {margherita.castellano,flore.nabet} @polytechnique.edu</footnote>

<footnote>Ludovic Goudenège
Laboratoire de Mathématiques et Modélisation d'Évry, UMR CNRS 8071, Université Évry-Paris-Saclay, 23 Bd de France, 91037 Évry, France
e-mail: ludovic.goudenege@math.cnrs.fr</footnote>
local volume is only occupied by surfactants. We introduce the total free energy of the system,

$$ \mathcal{E}: W^{1,4}(\Omega) \times H^1(\Omega) \to \mathbb{R} \\ (\phi, c) \mapsto \int_{\Omega} \left( \frac{\epsilon_{\phi}}{2} |\nabla \phi|^2 + \frac{\epsilon_c}{2} |\nabla c|^2 \right) + \int_{\Omega} \left( \frac{1}{\epsilon_{\phi}} f_{\phi}(\phi) + \frac{1}{\epsilon_c} f_c(c) + F(\phi, c) \right) \quad (1) $$

which incorporates the classical Cahn-Hilliard energy contributions for both the phase-field variable $\phi$ and the surfactant concentration $c$, and a coupling potential $F$ ensuring some required physical properties for the interaction air-water-surfactants. The first integral represents the interfacial energy, encouraging smooth transitions between phases. We denoted by $\epsilon_{\phi}$ and $\epsilon_c$ the interface thickness for $\phi$ (interface between water and air) and for $c$ (interface between regions with and without surfactants), respectively. The second integral contains the classical double-well potentials $f_{\phi}$ and $f_c$, ensuring phase separation, defined by $f_{\phi}(\phi) = \frac{1}{4}(\phi^2 - 1)^2$ and $f_c(c) = c^2(1-c)^2$, respectively, and having minima at the pure phases $\phi = \pm 1$ and $c = \{0, 1\}$. In the third integral, the coupling potential $F$ is given by

$$ F(\phi, c) = -\alpha |\nabla \phi|^2 c + p(\phi, c) + \frac{\delta}{4} |\nabla \phi|^4, \quad \text{with } p(\phi, c) = \beta \phi^2 c - \gamma \phi^3 c, $$

and $\alpha, \beta, \gamma, \delta \in \mathbb{R}_+^*$, some positive constants. Here, the first term $-\alpha |\nabla \phi|^2 c$ is a coupling term between the interface and surfactants, encouraging a high presence of surfactants along the interfaces. The second term $\beta \phi^2 c$ penalizes the presence of surfactants in pure phases, while the term $-\gamma \phi^3 c$ represents the different solubility of surfactants in air and water, penalizing the presence of surfactants in the air. Finally, the last term $\frac{\delta}{4} |\nabla \phi|^4$ is the stabilization term firstly introduced by [5] to ensure the well-posedness of the problem. As in the classical Cahn-Hilliard equation, the system is the $H^{-1}(\Omega) \times H^{-1}(\Omega)$ gradient flow associated with the energy $\mathcal{E}$ defined in (1), which leads to the following coupled Cahn-Hilliard-type system: Find $c: (0, T) \times \Omega \to \mathbb{R}$ and $\phi: (0, T) \times \Omega \to \mathbb{R}$, such that for any $(t, x) \in (0, T) \times \Omega$,

$$ \left\{ \begin{array}{l} \partial_t \phi = \Delta \mu, \\ \mu = -\epsilon_{\phi} \Delta \phi + \frac{1}{\epsilon_{\phi}} f'_{\phi}(\phi) + 2\alpha \operatorname{div}(c \nabla \phi) + \partial_{\phi} p(\phi, c) - \delta \operatorname{div}(|\nabla \phi|^2 \nabla \phi), \\ \partial_t c = \Delta \eta, \\ \eta = -\epsilon_c \Delta c + \frac{1}{\epsilon_c} f'_{c}(c) - \alpha |\nabla \phi|^2 + \partial_c p(\phi, c), \end{array} \right. \quad (2) $$

where $\Omega$ is a bounded open set of $\mathbb{R}^2$ and $\mu$ and $\eta$ are intermediate variables, called chemical potentials, representing the variational derivatives of the total free energy (1) with respect to $\phi$ and $c$, respectively. The system is completed with initial conditions in $\Omega$: $\phi(0, \cdot) = \phi_0 \in H^1(\Omega)$ and $c(0, \cdot) = c_0 \in H^1(\Omega)$, and homogeneous Neumann boundary conditions:
$$ \nabla \phi \cdot \mathbf{n} = \nabla_c \cdot \mathbf{n} = \nabla \mu \cdot \mathbf{n} = \nabla \eta \cdot \mathbf{n} = 0 \quad \text{on } (0, T) \times \partial\Omega, $$

where $\mathbf{n}$ is the outward unit normal vector to the boundary $\partial\Omega$.

In Section 2 we give a general overview of the DDFV method and notations (we refer the reader to [1] for a complete description). In Section 2.3, we provide the discrete scheme and introduce two discrete operators, specifically designed to ensure the discrete energy estimate proved in Section 3.

# 2 DDFV framework

## 2.1 The DDFV meshes and notations

We call a DDFV mesh $\mathcal{T}$ of the domain $\Omega$ the triple $(\mathfrak{M}, \mathfrak{M}^*, \mathfrak{D})$ of meshes defined as follows (see Fig. 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1: A DDFV mesh $\mathcal{T} = (\mathfrak{M}, \mathfrak{M}^*, \mathfrak{D})$.</figcaption></figure>

We partition the domain $\Omega$ in a primal mesh $\mathfrak{M}$ composed by disjoint polygons $K$, called primal control volumes, such that their union covers the entire domain $\tilde{\Omega}$. To each primal cell $K$ we associate a center $x_K$ in its interior. Based on the primal mesh, we define the associated dual mesh $\mathfrak{M}^*$, whose control volumes $K^*$ are built around the vertices of the primal mesh. Each dual cell $K^*$ is defined by joining the centers of the primal cells surrounding the vertex. Such vertex is the center of the dual cell $K^*$, denoted by $x_{K^*}$. For each couple of edges $(\sigma, \sigma^*)$ (where $\sigma = (x_{K^*}, x_{L^*})$ separate the primal control volumes $K$ and $L$, and $\sigma^* = (x_K, x_L)$ separate the dual control volumes $K^*$ and $L^*$), we define the diamond cell $D$ as the quadrilateral having $\sigma$ and $\sigma^*$ as diagonals. The set of all diamonds covers the entire domain $\Omega$ and is denoted by $\mathfrak{D}$. If $\sigma$ is an edge on the boundary $\partial\Omega$, the quadrilateral $D \in \mathfrak{D}^{ext}$ degenerate in a triangle, having as vertices the two endpoints of $\sigma$ and the center $x_K$ of the primal cell $K$ adjacent to $\sigma$. We denote by $\mathfrak{D}_K$ (resp. $\mathfrak{D}_{K^*}$) the set of diamonds having $K$
(resp. $K^*$) as one of their primal (resp. dual) cells. Finally, for any polygon $U$, we denote by $m_U$ the 2-dimensional Lebesgue measure of $U$.

**Nota bene:** The mesh $\mathcal{T}$ is supposed to satisfy some regularity factors, which are important for the convergence analysis of the scheme. We do not detail them here, as we will not talk about convergence, but we refer the reader to [1] for a complete description. However, such regularity factors are not very restrictive, and only limit the flattening of the diamonds in the mesh refining process.

## 2.2 Discrete unknowns and discrete operators

We define two types of discrete unknowns. The concentration fields $\phi$ and $c$ are approximated by scalar unknowns $u_{\mathcal{T}} \in \mathbb{R}^{N_{\mathcal{T}}}$ (where $N_{\mathcal{T}}$ denotes the total number of primal and dual cells) constant over the primal and dual cells. On the diamond mesh, we will approximate the gradients of such concentrations. Therefore, we also introduce vector unknowns $\xi_{\mathcal{D}} \in (\mathbb{R}^2)^{N_{\mathcal{D}}}$ constant over the diamond cells (where $N_{\mathcal{D}}$ denotes the number of diamond cells). Hence,

$$u_{\mathcal{T}} = ((u_K)_{K \in \mathcal{M}}, (u_{K^*})_{K^* \in \mathcal{M}^*}) \in \mathbb{R}^{N_{\mathcal{T}}} \quad \text{and} \quad \xi_{\mathcal{D}} = ((\xi_D)_{D \in \mathcal{D}}) \in (\mathbb{R}^2)^{N_{\mathcal{D}}}.$$

We define the following discrete dot products

$$\forall u_{\mathcal{T}}, v_{\mathcal{T}} \in \mathbb{R}^{N_{\mathcal{T}}}, \quad \llbracket u_{\mathcal{T}}, v_{\mathcal{T}} \rrbracket_{\mathcal{T}} = \frac{1}{2} \left( \sum_{K \in \mathcal{M}} m_K u_K v_K + \sum_{K^* \in \mathcal{M}^*} m_{K^*} u_{K^*} v_{K^*} \right), \\ \forall \xi_{\mathcal{D}}, \psi_{\mathcal{D}} \in (\mathbb{R}^2)^{N_{\mathcal{D}}}, \quad (\xi_{\mathcal{D}}, \psi_{\mathcal{D}})_{\mathcal{D}} = \sum_{D \in \mathcal{D}} m_D \xi_D \cdot \psi_D,$$

with the associated $L^2$-norms

$$\forall u_{\mathcal{T}} \in \mathbb{R}^{N_{\mathcal{T}}}, \quad \|u_{\mathcal{T}}\|_{\mathcal{T}} = \llbracket u_{\mathcal{T}}, u_{\mathcal{T}} \rrbracket_{\mathcal{T}}^{1/2}, \quad \text{and} \quad \forall \xi_{\mathcal{D}} \in (\mathbb{R}^2)^{N_{\mathcal{D}}}, \quad \|\xi_{\mathcal{D}}\|_{\mathcal{D}} = (\xi_{\mathcal{D}}, \xi_{\mathcal{D}})_{\mathcal{D}}^{1/2}.$$

In the DDFV setting the discrete gradient of a scalar field $\nabla^{\mathcal{D}} : \mathbb{R}^{N_{\mathcal{T}}} \to (\mathbb{R}^2)^{N_{\mathcal{D}}}$ and the discrete divergence of a vector field $\text{div}^{\mathcal{T}} : (\mathbb{R}^2)^{N_{\mathcal{D}}} \to \mathbb{R}^{N_{\mathcal{T}}}$ (see [1] for a precise definition) are in duality, mimiquing the continuous setting. This property gives a meaning to the name of the method and is established through the following discrete Green's theorem (where we have omitted the boundary term since we work with homogeneous Neumann boundary conditions).

**Theorem 1 (Green's theorem with Neumann homogeneous boundary conditions)**
For all $(\xi_{\mathcal{D}}, u_{\mathcal{T}}) \in (\mathbb{R}^2)^{N_{\mathcal{D}}} \times \mathbb{R}^{N_{\mathcal{T}}}$ such that $\xi_{\mathcal{D}} \cdot \mathbf{n} = 0$ on $\partial\Omega$, we have

$$\llbracket \text{div}^{\mathcal{T}} \xi_{\mathcal{D}}, u_{\mathcal{T}} \rrbracket_{\mathcal{T}} = -(\xi_{\mathcal{D}}, \nabla^{\mathcal{D}} u_{\mathcal{T}})_{\mathcal{D}}. \quad (3)$$
## 2.3 DDFV discretization of the model

We derive a DDFV approximation scheme for the coupled Cahn–Hilliard system (2) by integrating the equations over each primal and dual control volume, $K \in \mathfrak{M}$ and $K^* \in \mathfrak{M}^*$, respectively. For the time discretization, let $N \in \mathbb{N}^*$ be given, and define the time step $\Delta t = \frac{T}{N}$. The time interval $[0, T]$ is then uniformly partitioned by setting $t^n = n\Delta t$, for $n = \{0, \dots, N\}$. We employ the implicit Euler method for most terms. However, a semi-implicit treatment is applied to the derivatives of the double-well potentials $f_\phi$, $f_c$, as well as to the polynomial coupling potential $p$: for any $T = \{K \in \mathfrak{M}, K^* \in \mathfrak{M}^*\} \in \mathcal{T}$, we define them as

$$ 
\begin{aligned}
d_v^{f_v}(v_T^n, v_T^{n+1}) &= \frac{f_v(v_T^{n+1}) - f_v(v_T^n)}{(v_T^{n+1} - v_T^n)}, \quad \text{with } v = \{\phi, c\} \\
d_\phi^p(\phi_T^n, \phi_T^{n+1}, c_T^{n+1})(\phi_T^{n+1} - \phi_T^n) + d_c^p(\phi_T^{n+1})(c_T^{n+1} - c_T^n) \\
&= p(\phi_T^{n+1}, c_T^{n+1}) - p(\phi_T^n, c_T^n).
\end{aligned} 
\quad (4) $$

Furthermore, a scheme based on the midpoint rule is used for the term involving $\partial_\phi(-\alpha|\nabla\phi|^2c) = 2\alpha \operatorname{div}(c\nabla\phi)$ in the equation for $\mu$, while an explicit discretization is adopted for the term involving $\partial_c(-\alpha|\nabla\phi|^2c) = -\alpha|\nabla\phi|^2$ in the equation for $\eta$. As we shall see in Section 3, this choice is made in order to obtain an unconditionally stable energy estimate. We initialize the problem by taking the average projection of the initial data on the primal and dual meshes: $u_T^0 = \left(\left(\frac{1}{m_T} \int_T u^0(x)dx\right)_{T \in \mathfrak{M} \cup \mathfrak{M}^*}\right)$ for $u_T^0 = \{\phi_T^0, c_T^0\}$. The discrete problem reads as follows. Knowing $(\phi_T^n, c_T^n) \in \mathbb{R}^{N_T} \times \mathbb{R}^{N_T}$, find $(\phi_T^{n+1}, \mu_T^{n+1}, c_T^{n+1}, \eta_T^{n+1}) \in \mathbb{R}^{N_T} \times \mathbb{R}^{N_T} \times \mathbb{R}^{N_T} \times \mathbb{R}^{N_T}$, such that

$$ 
\left\{
\begin{aligned}
& \phi_T^{n+1} - \phi_T^n = \Delta t \operatorname{div}^T (\nabla^\flat \mu_T^{n+1}), \\
& \mu_T^{n+1} = -\epsilon_\phi \operatorname{div}^T (\nabla^\flat \phi_T^{n+1}) + \frac{1}{\epsilon_\phi} d_\phi^{f_\phi} (\phi_T^n, \phi_T^{n+1}) \\
& \qquad + 2\alpha \operatorname{div}^T \left( I_\flat (c_T^{n+1}) \nabla^\flat \left( \frac{\phi_T^n + \phi_T^{n+1}}{2} \right) \right) + d_\phi^p (\phi_T^n, \phi_T^{n+1}, c_T^{n+1}) \\
& \qquad - \delta \operatorname{div}^T (\nabla^\flat \phi_T^{n+1} |\nabla^\flat \phi_T^{n+1}|^2), \\
& c_T^{n+1} - c_T^n = \Delta t \operatorname{div}^T (\nabla^\flat \eta_T^{n+1}), \\
& \eta_T^{n+1} = -\epsilon_c \operatorname{div}^T (\nabla^\flat c_T^{n+1}) + \frac{1}{\epsilon_c} d_c^{f_c} (c_T^n, c_T^{n+1}) - \alpha \mathcal{P}_T (\nabla^\flat \phi_T^n) + d_c^p (\phi_T^n).
\end{aligned}
\right. 
\quad (5) $$

To take into account the homogeneous Neumann boundary condition, we impose for all unknowns $\nabla^D u_T^{n+1} \cdot \mathbf{n}_{\sigma K} = 0$, for any $D \in \mathfrak{D}^{ext}$, where $\mathbf{n}_{\sigma K}$ is the unit normal vector to the boundary edge $\sigma$ outward to the associated primal cell $K$. Note that the concentration $c$ is defined on the primal and dual cells, but its values on the diamond cells are required to discretize $\operatorname{div}(c\nabla\phi)$: the divergence operator acts on the diamond mesh, and $\nabla\phi$ is itself defined on the diamond mesh, so it would
be mathematically inconsistent to multiply quantities defined on different meshes. Similarly, although $\nabla\phi$ is naturally defined on the diamond mesh, the equation for $\eta$ is discretized on the primal and dual cells, requiring its values on these meshes. This motivates the construction of operators $\mathcal{I}_{\mathcal{D}}$ and $\mathcal{P}_{\mathcal{T}}$ so as to be consistent with the mesh structure, and to ensure the discrete energy estimate.

Firstly, the operator $\mathcal{P}_{\mathcal{T}}$ acts as a discrete $L^2$ projection from the space of vector fields $(\mathbb{R}^2)^{N_{\mathcal{D}}}$ onto the space of scalar fields $\mathbb{R}^{N_{\mathcal{T}}}$, we define it as follows

$$ \mathcal{P}_{\mathcal{T}}(\xi_{\mathcal{D}}) = \frac{1}{m_{\mathcal{T}}} \sum_{D \in \mathcal{D}_{\mathcal{T}}} m_{D \cap T} |\xi_{D}|^2, \quad \forall \mathcal{T} = \{K \in \mathfrak{M}, K^* \in \mathfrak{M}^*\} \in \mathcal{T}. \quad (6) $$

Conversely, the operator $\mathcal{I}_{\mathcal{D}}$ maps a scalar field of $\mathbb{R}^{N_{\mathcal{T}}}$ to a scalar field of $\mathbb{R}^{N_{\mathcal{D}}}$, associating with each diamond a weighted average of the values of the primal and dual cells within that diamond. We define the operator $\mathcal{I}_{\mathcal{D}}$ as in [2]:

$$ \mathcal{I}_{\mathcal{D}} : u_{\mathcal{T}} \in \mathbb{R}^{N_{\mathcal{T}}} \mapsto \left( \frac{1}{m_D} \int_D u_{\mathcal{T}}(x) dx \right)_{D \in \mathcal{D}} \in \mathbb{R}^{N_{\mathcal{D}}}. \quad (7) $$

We are now also able to define the natural discrete counterpart of the continuous free energy (1) by replacing the continuous functions and operators with their discrete counterparts. The discrete free energy $\mathcal{E}_{\mathcal{T}}: \mathbb{R}^{N_{\mathcal{T}}} \times \mathbb{R}^{N_{\mathcal{T}}} \to \mathbb{R}$ is then defined as

$$ \begin{aligned} \mathcal{E}_{\mathcal{T}}(\phi_{\mathcal{T}}, c_{\mathcal{T}}) := & \frac{\epsilon_{\phi}}{2} \|\nabla_{\mathcal{D}}\phi_{\mathcal{T}}\|_{\mathcal{D}}^2 + \frac{\epsilon_{c}}{2} \|\nabla_{\mathcal{D}}c_{\mathcal{T}}\|_{\mathcal{D}}^2 + \frac{1}{\epsilon_{\phi}} \left[ \|f_{\phi}(\phi_{\mathcal{T}}), \mathbf{1}_{\mathcal{T}}\|_{\mathcal{T}} \right. \\ & \left. + \frac{1}{\epsilon_{c}} \left[ \|f_{c}(c_{\mathcal{T}}), \mathbf{1}_{\mathcal{T}}\|_{\mathcal{T}} - \alpha \sum_{D \in \mathcal{D}} m_{D} \mathcal{I}_{D}(c_{\mathcal{T}}) |\nabla^{D}\phi_{\mathcal{T}}|^{2} \right] \right. \\ & \left. + \left[ \|p(\phi_{\mathcal{T}}, c_{\mathcal{T}}), \mathbf{1}_{\mathcal{T}}\|_{\mathcal{T}} + \frac{\delta}{4} \sum_{D \in \mathcal{D}} m_{D} |\nabla^{D}\phi_{\mathcal{T}}|^{4} \right] \end{aligned} \quad (8) $$

where $\mathbf{1}_{\mathcal{T}}$ denotes the piece-wise constant function equal to one on each interior primal cell $K \in \mathfrak{M}$, and on each dual cell $K^* \in \mathfrak{M}^*$, and zero elsewhere. Note that $\mathcal{E}_{\mathcal{T}}$ is a consistent approximation of the continuous energy $\mathcal{E}$ defined in (1), in the following sense. For sufficiently regular functions $\phi$ and $c$, if we consider the average projection of such functions on the disks $B_T$ centered in $x_T$ and included in the cell $T = \{K \in \mathfrak{M}, K^* \in \mathfrak{M}^*\}$ i.e. $\mathbb{P}_m^T u = \left(\frac{1}{m_{B_T}} \int_{B_T} u(z) dz\right)_{T \in \mathfrak{M} \cup \mathfrak{M}^*}$. Then the difference $|\mathcal{E}(\phi, c) - \mathcal{E}_{\mathcal{T}}(\mathbb{P}_m^T \phi, \mathbb{P}_m^T c)|$ tends to zero when the mesh size goes to zero. This result follows from [1, Proposition 3.6].

## 3 Discrete energy estimate and numerical results

In this section, we provide a discrete energy estimate for the DDFV scheme (5).
**Proposition 1 (Discrete energy estimation)** Let $\phi_{\tau}^n \in \mathbb{R}^{N_{\tau}}$, $c_{\tau}^n \in \mathbb{R}^{N_{\tau}}$. Suppose there exists a solution $(\phi_{\tau}^{n+1}, \mu_{\tau}^{n+1}, c_{\tau}^{n+1}, \eta_{\tau}^{n+1}) \in \mathbb{R}^{N_{\tau}} \times \mathbb{R}^{N_{\tau}} \times \mathbb{R}^{N_{\tau}} \times \mathbb{R}^{N_{\tau}}$ to problem (5). Then, for all time step $\Delta t$, the following equality holds

$$ 
\begin{aligned}
& \Delta t \|\nabla^{\mathcal{D}} \mu_{\tau}^{n+1}\|_{\mathcal{D}}^2 + \Delta t \|\nabla^{\mathcal{D}} \eta_{\tau}^{n+1}\|_{\mathcal{D}}^2 + \mathcal{E}_{\tau}(\phi_{\tau}^{n+1}, c_{\tau}^{n+1}) - \mathcal{E}_{\tau}(\phi_{\tau}^n, c_{\tau}^n) \\
& \qquad + \frac{\epsilon_{\phi}}{2} \|\nabla^{\mathcal{D}} (\phi_{\tau}^{n+1} - \phi_{\tau}^n)\|_{\mathcal{D}}^2 + \frac{\delta}{4} \sum_{D \in \mathcal{D}} m_D (|\nabla^D \phi_{\tau}^{n+1}|^2 - |\nabla^D \phi_{\tau}^n|^2)^2 \\
& \qquad + \frac{\epsilon_c}{2} \|\nabla^{\mathcal{D}} (c_{\tau}^{n+1} - c_{\tau}^n)\|_{\mathcal{D}}^2 + \frac{\delta}{2} \sum_{D \in \mathcal{D}} m_D |\nabla^D \phi_{\tau}^{n+1}|^2 |\nabla^D \phi_{\tau}^{n+1} - \nabla^D \phi_{\tau}^n|^2 = 0.
\end{aligned}
\quad (9) $$

*Proof.* The idea of the proof is to mimick the classical approach of the continuous setting, in the discrete framework. We compute the dot product between the first equation of (5) and $\mu_{\tau}^{n+1}$, between the second one and $(\phi_{\tau}^{n+1} - \phi_{\tau}^n)$, the third one and $\eta_{\tau}^{n+1}$, the fourth one and $(c_{\tau}^{n+1} - c_{\tau}^n)$. Then, we apply the Green theorem (3) for all terms involving the divergence operator $\text{div}^{\mathcal{T}}$, and then use the relation $b(b-a) = \frac{1}{2}(b^2 - a^2 + (b-a)^2)$ for the terms of involving the discrete gradient of $\phi_{\tau}$ and $c_{\tau}$, with $a = \{\nabla^{\mathcal{D}}\phi_{\tau}^n, \nabla^{\mathcal{D}}c_{\tau}^n\}$ and $b = \{\nabla^{\mathcal{D}}\phi_{\tau}^{n+1}, \nabla^{\mathcal{D}}c_{\tau}^{n+1}\}$ (first and second line in the equation below). Adding everything together gives

$$ 
\begin{aligned}
0 = \Delta t \|\nabla^{\mathcal{D}} \mu_{\tau}^{n+1}\|_{\mathcal{D}}^2 &+ \frac{\epsilon_{\phi}}{2} (\|\nabla^{\mathcal{D}} \phi_{\tau}^{n+1}\|_{\mathcal{D}}^2 - \|\nabla^{\mathcal{D}} \phi_{\tau}^n\|_{\mathcal{D}}^2 + \|\nabla^{\mathcal{D}} (\phi_{\tau}^{n+1} - \phi_{\tau}^n)\|_{\mathcal{D}}^2) \\
&+ \Delta t \|\nabla^{\mathcal{D}} \eta_{\tau}^{n+1}\|_{\mathcal{D}}^2 + \frac{\epsilon_c}{2} (\|\nabla^{\mathcal{D}} c_{\tau}^{n+1}\|_{\mathcal{D}}^2 - \|\nabla^{\mathcal{D}} c_{\tau}^n\|_{\mathcal{D}}^2 + \|\nabla^{\mathcal{D}} (c_{\tau}^{n+1} - c_{\tau}^n)\|_{\mathcal{D}}^2) \\
&+ \frac{1}{\epsilon_{\phi}} [d^{f_{\phi}}(\phi_{\tau}^n, \phi_{\tau}^{n+1}), \phi_{\tau}^{n+1} - \phi_{\tau}^n]_{\mathcal{T}} + \frac{1}{\epsilon_c} [d^{f_c}(c_{\tau}^n, c_{\tau}^{n+1}), c_{\tau}^{n+1} - c_{\tau}^n]_{\mathcal{T}} \\
&+ [d^{p_{\phi}}(\phi_{\tau}^n, \phi_{\tau}^{n+1}, c_{\tau}^n), \phi_{\tau}^{n+1} - \phi_{\tau}^n]_{\mathcal{T}} + [d^{p_c}(\phi_{\tau}^{n+1}), c_{\tau}^{n+1} - c_{\tau}^n]_{\mathcal{T}} \\
&- 2\alpha \left( \mathcal{I}_{\mathcal{D}}(c_{\tau}^{n+1}) \nabla^{\mathcal{D}} \left( \frac{\phi_{\tau}^n + \phi_{\tau}^{n+1}}{2} \right), \nabla^{\mathcal{D}} (\phi_{\tau}^{n+1} - \phi_{\tau}^n) \right)_{\mathcal{D}} \\
&- \alpha [\mathcal{P}_{\mathcal{T}}(\nabla^{\mathcal{D}} \phi_{\tau}^n), c_{\tau}^{n+1} - c_{\tau}^n]_{\mathcal{T}} + \delta (\nabla^{\mathcal{D}} \phi_{\tau}^{n+1} |\nabla^{\mathcal{D}} \phi_{\tau}^{n+1}|^2, \nabla^{\mathcal{D}} (\phi_{\tau}^{n+1} - \phi_{\tau}^n))_{\mathcal{D}}.
\end{aligned}
\quad (10) $$

The discrete operators $\mathcal{P}_{\mathcal{T}}$ and $\mathcal{I}_{\mathcal{D}}$ (see Definitions (6) and (7)) are then defined in such a way that we can obtain

$$ 
\begin{aligned}
& 2(\mathcal{I}_{\mathcal{D}}(c_{\tau}^{n+1})\nabla^{\mathcal{D}}\left(\frac{\phi_{\tau}^n + \phi_{\tau}^{n+1}}{2}\right), \nabla^{\mathcal{D}}(\phi_{\tau}^{n+1} - \phi_{\tau}^n))_{\mathcal{D}} + [\mathcal{P}_{\mathcal{T}}(\nabla^{\mathcal{D}}\phi_{\tau}^n), c_{\tau}^{n+1} - c_{\tau}^n]_{\mathcal{T}} \\
& = \frac{1}{2}\sum_{D \in \mathcal{D}} m_D \mathcal{I}_D(c_{\tau}^{n+1})|\nabla^D \phi_{\tau}^{n+1}|^2 - \frac{1}{2}\sum_{D \in \mathcal{D}} m_D \mathcal{I}_D(c_{\tau}^n)|\nabla^D \phi_{\tau}^n|^2.
\end{aligned}
$$

We also emphasize that the choice of time discretization — using the midpoint rule for the equation on $\mu$, and an explicit scheme for the equation on $\eta$ — allowed us to use the identity $(a+b)(a-b) = a^2 - b^2$ to eliminate the crossed terms
$I_D(c_T^{n+1})|\nabla^D \phi_T^n|^2$ in the above relation (take $a = \nabla^D \phi_T^{n+1}$ and $b = \nabla^D \phi_T^n$). This choice therefore plays a crucial role in obtaining this result and, consequently, in deriving the forthcoming energy estimate. Finally, using the identity $a^3(a-b) = \frac{1}{4}(a^4 - b^4 + (a^2-b^2)^2 + 2a^2(a-b)^2)$, with $a = \nabla^D \phi_T^{n+1}$ and $b = \nabla^D \phi_T^n$, the stabilization term satisfies

$$ (\nabla^D \phi_T^{n+1} |\nabla^D \phi_T^{n+1}|^2, \nabla^D (\phi_T^{n+1} - \phi_T^n))_{\mathcal{D}} = \frac{1}{4} \sum_{D \in \mathcal{D}} m_D (|\nabla^D \phi_T^{n+1}|^4 - |\nabla^D \phi_T^n|^4) \\ + \frac{1}{4} \sum_{D \in \mathcal{D}} m_D (|\nabla^D \phi_T^{n+1}|^2 - |\nabla^D \phi_T^n|^2)^2 + \frac{1}{2} \sum_{D \in \mathcal{D}} m_D |\nabla^D \phi_T^{n+1}|^2 |\nabla^D \phi_T^{n+1} - \nabla^D \phi_T^n|^2. $$

By combining all identities and applying the definition of the discrete energy (8), along with the semi-implicit discretizations (4) of the polynomial terms—which eliminates lines 3 and 4 in equality (10)—we obtain the desired estimate (9). □

To conclude, we performed some numerical simulations to validate the model. We consider the case study of a static air bubble, which remains static in time due to the mass conservation property. We initialize the problem with some random values around $c = 0.5$ for the surfactants concentration, exponentially decreasing around the boundary. Let $\epsilon_\phi = \epsilon_c = 0.044$, the time step $\Delta t = 10^{-4}$ and the mesh size be around $h = 0.0145$. For this numerical test, we have set $\alpha = \beta = 0.03$, $\gamma = 10.8$, $\delta = 1.35 \cdot 10^{-6}$. In Fig. 2, we have overlaid the contour lines of the surfactants concentration on the plot representing the static air bubble. In the figures, we can distinguish in blue the water, in red the air, and in white the interface between air and water. For the surfactants, the black contour lines represent the regions of high surfactants concentration, and in white the regions where there are none.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Surfactants concentration on the dual mesh at time $t_1, t_{50}, t_{1000}$ and $t_{5000}$.</figcaption></figure>

For increasing times, surfactants molecules aggregate to each other, distributing themselves at the interface between air and water.

To support the result stated in Proposition 1, we also present below (Fig. 3) the evolution of the discrete free energy (8), which decreases over time (Fig. 3a), together with the numerically computed values of estimation (9) over the iterations,
corresponding to the machine rounding errors (Fig. 3b). Both plots report the time $t^n = n\Delta t$ for $n \in [0, 500]$.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Validation of Proposition 1</figcaption></figure>

**Acknowledgements** The authors acknowledge the French National Agency of Research for its financial support of the FiVAST project n°ANR-24-CE40-1885

## References

1. Andreianov, B., Boyer, F., Hubert, F.: Discrete Duality Finite Volume schemes for Leray-Lions-type elliptic problems on general 2D meshes. Numerical Methods for Partial Differential Equations **23**(1), 145–195 (2007). DOI 10.1002/num.20170

2. Chainais-Hillairet, C., Krell, S., Mouton, A.: Convergence analysis of a DDFV scheme for a system describing miscible fluid flows in porous media. Numerical Methods for Partial Differential Equations **31**(3), 723–760 (2015). DOI 10.1002/num.21913

3. Laradji, M., Guo, H., Grant, M., Zuckermann, M.J.: The effect of surfactants on the dynamics of phase separation. Journal of Physics: Condensed Matter **4**(32), 6715 (1992). DOI 10.1088/0953-8984/4/32/006

4. van der Sman, R.G.M., Meinders, M.B.J.: Analysis of improved Lattice Boltzmann phase field method for soluble surfactants. Computer Physics Communications **199**, 12–21 (2016). DOI 10.1016/j.cpc.2015.10.002

5. Zhang, J., Chen, C., Wang, J., Yang, X.: Efficient, second order accurate, and unconditionally energy stable numerical scheme for a new hydrodynamics coupled binary phase-field surfactant system. Computer Physics Communications **251**, 107122 (2020). DOI 10.1016/j.cpc.2019.107122