# An energy-stable Smoothed Particle Hydrodynamics discretization of the
Navier-Stokes-Cahn-Hilliard model for incompressible two-phase flows

Xiaoyu Feng<EMAIL_ADDRESS>Zhonghua Qiao
<EMAIL_ADDRESS>Shuyu Sun<EMAIL_ADDRESS>Xiuping Wang
<EMAIL_ADDRESS>
###### Abstract

Varieties of energy-stable numerical methods have been developed for
incompressible two-phase flows based on the Navier-Stokes–Cahn–Hilliard (NSCH)
model in the Eulerian framework, while few investigations have been made in
the Lagrangian framework. Smoothed particle hydrodynamics (SPH) is a popular
mesh-free Lagrangian method for solving complex fluid flows. In this paper, we
present a pioneering study on the energy-stable SPH discretization of the NSCH
model for incompressible two-phase flows. We prove that this SPH method
inherits mass and momentum conservation and the energy dissipation properties
at the fully discrete level. With the projection procedure to decouple the
momentum and continuity equations, the numerical scheme meets the divergence-
free condition. Some numerical experiments are carried out to show the
performance of the proposed energy-stable SPH method for solving the two-phase
NSCH model. The inheritance of mass and momentum conservation and the energy
dissipation properties are verified numerically.

###### keywords:

Smoothed Particle Hydrodynamics , Energy stability , Two-phase Flow , Navier-
Stokes-Cahn-Hilliard model

††journal: Submitted to XXX Journal

The first study of energy-stable ODE discretization and the energy-stable
fully discrete scheme by the SPH method for two-phase problems

The scheme ensures the inheritance of momentum conservation and the energy
dissipation law from the PDE level to the ODE level, and then to the fully
discrete level.

This novel method based on NSCH model enjoys the physical consistency, and the
detailed mathematical proof is also provided.

The scheme helps increase the stability of the numerical method, which allows
much larger time step sizes than the traditional ISPH method.

This energy-stable scheme not only alleviates tensile instability, but it also
captures the behavior of the interface and energy variation well.

[1] 11footnotetext: https://ctpl.kaust.edu.sa
[inst1]organization=Computational Transport Phenomena Laboratory (CTPL),
Division of Physical Sciences and Engineering (PSE), King Abdullah University
of Science and Technology (KAUST),addressline=Thuwal, postcode=23955-6900,
country=Kingdom of Saudi Arabia [inst2]organization=Department of Applied
Mathematics and Research Institute for Smart Energy, The Hong Kong Polytechnic
University,addressline= Hung Hom, Kowloon, Hong Kong, city=Hong Kong

[inst3]organization=Computational Transport Phenomena Laboratory (CTPL),
Division of Computer, Electrical and Mathematical Sciences and Engineering
(CEMSE), King Abdullah University of Science and Technology
(KAUST),addressline=Thuwal, postcode=23955-6900, country=Kingdom of Saudi
Arabia

##  1 Introduction

Initially proposed [1, 2] for astrophysical problems, the smoothed particle
hydrodynamics (SPH) approach has become one of the most popular mesh-free
particle methods and it is now widely used in various applications, such as
multi-phase flow [3], flow in porous media [4], fracture and crack propagation
[5], lava flow [6], magneto-hydrodynamics [7], the evolution of planets and
stars [8], and even movie special effects [9]. Unlike the nodes in other mesh-
free methods [10], the SPH particles can be regarded not only mathematically
as interpolation points but also physically as material components, just like
the “embedded atoms" in molecular dynamics (MD) simulation [11].

SPH possesses many well-known appealing features. First of all, compared to
well-studied grid-based Eulerian numerical methods, SPH naturally takes
advantage of the Lagrangian framework, and it treats the convection term
effectively and stably. Free surfaces, material interfaces, and moving
boundaries can all be traced naturally through this method without tracking or
reconstructing. This mesh-free feature facilitates applications like multi-
phase flow and high-energy events like explosions, high-speed collisions, and
penetrations. Local mass conservation is automatically retained in SPH since
the mass is carried by each particle as a unique property of the particle. By
using the unique “color” property of each particle, it is much easier to
describe complex multi-component or multi-material phenomena, such as the
hydrocarbon components in oil and gas reservoirs. In addition, the SPH method
can handle large deformations and complicated geometry without causing mesh
distortion. With the progress of neighbor searching schemes and parallel
computing strategies [12], SPH enjoys even greater computational performance
and potential. Some review literature for the SPH method can be referred to
[13, 14].

In this paper, we consider incompressible two-phase flows based on the Navier-
Stokes–Cahn–Hilliard (NSCH) model. Even though possessing a lot of the
appealing features mentioned above, SPH still faces a number of open problems
and challenges. First of all, the energy-stability and physical consistency of
SPH for the NS equation have not yet been rigorously studied. Second, the
treatment of incompressibility in SPH remains tricky and challenging; in fact,
most incompressibility treatments in the literature are not energy stable.
Third, to the best of our knowledge, there is not an energy-stable SPH method
for the NSCH system proposed in the literature.

Two common perspectives for the stability analysis of the SPH method include:
(1) tensile instability caused by disordered distributions of particles; (2)
the maximum time step that keeps the simulation stable. For the former, a
novel technique called “particle shifting” has been employed to alleviate
particle gathering [15]. For the second, to date, several attempts have been
made to study the energy variation of the SPH methods, which reveal that the
entropy-increasing (energy dissipation) methods possess better stability [7].
In many approaches, the energy dissipation was explicitly added by an
artificial viscosity term, but the artificial viscosity term can introduce
additional numerical errors (commonly known as numerical diffusion). For the
treatment of the incompressibility of the fluid, two common approaches are
used: (1) approximately simulating incompressible flows with small
compressibility, known as Weakly Compressible SPH (WCSPH); (2) simulating
incompressible flows by enforcing the incompressibility, known as
Incompressible SPH (ISPH)[16]. However, energy stable treatment of
incompressibility has not fully been addressed in the literature. Recently,
Sun and Zhu [17] proposed an energy-stable SPH method for incompressible
single-phase flow that does not use artificial viscosity terms, and it serves
as the basis and inspiration for this paper.

The interface, as a critical element of the two-phase fluid system, can be
modeled either as a sharp interface [18, 19, 20] or as a diffuse interface
[21]. The NSCH model belongs to the diffuse interface approach, and it obeys
thermodynamically consistent energy dissipation laws. The same energy law is
also desired to be retained in the discretized equations. Some techniques,
like convex-concave splitting [22, 23] and the stabilizing approach, are used
to construct energy-stable schemes. Several advances in energy-stable schemes
include the scalar auxiliary variable method (SAV), the invariant energy
quadratization method (IEQ) [24], the exponential time differencing method
(ETD) [25], and the linear energy-factorization method (EF) [26, 27]. The
diffuse interface model based on the EoS is also proposed for modeling complex
two-phase and multi-component mixtures [28, 29, 30]. However, the above
studies are all based on a mesh and the Euler framework. Most of the SPH two-
phase methods are based on the sharp interface model [31]. Researchers rarely
looked into how to combine the SPH method and the diffuse interface model
together [32] and how to design an energy-stable scheme for the system.

The rest of this paper is organized as follows. In Section 2, we review the
NSCH system and give a brief proof of energy law at the PDE level. In Section
3, some basics of the SPH method are given and an energy-stable ODE model is
proposed based on the SPH framework. In Section 4, an energy-stable fully
discrete scheme is well developed. We provide the detailed proof and
derivation of the ODE system and fully discrete scheme. In Section 5, three
numerical examples are presented to validate the scheme. Finally, some
concluding remarks are given in Section 6.

##  2 PDE model and its energy law

We now study a mixture of two immiscible, incompressible fluids in a confined
domain $\Omega\subset\mathbb{R}^{d}(d=2,3)$ and the NSCH model is introduced.
To identify the regions occupied by the two fluids, we introduce a phase
function $\phi$, such that

$\phi(\mathbf{x},t)=\begin{cases}1&\text{ fluid 1 }\\\ -1&\text{ fluid 2
}\end{cases}$

According to the idea of the diffuse interface and gradient flow theory, one
of the most popular thermodynamic theories for inhomogeneous fluids, the total
mixing energy has two contributions: $F_{b}(\phi)$ from the homogeneous part
of the fluid and $F_{\nabla}(\phi)$ from the inhomogeneity of the fluid. The
thermodynamic behavior of the entire two-phase system is governed by the
mixing energy functional $F(\phi,\nabla\phi)$,

$\begin{split}F(\phi,\nabla\phi)&=\int_{\Omega}f(\phi,\nabla\phi)d\mathbf{x}=F_{b}(\phi)+F_{\nabla}(\phi)=\int_{\Omega}f_{b}(\phi)d\mathbf{x}+\int_{\Omega}f_{\nabla}(\phi)d\mathbf{x}\\\
&=\int_{\Omega}\left\\{f_{b}(\phi)+\frac{\lambda}{2}\|\nabla\phi\|^{2}\right\\}d\mathbf{x},\end{split}$
(2.1)

where $f$ represents the energy density, and $\lambda$ denotes the
characteristic strength of the phase mixing energy with respect to $\phi$.
$\lambda$ has a relation with the surface tension coefficient $\sigma$ at the
equilibrium state: $\lambda=\frac{3\sigma}{2\sqrt{2}}\varepsilon$, where
$\varepsilon$ is the capillary width of the interface thickness. The energy
density function from homogeneity, $f_{b}(\phi)$, only depends on $\phi$
locally and
$f_{b}(\phi)=\lambda\left(\phi^{2}-1\right)^{2}/\left(4\varepsilon^{2}\right)$.
The equilibrium process of a two-phase system minimizes the total mixing
energy (2.1). Through the variational derivative of the energy functional $F$
with respect to $\phi$, we obtain the chemical potential $\mu$ and
$\mu=\mu_{b}+\mu_{\nabla}$. By setting the coefficient
$\lambda_{b}=\lambda/\varepsilon^{2}$, we have
$\displaystyle\mu_{b}:=\frac{\delta
F_{b}}{\delta\phi}=f_{b}^{\prime}\left(\phi\right)=\lambda_{b}(\phi^{3}-\phi)$,
and $\displaystyle\mu_{\nabla}:=\frac{\delta
F_{\nabla}}{\delta\phi}=-\lambda\Delta\phi.$

The dynamics of the phase function $\phi$ can be determined by an $H^{-1}$
gradient flow of (2.1), which leads to the following Cahn-Hilliard (CH)
equations:

$\displaystyle\frac{D\phi}{Dt}=\frac{\partial\phi}{\partial
t}+(\mathbf{v}\cdot\nabla)\phi=\nabla\cdot(M\nabla\mu),$ (2.2a)
$\displaystyle\mu=\lambda_{b}(\phi^{3}-\phi)-\lambda\Delta\phi.$ (2.2b)

Here, $M$ is a mobility parameter related to the relaxation time scale, and
$\mathbf{v}$ is the velocity of the fluids. When it turns to the Navier-Stokes
(NS) part, for the two-phase system, the momentum equation takes the usual
form below.

$\rho\frac{D\mathbf{v}}{Dt}=\rho\left(\frac{\partial\mathbf{v}}{\partial
t}+(\mathbf{v}\cdot\nabla)\mathbf{v}\right)=-\nabla p+\nabla\cdot\left(\eta
D(\mathbf{v})\right)-\lambda\nabla\cdot(\nabla\phi\otimes\nabla\phi),$

where $\rho$ denotes the density, $p$ is the pressure, and the shear viscosity
is denoted by $\eta$. $-\lambda\nabla\phi\otimes\nabla\phi$ represents the
extra elastic stress induced by the interfacial tension. With the aid of
modified pressure and the divergence theorem, we arrive at the mathematical
formula for the matched density and viscosity case:

$\rho\frac{D\mathbf{v}}{Dt}=-\nabla p+\eta\Delta\mathbf{v}+\mu\nabla\phi.$
(2.3)

with the continuity equation characterizing the mass conservation:

$\frac{D\rho}{Dt}=\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathbf{v})=0.$ (2.4)

For incompressible flows, (2.4) can be converted into the divergence-free
condition:

$\nabla\cdot\mathbf{v}=0.$ (2.5)

The periodic boundary condition is applied here:

$\Phi(x)=\Phi(x+nL),\quad n\in Z,$ (2.6)

where $\Phi$ denotes a general variable, and it can be $\phi$, $p$,
$\mathbf{v}$, $\rho$ and so on. $L$ is a period of function. This periodic
boundary condition makes
$\oint_{\partial\Omega}\Phi(\mathbf{x})\cdot\mathbf{n}ds=0$ hold.

We now derive the energy law at the PDE level for the NSCH system includes
(2.2a) , (2.2b), (2.3), (2.5) and (2.6). Here and after, for any functions
$f,g\in L^{2}(\Omega)$, we use $(f,g)$ to denote
$\int_{\Omega}fg\mathrm{~{}d}\mathbf{x}$ and $\|f\|^{2}=(f,f)$.

######  Proposition 1.

The NSCH system satisfies the following energy dissipation law:

$\frac{dE_{\text{total
}}}{dt}=\frac{dE_{k}}{dt}+\frac{dF}{dt}=-\eta\|\nabla\mathbf{v}\|^{2}-M\|\nabla\mu\|^{2}\leq
0,$ (2.7)

where the kinetic energy $\displaystyle
E_{k}:=\frac{\rho}{2}\|\mathbf{v}\|^{2}$ and the total energy $E_{\text{total
}}=E_{k}+F$.

Proof By taking the inner product of equation (2.3) with $\mathbf{v}$ and
applying the divergence-free condition (2.5), we get

$\frac{dE_{k}}{dt}=\frac{d(\frac{1}{2}\rho\|\mathbf{v}\|^{2})}{dt}=(\rho\frac{D\mathbf{v}}{Dt},\mathbf{v})=(-\nabla
p,\mathbf{v})+(\eta\Delta\mathbf{v},\mathbf{v})+(\mu\nabla\phi,\mathbf{v}).$

Combining the inner product of equation (2.2a) with $\mu$ and the inner
product of equation (2.2b) with $\frac{\partial\phi}{\partial t}$, we have

$\frac{dF}{dt}=\frac{d}{dt}\left\\{\left(F_{b}(\phi),1\right)+\frac{\lambda}{2}\|\nabla\phi\|^{2}\right\\}=-(\mathbf{v}\cdot\nabla\phi,\mu)-M\|\nabla\mu\|^{2}.$
(2.8)

Based on the divergence theorem and the boundary condition, the inner product
of the interfacial term $(\mu\nabla\phi,\mathbf{v})$ in the NS system and the
inner product of the convection term $(\mathbf{v}\cdot\nabla\phi,\mu)$ in the
CH system is equal. These two terms characterize the energies transfer between
the kinetic energy $E_{k}$ and the free energy $F$. We define
$\text{transfer}_{\mathrm{F-EK}}:=\int_{\Omega}(\mu\nabla\phi,\mathbf{v})d\mathbf{x}$
as the energy quantity transferred from the free energy to the kinetic energy
per unit time. The subscript $\mathrm{F-EK}$ means the energy transfer from
$F$ to $E_{k}$ and vice versa. In the system, we have the relation

$\text{transfer}_{\mathrm{F-EK}}=-\text{transfer}_{\mathrm{EK-F}}.$

For the NS system, we have $\displaystyle(-\nabla
p,\mathbf{v})=\int_{\partial\Omega}p\mathbf{v}\cdot\mathbf{n}ds+\int_{\Omega}p\nabla\cdot\mathbf{v}d\mathbf{x}=0$.
The energy conservation laws for the CH system and the NS system can be
concluded, respectively, as follows,

$\frac{dF}{dt}=\operatorname{transfer}_{\mathrm{EK}-\mathrm{F}}-M\|\nabla\mu\|^{2},\quad\frac{dE_{k}}{dt}=\operatorname{transfer}_{\mathrm{F}-\mathrm{EK}}-\eta\|\nabla\mathbf{v}\|^{2}.$

Consequently, the desired energy law (2.7) for the NSCH system is obtained.

##  3 SPH discretization and the NSCH ODE system

In this section, some basics of the SPH method will be given first. Then a
physical consistent NSCH ODE system will be derived based on the SPH
discretization in space.

There are two main steps to the SPH method: the integral representation and
the particle approximation. At the integral representation step, a physical
field function $f(x)$ is converted into its integral representation based on
the convolution of variables through a smoothed and weighted function, which
gives

$\langle
f(\mathbf{x})\rangle_{h}=\int_{\Omega}f\left(\mathbf{x}^{\prime}\right)W\left(\mathbf{x}-\mathbf{x}^{\prime},h\right)d\mathbf{x}^{\prime}.$
(3.1)

Here, $\Omega$ is the supporting domain, and with a smoothing length $h$, a
smoothing kernel function $W\left(\mathbf{x}-\mathbf{x}^{\prime},h\right)$ is
used to approximate the $\delta$ function in the convolution. The SPH
approximation is represented by the symbol $\langle\rangle_{h}$, or by the
subscript $h$, as in $f_{h}(\mathbf{x})$. The expression specifies the
contribution to any field variable at position $\mathbf{x}$ by a particle at
$\mathbf{x}^{\prime}$ that lies within the compact support of the kernel
function. The kernel function becomes a Dirac delta function when
$h\rightarrow 0$, namely $\lim_{h\rightarrow
0}W\left(\mathbf{x}-\mathbf{x}^{\prime},h\right)=\delta\left(\mathbf{x}-\mathbf{x}^{\prime}\right)$,
and it must satisfy the unity condition:

$\int_{\Omega}W\left(\mathbf{x}-\mathbf{x}^{\prime},h\right)d\mathbf{x}^{\prime}=1.$

Particle approximation: it is achieved by replacing the integral
representation of the field function with the sum of all particle values. It
also means that these material particles interact with each other, and the
kernel function controls the range of their effects. Function values at any
position in the domain can be approximated through

$f_{h}(\mathbf{x})=\sum_{j}f\left(\mathbf{x}_{j}\right)W\left(\mathbf{x}-\mathbf{x}_{j},h\right)V_{j}=\sum_{j}\frac{m_{j}}{\rho_{j}}f\left(\mathbf{x}_{j}\right)W\left(\mathbf{x}-\mathbf{x}_{j},h\right),$
(3.2)

where $V_{j}$ is the volume element of one particle $j$ at $\mathbf{x}_{j}$,
and it has been replaced by the ratio between the mass and density:
$V_{j}=m_{j}/\rho_{j}$. If we consider the location of one particle $i$:

$f_{h}(\mathbf{x}_{i})=\sum_{j}\frac{m_{j}}{\rho_{j}}f\left(\mathbf{x}_{j}\right)W\left(\mathbf{x}_{i}-\mathbf{x}_{j},h\right).$
(3.3)

Here, $i$ and $j$ are particle indexes, and particle $j$ is one of neighboring
particles for $i$. The physical variable $f\left(\mathbf{x}_{i}\right)$ can be
approximated using the interpolation formula with the corresponding physical
variables of neighboring particles $f\left(\mathbf{x}_{j}\right)$. We also use
$\mathbf{x}_{ij}$ and $r_{ij}$ to replace $\mathbf{x}_{i}-\mathbf{x}_{j}$ and
$\|\mathbf{x}_{i}-\mathbf{x}_{j}\|$ respectively in the following content. It
must be noted that the kernel function applied in this paper is the Quintic
Spline SPH kernel.

$W\left(\mathbf{x}_{ij},h\right)=\sigma_{d}\begin{cases}(3-q)^{5}-6(2-q)^{5}+15(1-q)^{5},&0\leq
q<1,\\\ (3-q)^{5}-6(2-q)^{5},&1\leq q<2,\\\ (3-q)^{5},&2\leq q<3,\\\ 0,&3\leq
q.\end{cases}$

If we set the support radius of the domain of influence to be $h$, we have
$\hat{h}=h/3$ and $q=\left\|\mathbf{x}_{ij}\right\|/\hat{h}$. $\sigma_{d}$ is
the dimensional normalization factor, where
$\sigma_{1}=1/120\hat{h},\sigma_{2}=7/478\pi\hat{h}^{2}$ and
$\sigma_{3}=3/359\pi\hat{h}^{3}$ in 1D, 2D and 3D domains respectively.

Table 1: Some common SPH operators for discretizing PDE formulas Divergence | $\nabla\cdot\mathbf{A}(\mathbf{x}_{i})$ | $\nabla_{h}\cdot\mathbf{A}(\mathbf{x}_{i})=-\sum_{j}\frac{m_{j}}{\rho_{i}}\left(\mathbf{A}_{i}-\mathbf{A}_{j}\right)\cdot\nabla_{i}W_{ij}$
---|---|---
Gradient | $\nabla f(\mathbf{x}_{i})$ | $\nabla_{h}f\left(\mathbf{x}_{i}\right)=-\sum_{j}\frac{m_{j}}{\rho_{i}}\left(f_{i}-f_{j}\right)\nabla_{i}W_{ij}$
Gradient (symmetric) | $\nabla\frac{f}{\rho}|_{\mathbf{x}=\mathbf{x}_{i}}$ | $\nabla_{h}\frac{f}{\rho}|_{\mathbf{x}=\mathbf{x}_{i}}=\sum_{j}m_{j}\left(\frac{f_{i}}{\rho_{i}^{2}}+\frac{f_{j}}{\rho_{j}^{2}}\right)\nabla_{i}W_{ij}$
Laplace (scalar/vector) | $\nabla^{2}f(\mathbf{x}_{i})$ | $\nabla_{h}^{2}f(\mathbf{x}_{i})=2\sum_{j}\frac{m_{j}}{\rho_{j}}\frac{f(\mathbf{x}_{i})-f(\mathbf{x}_{j})}{\left\|\mathbf{x}_{ij}\right\|^{2}}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right)$
Laplace (vector) | $\nabla^{2}\mathbf{A}(\mathbf{x}_{i})$ | $\nabla_{h}^{2}\mathbf{A}(\mathbf{x}_{i})=2(d+2)\sum_{j}\frac{m_{j}}{\rho_{j}}\frac{\left(\mathbf{A}_{i}-\mathbf{A}_{j}\right)\cdot\mathbf{x}_{ij}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\nabla_{i}W_{ij}$

It is essential to discretize the PDE equations through appropriate SPH
operators to maintain the inheritance of physical properties like mass
conservation, momentum conservation and energy dissipation law at the ODE and
fully discrete level. Some common SPH operators are introduced in Table 1.
$W\left(\mathbf{x}_{ij},h\right)$ is denoted by $W_{ij}$ as well and
$\nabla_{i}W_{ij}$ is used to represent the gradient of the kernel function
$W$ at the position of particle $i$. Therefore, the anti-symmetric property of
gradient can be obtained:

$\nabla_{j}W_{ij}=-\nabla_{i}W_{ij}.$ (3.4)

The number of dimensions is denoted by $d$ in the SPH Laplace operator. In the
following, we will explain why the symmetric or anti-symmetric properties of
the SPH operators or the SPH mathematical expression with physical consistency
are vital. In order to avoid singularities, a minor term of $0.01h^{2}$ is
added to the denominator of the SPH Laplace operator, where
$\left\|\mathbf{x}_{ij}\right\|^{2}$ may vanish when $i$ is equal to $j$ as
follows:

$\nabla_{h}^{2}f(\mathbf{x}_{i})=2\sum_{j}\frac{m_{j}}{\rho_{j}}\frac{f(\mathbf{x}_{i})-f(\mathbf{x}_{j})}{\left\|\mathbf{x}_{ij}\right\|^{2}+0.01h^{2}}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right).$

Based on the anti-symmetric property of gradient as (3.4), the important
negative property of the kernel function can be deduced:

$\nabla_{i}W_{ij}=\mathbf{x}_{ij}\omega\left(r_{ij}\right),$ (3.5)

where $\omega$ is a scalar function and $\omega\leq 0$ always. This negative
property can be easily seen from the Figure 1 as well.

Figure 1: (a) Contour graph for Quintic kernel in 2D; (b) First derivative of
kernel function; (c) Plot of Quintic kernel function and its first derivative
in 1D

From (3.4) and (3.5), we can derive the following lemma. We omit the proof and
leave it to the interested readers.

######  Lemma 3.1.

Assume the SPH operator $\sigma_{ij}$ is symmetric
$(\sigma_{ij}=\sigma_{ji})$, we obtain

$\sum_{i}\sum_{j}m_{i}m_{j}\sigma_{ij}\nabla_{i}W_{ij}=0,$ (3.6)

and if $\sigma_{ij}\geq 0$, we have

$\left\\{\begin{array}[]{l}\sum_{i}\sum_{j}m_{i}m_{j}\sigma_{ij}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right)\leq
0,\\\
\sum_{i}\sum_{j}m_{i}m_{j}\sigma_{ij}\left(\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}\right)\mathbf{v}_{i}\cdot\nabla_{i}W_{ij}\leq
0.\end{array}\right.$ (3.7)

###  3.1 The NSCH ODE system

Firstly, we review the NSCH PDE system in Lagrangian form. Using the identity
$\nabla\cdot(\nabla\phi\otimes\nabla\phi)=\Delta\phi\nabla\phi+\frac{1}{2}\nabla\|\nabla\phi\|^{2}$,
the PDE system can be rewritten as

$\displaystyle\frac{D\phi}{Dt}=\nabla\cdot(M\nabla\mu),$ (3.8a)
$\displaystyle\mu=\lambda_{b}\left(\phi^{3}-\phi\right)-\lambda\Delta\phi,$
(3.8b)
$\displaystyle\frac{D\mathbf{v}}{Dt}=-\nabla\frac{p}{\rho}+\frac{\eta}{\rho}\Delta\mathbf{v}+\frac{\mu_{\nabla}\nabla\phi}{\rho},$
(3.8c) $\displaystyle\nabla\cdot\mathbf{v}=0.$ (3.8d)

The design of ODE expression starts from the definition of free energy at the
ODE level. The homogenous and non-homogenous parts of total free energy are

$\left\\{\begin{array}[]{l}F_{h}^{b}=F_{h}^{b}\left(\phi_{1},\phi_{2},\ldots,\phi_{n}\right),\\\
F_{h}^{\nabla}=F_{h}^{\nabla}\left(\phi_{1},\phi_{2},\ldots,\phi_{n},\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n}\right),\\\
F_{h}=F_{h}^{b}+F_{h}^{\nabla}.\end{array}\right.$

To keep the local and non-local properties respectively, the energy at the ODE
level can be defined through the summation of particles one by one as follows:

$\left\\{\begin{array}[]{l}F_{h}^{b}:=\sum_{i}\frac{m_{i}}{\rho_{i}}f_{b}\left(\phi_{i}\right),\\\
F_{h}^{\nabla}:=\sum_{i}\frac{m_{i}}{\rho_{i}}\frac{\lambda}{2}\left(\nabla_{h}\phi_{i}\right)\cdot\left(\nabla_{h}\phi_{i}\right).\end{array}\right.$
(3.9)

Then, the free energy change with time in the ODE form is given as

$\frac{dF_{h}}{dt}=\frac{dF_{h}^{b}}{dt}+\frac{dF_{h}^{\nabla}}{dt}=\sum_{i}\frac{m_{i}}{\rho_{i}}\frac{\partial
f_{b}\left(\phi_{i}\right)}{\partial\phi_{i}}\frac{d\phi_{i}}{dt}+\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\phi_{i}}\frac{d\phi_{i}}{dt}+\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}\frac{\partial\mathbf{x}_{i}}{\partial
t}.$ (3.10)

The term $\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}\frac{\partial\mathbf{x}_{i}}{\partial
t}$ is responsible for energy transferred from kinetic energy $E_{k}$ to free
energy $F$ at the ODE level, which gives

$\langle\text{transfer}_{\mathrm{EK-F}}\rangle_{h}=\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}\frac{\partial\mathbf{x}_{i}}{\partial
t}.$ (3.11)

The energy-transfer term is originated from the inhomogeneity (the interface
region) and other two terms control the energy dissipation coming from the
diffusive effect of the CH equation. Since the chemical potential
$\mu_{i}^{b}$ locally depends on $\phi_{i}$ as $\mu_{i}^{b}:=\frac{\partial
f_{b}\left(\phi_{i}\right)}{\partial\phi_{i}}$ and $\mu_{i}^{\nabla}$ is a
non-local function given by
$\mu_{i}^{\nabla}:=\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}}{\partial\phi_{i}}$, $\mu$ is defined at the ODE level as
follows:

$\mu_{i}=\frac{\partial
f_{b}\left(\phi_{i}\right)}{\partial\phi_{i}}+\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}}{\partial\phi_{i}}.$ (3.12)

When we turn to the momentum equation of the NS system, in order to keep the
physical consistency of energy transfer between the CH system and the NS
system, the interfacial force should be defined as

$\boldsymbol{\mathcal{T}}_{i}:=\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}.$ (3.13)

This NSCH ODE system guarantees these physical properties and it is finalized
as

$\displaystyle\begin{aligned}
\frac{d\phi_{i}}{dt}=2M\sum_{j}\frac{m_{j}}{\rho_{j}}\frac{\mu_{i}-\mu_{j}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right),\end{aligned}$
(3.14a) $\displaystyle\begin{aligned}
\mu_{i}=\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}}{\partial\phi_{i}}+f_{b}^{\prime}\left(\phi_{i}\right),\end{aligned}$
(3.14b) $\displaystyle\begin{aligned}
\rho_{i}=\sum_{j}m_{j}W_{ij},\end{aligned}$ (3.14c)
$\displaystyle\begin{aligned}
\frac{d\mathbf{v}_{i}}{dt}=&-\frac{\boldsymbol{\mathcal{T}}_{i}}{m_{i}}-\sum_{j}m_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}\right)\nabla_{i}W_{ij}+(d+2)\sum_{j}\frac{m_{j}}{\rho_{i}\rho_{j}}\frac{\left(\eta_{i}+\eta_{j}\right)\left(\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}\right)}{\left\|\mathbf{x}_{ij}\right\|^{2}}\nabla_{i}W_{ij},\end{aligned}$
(3.14d) $\displaystyle\begin{aligned}
-\frac{1}{\rho_{i}}\sum_{j}m_{j}\left(\mathbf{v}_{i}-\mathbf{v}_{j}\right)\cdot\nabla_{i}W_{ij}=0.\end{aligned}$
(3.14e)

Here and after, the notation $\mathbf{v}_{ij}$ is equivalent to
$\mathbf{v}_{i}-\mathbf{v}_{j}$. The appearance of the negative sign of the
interfacial term $-\frac{\boldsymbol{\mathcal{T}}_{i}}{m_{i}}$ in the ODE
system is consistent with the variational derivative of PDE in $H^{-1}$ space.
The viscosity is chosen an average to ensure the symmetric form.

###  3.2 Energy dissipation law at the ODE level

The energy dissipation property of the given ODE system will be proved.
Similar to the verification of the PDE model’s energy dissipation law, the
inner product is utilized here.

######  Theorem 3.1.

Assume that the NSCH ODE system based on SPH discretization as (3.14a),
(3.14b), (3.14c), (3.14d) and (3.14e) is employed, the energy dissipation law
at the ODE level can be maintained as follows:

$\displaystyle\left\langle\frac{dE_{\text{total }}}{dt}\right\rangle_{h}$
$\displaystyle=\left\langle\frac{dE_{k}}{dt}\right\rangle_{h}+\left\langle\frac{dF}{dt}\right\rangle_{h}\leq
0.$ (3.15)

Proof Conservation of energy transfer: We have defined the interfacial force
as (3.13) and the chemical potential that comes from inhomogeneity as (3.12).
Then by taking the inner product of the time derivative term with
$m_{i}\mathbf{v}_{i}$, we get the change of kinetic energy $E_{k}$ per time in
the system:

$\left\langle\frac{dE_{k}}{dt}\right\rangle_{h}=\sum_{i}\left(m_{i}\mathbf{v}_{i},\frac{d\mathbf{v}_{i}}{dt}\right)=\sum_{i}\frac{d\left(\frac{1}{2}m_{i}\left\|\mathbf{v}_{i}\right\|^{2}\right)}{dt}.$
(3.16)

The interfacial term $-\frac{\boldsymbol{\mathcal{T}}_{i}}{m_{i}}$ in the
momentum equation governs the energy transferred from free energy $F$ and the
kinetic energy $E_{k}$. By taking the inner product of this term with
$m_{i}\mathbf{v}_{i}$ and adding up the total for all particles, we obtain

$\langle\text{transfer}_{\mathrm{F-EK}}\rangle_{h}=\sum_{i}\left(-\frac{\boldsymbol{\mathcal{T}}_{i}}{m_{i}},m_{i}\mathbf{v}_{i}\right)=-\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}\frac{\partial\mathbf{x}_{i}}{\partial
t}.$ (3.17)

Thus, we can conclude that

$\langle\text{transfer}_{\mathrm{EK-F}}\rangle_{h}=-\langle\text{transfer}_{\mathrm{F-EK}}\rangle_{h}=\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}\frac{\partial\mathbf{x}_{i}}{\partial
t}.$ (3.18)

The conservation of energy transfer between kinetic energy $E_{k}$ and free
energy $F$ is examined.

Conservation of energy in the CH system: If we take the inner product of
(3.14a) with $\mu_{i}$ and integrate it throughout the whole domain, taking
into account the particle’s control volume $\frac{m_{i}}{\rho_{i}}$, the
results turn out easily by swapping the index $i$ and $j$ and taking advantage
of the SPH negative property as (3.5).

$\begin{split}\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{d\phi_{i}}{dt},\mu_{i}\right)&=\sum_{i}\frac{m_{i}}{\rho_{i}}\left(2M\sum_{j}\frac{m_{j}}{\rho_{j}}\frac{\mu_{i}-\mu_{j}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right),\mu_{i}\right)\\\
&=M\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\frac{\left(\mu_{i}-\mu_{j}\right)^{2}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right)\\\
&=M\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\left(\mu_{i}-\mu_{j}\right)^{2}\omega\left(r_{ij}\right)\leq
0.\end{split}$ (3.19)

Then, by taking the inner product of (3.14b) with $\frac{d\phi_{i}}{dt}$ and
summing for all particles and based on the equation (3.10), we find

$\begin{split}\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\mu_{i},\frac{d\phi_{i}}{dt}\right)&=\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}}{\partial\phi_{i}},\frac{d\phi_{i}}{dt}\right)+\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{\partial
f_{b}\left(\phi_{i}\right)}{\partial\phi_{i}},\frac{d\phi_{i}}{dt}\right)\\\
&=\left\langle\frac{dF_{b}}{dt}\right\rangle_{h}+\left\langle\frac{dF_{\nabla}}{dt}\right\rangle_{h}=\left\langle\frac{dF}{dt}\right\rangle_{h}-\langle\text{transfer}_{\mathrm{F-EK}}\rangle_{h}.\end{split}$
(3.20)

Combining equations (3.19) and (3.20), it reaches the conclusion:

$\left\langle\frac{dF}{dt}\right\rangle_{h}=\langle\text{transfer}_{\mathrm{EK-F}}\rangle_{h}+M\sum_{i}\sum_{i}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\left(\mu_{i}-\mu_{j}\right)^{2}\omega\left(r_{ij}\right).$
(3.21)

Conservation of energy in the NS system: We take the inner product of the
momentum equation (3.14d) with $m_{i}\mathbf{v}_{i}$ and do summation for all
particles. The change of kinetic energy $E_{k}$ at the ODE level is given in
(3.16). The divergence free condition (3.14e) yields

$-\sum_{i}m_{i}\frac{p_{i}}{\rho_{i}^{2}}\sum_{j}m_{j}\left(\mathbf{v}_{i}-\mathbf{v}_{j}\right)\cdot\nabla_{i}W_{ij}=-\sum_{i}\sum_{j}m_{i}m_{j}\frac{p_{i}}{\rho_{i}^{2}}\left(\mathbf{v}_{i}-\mathbf{v}_{j}\right)\cdot\nabla_{i}W_{ij}=0.$

Let $I_{1},I_{2}$ denote the pressure term
$\sum_{i}\left(m_{i}\mathbf{v}_{i},-\sum_{j}m_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}\right)\nabla_{i}W_{ij}\right)$
and the viscosity term
$\sum_{i}\left(m_{i}\mathbf{v}_{i},\frac{\eta_{i}+\eta_{j}}{2\rho_{i}}\nabla^{2}_{h}\mathbf{v}_{i}\right)$,
respectively. Then we have

$I_{1}=-\frac{1}{2}\sum_{i}\sum_{j}m_{i}m_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}\right)\left(\mathbf{v}_{i}-\mathbf{v}_{j}\right)\cdot\nabla_{i}W_{ij}=0,$
(3.22)

and

$\begin{split}I_{2}&=\frac{(d+2)}{2}\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\frac{\left(\eta_{i}+\eta_{j}\right)\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\left(\mathbf{v}_{i}-\mathbf{v}_{j}\right)\cdot\nabla_{i}W_{ij}\\\
&=\frac{(d+2)}{2}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\frac{\left(\eta_{i}+\eta_{j}\right)}{\left\|\mathbf{x}_{ij}\right\|^{2}}\|\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}\|^{2}\omega\left(r_{ij}\right)\leq
0.\end{split}$ (3.23)

Combining (3.16), (3.17), (3.22) and (3.23), we obtain

$\left\langle\frac{dE_{k}}{dt}\right\rangle_{h}=\langle\text{transfer}_{\mathrm{F-EK}}\rangle_{h}+\frac{(d+2)}{2}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\frac{\left(\eta_{i}+\eta_{j}\right)}{\left\|\mathbf{x}_{ij}\right\|^{2}}\|\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}\|^{2}\omega\left(r_{ij}\right).$
(3.24)

Finally, combining (3.18), (3.21) and (3.24), we complete the proof.

###  3.3 Momentum conservation of the ODE system

Next, we turn to the proof of the momentum conservation holds as well for the
designed ODE system.

######  Theorem 3.2.

The NSCH ODE system which consists of equations (3.14a), (3.14b), (3.14c),
(3.14d) and (3.14e) can maintain momentum conservation, namely

$\sum_{i}m_{i}\frac{d\mathbf{v}_{i}}{dt}=0.$ (3.25)

Proof In the first place, by taking the inner product of momentum equation
(3.14d) with $m_{i}$ and doing summation for all particles, we have

$\begin{split}\sum_{i}m_{i}\frac{d\mathbf{v}_{i}}{dt}=&-\sum_{i}\sum_{j}m_{i}m_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{i}^{2}}\right)\nabla_{i}W_{ij}\\\
&+(d+2)\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\frac{\left(\eta_{i}+\eta_{j}\right)\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\nabla_{i}W_{ij}-\sum_{i}\boldsymbol{\mathcal{T}}_{i}.\end{split}$
(3.26)

Two above SPH operators can be treated as a whole symmetric one
($\sigma_{ij}=\sigma_{ji}$), like
$\sigma_{ij}:=-\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}\right)+\frac{(d+2)}{\rho_{i}\rho_{j}}\frac{\left(\eta_{i}+\eta_{j}\right)\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}}{\left\|x_{ij}\right\|^{2}}$.
Then, based on the property (3.6), we find

$\sum_{i}\sum_{j}m_{i}m_{j}\sigma_{ij}\nabla_{i}W_{ij}=0.$ (3.27)

By proving $\sum_{i}\mathcal{T}_{i}=\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}=0$, we can conclude momentum
conservation. It must be noted that the selected SPH gradient operator for
$\phi$ is

$\nabla_{h}\phi_{i}=-\sum_{j}\frac{m_{j}}{\rho_{i}}\left(\phi_{i}-\phi_{j}\right)\nabla_{i}W_{ij}.$

Based on the definition of $F_{h}^{\nabla}$ as (3.9), it becomes

$\begin{split}F^{\nabla}_{h}&=\sum_{k}\frac{m_{k}}{\rho_{k}}\frac{\lambda}{2}\left(\nabla_{h}\phi_{k}\right)\cdot\left(\nabla_{h}\phi_{k}\right)\\\
&=\sum_{k}\frac{m_{k}}{\rho_{k}}\frac{\lambda}{2}\left[-\sum_{l}\frac{m_{l}}{\rho_{l}}\left(\phi_{k}-\phi_{l}\right)\nabla_{k}W_{kl}\right]\cdot\left[-\sum_{q}\frac{m_{q}}{\rho_{q}}\left(\phi_{k}-\phi_{q}\right)\nabla_{k}W_{kq}\right].\end{split}$
(3.28)

We consider the compressible flow here. Hence, we can set: $m_{i}\equiv m$ and
$\rho_{i}\equiv\rho$ and define
$f_{klq}\equiv\left(\phi_{k}-\phi_{l}\right)\left(\phi_{k}-\phi_{q}\right)\nabla_{k}W_{kl}\cdot\nabla_{k}W_{kq}$,
then

$F_{h}^{\nabla}=\frac{\lambda
m^{3}}{2\rho^{3}}\sum_{k}\sum_{l}\sum_{q}f_{klq}.$ (3.29)

Correspondingly, we get

$\sum_{i}\boldsymbol{\mathcal{T}}_{i}=\sum_{i}\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}=\frac{\lambda
m^{3}}{2\rho^{3}}\sum_{k}\sum_{l}\sum_{q}\sum_{i}\frac{\partial
f_{klq}}{\partial\mathbf{x}_{i}}.$ (3.30)

Using the negative property of the kernel function, $f_{klq}$ can be expanded:

$\displaystyle\frac{\partial f_{klq}}{\partial\mathbf{x}_{i}}$
$\displaystyle=\left(\phi_{k}-\phi_{l}\right)\left(\phi_{k}-\phi_{q}\right)\left\\{\frac{\partial\left(\mathbf{x}_{kl}\omega\left(r_{kl}\right)\right)}{\partial\mathbf{x}_{i}}\mathbf{x}_{kq}\omega\left(r_{kq}\right)+\frac{\partial\left(\mathbf{x}_{kq}\omega\left(r_{kq}\right)\right)}{\partial\mathbf{x}_{i}}\mathbf{x}_{kl}\omega\left(r_{kl}\right)\right\\}$
$\displaystyle=\left(\phi_{k}-\phi_{l}\right)\left(\phi_{k}-\phi_{q}\right)\left(\delta_{ki}-\delta_{li}\right)\left\\{\omega\left(r_{kl}\right)+r_{kl}\omega^{\prime}\left(r_{kl}\right)\right\\}\omega\left(r_{kq}\right)\mathbf{x}_{kq}$
$\displaystyle\quad...\quad+\left(\phi_{k}-\phi_{l}\right)\left(\phi_{k}-\phi_{q}\right)\left(\delta_{ki}-\delta_{qi}\right)\left\\{\omega\left(r_{kq}\right)+r_{kq}\omega^{\prime}\left(r_{kq}\right)\right\\}\omega\left(r_{kl}\right)\mathbf{x}_{kl}.$

With the following equations of summation for the $i$ index:

$\left\\{\begin{array}[]{l}\sum_{i}\delta_{ki}=1,\\\
\sum_{i}\left(\delta_{ki}-\delta_{li}\right)=0\quad\text{and}\quad\sum_{i}\left(\delta_{ki}-\delta_{qi}\right)=0.\end{array}\right.$

We reach $\sum_{i}\frac{\partial f_{klq}}{\partial x_{i,\alpha}}=0.$ Finally,
we get $\sum_{i}\boldsymbol{\mathcal{T}}_{i}=0$ and (3.25) holds naturally.

######  Remark 3.1.

Through appropriate selections of SPH operators and the derivation based on
the original energy and physical consistency, the proposed NSCH ODE system
inherits the physical consistency, including mass conservation, momentum
conservation, and the energy dissipation law. The spatial discretization has
been completed based on the SPH framework.

##  4 Fully discrete scheme

Based on the energy-stable ODE system we proposed in the last section, we now
focus on the temporal scheme. In this section, we decouple the CH system and
the NS system and solve them separately. The projection method is adopted for
the solution of the NS system. Some techniques like the convex-concave
splitting approach, stabilizing term, and energy-factorization method are also
utilized to construct the energy-stable fully discrete scheme.

###  4.1 Temporal discretization

Now the fully discrete scheme of NSCH system reads as follows:

$\displaystyle\begin{aligned} \frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}=2M\sum_{j}\frac{m_{j}}{\rho_{j}}\frac{\mu_{i}^{k+1}-\mu_{j}^{k+1}}{\left\|\mathbf{x}_{ij}\right\|^{2}}\left(\mathbf{x}_{ij}\cdot\nabla_{i}W_{ij}\right),\end{aligned}$
(4.1a) $\displaystyle\begin{aligned}
\mu_{i}^{k+1}=\mu_{i}^{b}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)+\mu_{i}^{\nabla}\left(\phi^{k+1}\right)=f_{b}^{\prime}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)\end{aligned}+\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}\left(\phi^{k+1}\right)}{\partial\phi_{i}^{k+1}},$ (4.1b)
$\displaystyle\begin{aligned}
\rho_{i}^{k+1}=\sum_{j}m_{j}W_{ij},\end{aligned}$ (4.1c)
$\displaystyle\begin{aligned}
\frac{\tilde{\mathbf{v}}_{i}^{k+1}-\mathbf{v}_{i}^{k}}{\Delta
t}=-\frac{\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\tilde{\mathbf{x}}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}}{m_{i}}+(d+2)\sum_{j}\frac{m_{j}}{\rho_{i}^{k+1}\rho_{j}^{k+1}}\frac{\left(\eta_{i}+\eta_{j}\right)\left(\tilde{\mathbf{v}}_{ij}^{k+1}\cdot\mathbf{x}_{ij}\right)}{\left\|\mathbf{x}_{ij}\right\|^{2}}\nabla_{i}W_{ij},\end{aligned}$
(4.1d) $\displaystyle\begin{aligned}
\frac{\mathbf{v}_{i}^{k+1}-\tilde{\mathbf{v}}_{i}^{k+1}}{\Delta
t}=-\sum_{j}m_{j}\left(\frac{p_{i}^{k+1}}{\left(\rho_{i}^{k+1}\right)^{2}}+\frac{p_{j}^{k+1}}{\left(\rho_{j}^{k+1}\right)^{2}}\right)\nabla_{i}W_{ij},\end{aligned}$
(4.1e) $\displaystyle\begin{aligned}
\nabla_{h}\cdot\mathbf{v}_{i}^{k+1}=-\frac{1}{\rho_{i}^{k+1}}\sum_{j}m_{j}\left(\mathbf{v}_{i}^{k+1}-\mathbf{v}_{j}^{k+1}\right)\cdot\nabla_{i}W_{ij}=0,\end{aligned}$
(4.1f) $\displaystyle\begin{aligned}
\mathbf{x}_{i}^{k+1}-\mathbf{x}_{i}^{k}=\mathbf{v}_{i}^{k+1}\Delta
t,\end{aligned}$ (4.1g)

where $\phi$ represents the collection of phase function $\phi_{i}$ on each
particle and $\mathbf{x}$ represents the collection of position information
$\mathbf{x}_{i}$ on each particle.

Some terms in the fully discrete system require special treatment to ensure
energy stability. The details of them are presented here. A semi-implicit
scheme applied to term $\mu_{i}^{b}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)$
reads:

$\mu_{i}^{b}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)=f_{b}^{\prime}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)=\lambda_{b}\left[\left(\xi+\left(\phi_{i}^{k}\right)^{2}\right)\phi_{i}^{k+1}-(\xi+1)\phi_{i}^{k}\right].$
(4.2)

where $\xi$ is the energy-factorization parameter in the EF method. It results
in a linear scheme for determining $\phi$. Based on the expression of
$F^{\nabla}_{h}$ as (3.28), a fully implicit scheme applied to term
$\mu_{i}^{\nabla}(\phi^{k+1})$ reads:

$\mu_{i}^{\nabla}\left(\phi^{k+1}\right)=\frac{\rho_{i}}{m_{i}}\sum_{j}\frac{m_{j}}{\rho_{j}}\lambda\left[\sum_{p}\frac{m_{p}}{p_{j}}\left(\delta_{ij}-\delta_{ip}\right)\nabla_{j}W_{jp}\right]\left[\sum_{l}\frac{m_{l}}{\rho_{j}}\left(\phi_{j}^{k+1}-\phi_{l}^{k+1}\right)\nabla_{j}W_{jl}\right].$
(4.3)

Then, we utilize the stabilizing method for the interfacial term
$\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\tilde{\mathbf{x}}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}$.
Based on the convex-concave splitting approach, we have

$\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\tilde{\mathbf{x}}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}=\frac{2\mathcal{M}}{\Delta
t}\tilde{\mathbf{x}}_{i}^{k+1}+\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k}\right)-\frac{2\mathcal{M}}{\Delta
t}\mathbf{x}_{i}^{k}=2\mathcal{M}\tilde{\mathbf{v}}_{i}^{k+1}+\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k}\right),$
(4.4)

where $\mathcal{M}$ is a stabilization parameter.

###  4.2 Energy dissipation law at the fully discrete level

Now we set out to prove the energy dissipation law of this full discrete
system. Since the entire fully discrete scheme is decoupled into two steps, we
update the $\phi$ by solving the CH system in the first step. The free energy
of the system is changed as

$\frac{F_{h}\left(\phi^{k+1}\right)-F_{h}\left(\phi^{k}\right)}{\Delta
t}=\frac{F_{h}\left(\phi^{k+1},\mathbf{x}^{k}\right)-F_{h}\left(\phi^{k},\mathbf{x}^{k}\right)}{\Delta
t}.$

Then, at the second step, particles are moved with the velocity we determined
form the NS system. Both the kinetic energy and the free energy are changed.

$\frac{E_{k,h}\left(\mathbf{v}^{k+1}\right)-E_{k,h}\left(\mathbf{v}^{k}\right)}{\Delta
t}+\frac{F_{h}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)}{\Delta
t}.$

We give two lemmas to support the proof of the energy dissipation law at the
fully discrete level for the CH system and the NS system, respectively.

######  Lemma 4.1.

Assume that the numerical schemes for the CH system are (4.1a), (4.1b) and
schemes as (4.2), (4.3) are used for chemical potential. The free energy will
dissipate with time as follows:

$\frac{F_{h}\left(\phi^{k+1}\right)-F_{h}\left(\phi^{k}\right)}{\Delta t}\leq
M\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\left(\mu_{i}^{k+1}-\mu_{j}^{k+1}\right)^{2}\omega\left(r_{ij}\right)\leq
0.$ (4.5)

Proof For homogeneous part: By introducing an intermediate function
$c\left(\phi_{i}\right)=\frac{1}{2}\left(\xi+1-\phi_{i}^{2}\right)$ as [26],
we have energy density

$f_{b}\left(\phi_{i}\right)=\lambda_{b}\left(c\left(\phi_{i}\right)^{2}+\frac{1}{2}\xi\phi_{i}^{2}-\frac{1}{2}\xi-\frac{1}{4}\xi^{2}\right).$

A sufficient large $\xi>0$ ensures $\phi_{i}\in
I_{\xi}=[-\sqrt{\xi+1},\sqrt{\xi+1}]$. Obviously, $c\left(\phi_{i}\right)$ is
a non-negative concave function and the quadratic function $\phi_{i}^{2}$ is a
convex function. We note that
$\bar{c}=c\left(\phi_{i}^{k}\right)+c^{\prime}\left(\phi_{i}^{k}\right)\left(\phi_{i}^{k+1}-\phi_{i}^{k}\right)$
and $0\leq c\left(\phi_{i}^{k+1}\right)\leq\bar{c}$ for
$\phi_{i}^{k},\phi_{i}^{k+1}\in I_{\xi}$. The energy density change can be
estimated by

$\displaystyle
f_{b}\left(\phi_{i}^{k+1}\right)-f_{b}\left(\phi_{i}^{k}\right)$
$\displaystyle=\lambda_{b}\left[c\left(\phi_{i}^{k+1}\right)^{2}-c\left(\phi_{i}^{k}\right)^{2}+\frac{1}{2}\xi\left(\left(\phi_{i}^{k+1}\right)^{2}-\left(\phi_{i}^{k}\right)^{2}\right)\right]$
$\displaystyle\leq\lambda_{b}\left[\bar{c}^{2}-c\left(\phi_{i}^{k}\right)^{2}+\xi\phi_{i}^{k+1}\left(\phi_{i}^{k+1}-\phi_{i}^{k}\right)\right]$
$\displaystyle=\lambda_{b}\left[\left(\bar{c}+c\left(\phi_{i}^{k}\right)\right)c^{\prime}\left(\phi_{i}^{k}\right)+\xi\phi_{i}^{k+1}\right]\left(\phi_{i}^{k+1}-\phi_{i}^{k}\right).$

By substituting values of $c^{\prime}\left(\phi_{i}^{k}\right)$ and
$c\left(\phi_{i}^{k}\right)$, the semi-implicit scheme for $\mu_{i}^{b}$
becomes

$\begin{split}\mu_{i}^{b}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)&=\lambda_{b}\left[\left(\bar{c}+c\left(\phi_{i}^{k}\right)\right)c^{\prime}\left(\phi_{i}^{k}\right)+\xi\phi_{i}^{k+1}\right]\\\
&=\lambda_{b}\left[\left(\xi+\left(\phi_{i}^{k}\right)^{2}\right)\phi_{i}^{k+1}-(\xi+1)\phi_{i}^{k}\right].\end{split}$

It naturally yields

$f_{b}\left(\phi_{i}^{k+1}\right)-f_{b}\left(\phi_{i}^{k}\right)\leq
f_{b}^{\prime}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right)\left(\phi_{i}^{k+1}-\phi_{i}^{k}\right).$
(4.6)

Therefore, we conclude that there is always a sufficient large enough $\xi>0$
to make inequality (4.6) hold. For the entire system

$\frac{F_{h}^{b}\left(\phi^{k+1}\right)-F_{h}^{b}\left(\phi^{k}\right)}{\Delta
t}\leq\sum_{i}\frac{m_{i}}{\rho_{i}}\left(f_{b}^{\prime}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right),\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}\right).$ (4.7)

######  Remark 4.1.

This kind of energy-factorization method finally results in a linear energy-
stable discrete scheme for $\mu_{i}^{b}$, which benefits the implementation.
It must be noted that other well-known schemes like SAV and IEQ can be applied
to this part as well. Also, the numerical value of $\phi_{i}$ can be bounded,
which validates our conclusion. The relevant proof can be referred to [33, 34,
35].

For inhomogeneous part: We turn to the scheme of chemical potential
$\mu_{\nabla}$ caused by inhomogeneity. Since the corresponding energy
$F_{\nabla}$ is not a local function, it depends not only on the $\phi$ value
but also on the position $\mathbf{x}$ of all particles. In the current step
for solving the CH system, positions of the particles will not change. An
intermediate variable $F_{i,h}^{\nabla}$ for a single particle is defined for
our proof as follows,

$F_{i,h}^{\nabla}=\frac{\lambda}{2}\left(\nabla_{h}\phi_{i}\right)\cdot\left(\nabla_{h}\phi_{i}\right).$
(4.8)

The free energy caused by the inhomogeneity of
system$F_{h}^{\nabla}=\sum_{i}\frac{m_{i}}{\rho_{i}}F_{i,h}^{\nabla}$.Since
$F_{i,h}^{\nabla}$ is not a local function as well, the convex-concave
property is related to its Hessian matrix. Taking the first order derivative
with respect to,

$\frac{\partial
F_{i,h}^{\nabla}}{\partial\phi_{m}}=\lambda\sum_{j}\frac{m_{i}}{\rho_{j}}\left(\delta_{im}-\delta_{jm}\right)\nabla_{i}W_{ij}\sum_{k}\frac{m_{i}}{\rho_{k}}\left(\phi_{i}-\phi_{k}\right)\nabla_{i}W_{ik}.$

Then, taking the derivative with respect to $\phi_{n}$, we have elements of
its Hessian matrix,

$\frac{\partial^{2}F_{i,h}^{\nabla}}{\partial\phi_{m}\phi_{n}}=\lambda\sum_{m}\sum_{j}\frac{m_{i}}{\rho_{j}}\left(\delta_{im}-\delta_{jm}\right)\nabla_{i}W_{ij}\sum_{n}\sum_{k}\frac{m_{i}}{\rho_{k}}\left(\delta_{in}-\delta_{kn}\right)\nabla_{i}W_{ik}.$

For any vector $\mathbb{b}\neq\mathbb{0}$, we have

$\sum_{m}\sum_{n}b_{m}\frac{\partial^{2}F_{i,h}^{\nabla}}{\partial\phi_{m}\phi_{n}}b_{n}=\lambda\left(\sum_{n}\sum_{k}b_{n}\frac{m_{i}}{\rho_{k}}\left(\delta_{in}-\delta_{kn}\right)\nabla_{i}W_{ik}\right)^{2}\geq
0.$ (4.9)

The Hessian matrix
$\frac{\partial^{2}F_{i,h}^{\nabla}}{\partial\phi_{m}\partial\phi_{n}}$ is
positive-definition. Thus, $F_{i,h}^{\nabla}$ is convex with respect to
$\phi$. If the chemical potential $\mu_{i}^{\nabla}$ is treated with a fully
implicit scheme as (4.3), it will make the energy inequality hold as

$\frac{F_{h}^{\nabla}\left(\phi^{k+1}\right)-F_{h}^{\nabla}\left(\phi^{k}\right)}{\Delta
t}\leq\sum_{i}\frac{m_{i}}{\rho_{i}}\sum_{j}\left(\frac{\partial
F_{i,h}^{\nabla}\left(\phi^{k+1}\right)}{\partial\phi_{j}^{k+1}},\frac{\phi_{j}^{k+1}-\phi_{j}^{k}}{\Delta
t}\right).$ (4.10)

Since we consider the incompressible fluid, it yields the following identity:

$\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}\left(\phi^{k+1}\right)}{\partial\phi_{i}^{k+1}},\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}\right)=\sum_{i}\frac{m_{i}}{\rho_{i}}\sum_{j}\left(\frac{\partial
F_{i,h}^{\nabla}\left(\phi^{k+1}\right)}{\partial\phi_{j}^{k+1}},\frac{\phi_{j}^{k+1}-\phi_{j}^{k}}{\Delta
t}\right).$

Then we find the inequality

$\frac{F_{h}^{\nabla}\left(\phi^{k+1}\right)-F_{h}^{\nabla}\left(\phi^{k}\right)}{\Delta
t}\leq\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}\left(\phi^{k+1}\right)}{\partial\phi_{i}^{k+1}},\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}\right).$ (4.11)

Similarly, we use the inner product to get the energy dissipation law at the
fully discrete level. By taking the inner product of equation (4.1a) with
$\mu_{i}^{k+1}$

$\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t},\mu_{i}^{k+1}\right)=M\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}\rho_{j}}\left(\mu_{i}^{k+1}-\mu_{j}^{k+1}\right)^{2}\omega\left(r_{ij}\right).$
(4.12)

By taking the inner product of (4.1b) with
$\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta t}$, we have

$\begin{split}\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\mu_{i}^{k+1},\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}\right)=&\sum_{i}\frac{m_{i}}{\rho_{i}}\left(\frac{\rho_{i}}{m_{i}}\frac{\partial
F_{h}^{\nabla}\left(\phi^{k+1}\right)}{\partial\phi_{i}^{k+1}},\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}\right)\\\
&+\sum_{i}\frac{m_{i}}{\rho_{i}}\left(f_{b}^{\prime}\left(\phi_{i}^{k},\phi_{i}^{k+1}\right),\frac{\phi_{i}^{k+1}-\phi_{i}^{k}}{\Delta
t}\right).\end{split}$ (4.13)

Owing to
$F_{h}\left(\phi\right)=F_{h}^{b}\left(\phi\right)+F_{h}^{\nabla}\left(\phi\right)$,
with the above equations (4.12), (4.13) and energy inequalities (4.7), (4.11)
derived from the homogeneous and inhomogeneous parts, we get the energy
dissipation law of the CH system as (4.5)

######  Lemma 4.2.

Assume the fully discrete scheme for NS system as (4.1d), (4.1e), (4.1f),
(4.1g) is used. Among them, a stabilizing semi-implicit scheme is applied to
the interfacial term as (4.4). Then a sufficient great $\mathcal{M}$ always
makes the energy dissipation law of the NS system at the full discrete level
hold well.

$\frac{E_{k,h}\left(\mathbf{v}^{k+1}\right)-E_{k,h}\left(\mathbf{v}^{k}\right)}{\Delta
t}+\frac{F_{h}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)}{\Delta
t}\leq 0.$ (4.14)

Proof It has been mentioned that the change of energy $F^{\nabla}_{h}$ is
partially caused by the change of the position of particles $\mathbf{x}$. When
the position of particles is updated through (4.1g) and $\phi$ is unvarying,
this part of energy change characterizes the energy transfer from kinetic
energy $E_{k}$ to the free energy $F$. This relationship is consistent with
the definition in the ODE system. Using the symbol $\langle\rangle_{h}^{tem}$
to denote the fully discrete scheme, including temporal discretization. Thus,

$\langle\text{transfer}_{\mathrm{EK-F}}\rangle_{h}^{tem}=\frac{F_{h}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)}{\Delta
t}.$ (4.15)

Correspondingly, the energy transferred from free energy $F$ to kinetic energy
$E_{k}$ is characterized by the inner product of interfacial term with
$m_{i}\tilde{\mathbf{v}}_{i}^{k+1}$.

$\begin{split}\left\langle\operatorname{transfer}_{\mathrm{F}-\mathrm{EK}}\right\rangle_{h}^{\mathrm{tem}}&=-\sum_{i}\left(m_{i}\tilde{\mathbf{v}}_{i}^{k+1}\cdot\frac{\left\langle\mathcal{T}_{i}\left(\tilde{\mathbf{v}}^{k+1},\mathbf{v}^{k}\right)\right\rangle_{h}}{m_{i}}\right)\\\
&=-\sum_{i}\tilde{\mathbf{v}}_{i}^{k+1}\cdot\left\langle\mathcal{T}_{i}\left(\tilde{\mathbf{x}}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}\approx-\sum_{i}\left\langle\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}\frac{\partial\mathbf{x}_{i}}{\partial
t}\right\rangle_{h}.\end{split}$ (4.16)

Since the NS system is determined by projection methods, we split the entire
process into three steps: (1) solving the momentum equation; (2) solving the
Poisson equation; (3) particle movement. So the energy change can be expressed
as the summation of

$\frac{\left[E_{k,h}\left(\tilde{\mathbf{v}}^{k+1}\right)-E_{k,h}\left(\mathbf{v}^{k}\right)\right]}{\Delta
t},\frac{\left[E_{k,h}\left(\mathbf{v}^{k+1}\right)-E_{k,h}\left(\tilde{\mathbf{v}}^{k+1}\right)\right]}{\Delta
t}$, and
$\frac{\left[F_{h}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)\right]}{\Delta
t}$.

With the help of inner product operation to get the energy and based on the
inequality
$\frac{1}{2}\|\mathbf{a}\|^{2}-\frac{1}{2}\|\mathbf{b}\|^{2}\leq\mathbf{a}\cdot(\mathbf{a}-\mathbf{b})$,
we can get the energy relations for each step.

(1) momentum equation: By multiplying $m_{i}\tilde{\mathbf{v}}_{i}^{k+1}$ with
the momentum equation (4.1d), the left-hand side (LHS) becomes

$E_{k,h}\left(\tilde{\mathbf{v}}^{k+1}\right)-E_{k,h}\left(\mathbf{v}^{k}\right)=\frac{1}{2}\sum_{i}m_{i}\left\|\tilde{\mathbf{v}}_{i}^{k+1}\right\|^{2}-\frac{1}{2}\sum_{i}m_{i}\left\|\mathbf{v}_{i}^{k}\right\|^{2}\leq\sum_{i}m_{i}\tilde{\mathbf{v}}_{i}^{k+1}\cdot\left(\tilde{\mathbf{v}}_{i}^{k+1}-\mathbf{v}_{i}^{k}\right).$

Then, combining (4.16), it can concluded that

$\begin{split}\frac{E_{k,h}\left(\tilde{\mathbf{v}}^{k+1}\right)-E_{k,h}\left(\mathbf{v}^{k}\right)}{\Delta
t}&\leq\frac{(d+2)}{2}\sum_{i}\sum_{j}\frac{m_{i}m_{j}}{\rho_{i}^{k+1}\rho_{j}^{k+1}}\frac{\left(\eta_{i}+\eta_{j}\right)\left\|\tilde{\mathbf{v}}_{ij}^{k+1}\cdot\mathbf{x}_{ij}^{k}\right\|^{2}}{\left\|\mathbf{x}_{ij}^{k}\right\|^{2}}\omega\left(r_{ij}\right)\\\
&+\langle\text{transfer}_{\mathrm{F-EK}}\rangle_{h}^{tem}.\end{split}$ (4.17)

(2) Poisson equation: By multiplying $m_{i}\mathbf{v}_{i}^{k+1}$ with the
Poisson equation (4.1e), we deduce

$\sum_{i}\left(m_{i}\mathbf{v}_{i}^{k+1},\frac{\mathbf{v}_{i}^{k+1}-\tilde{\mathbf{v}}_{i}^{k+1}}{\Delta
t}\right)=\sum_{i}\left(m_{i}\mathbf{v}_{i}^{k+1},-\sum_{j}m_{j}\left(\frac{p_{i}^{k+1}}{\left(\rho_{i}^{k+1}\right)^{2}}+\frac{p_{j}^{k+1}}{\left(\rho_{j}^{k+1}\right)^{2}}\right)\nabla_{i}W_{ij}^{n}\right).$
(4.18)

Based on the divergence-free condition at the discrete level, it naturally
yields

$-\sum_{i}m_{i}\frac{p_{i}^{k+1}}{\left(\rho_{i}^{k+1}\right)^{2}}\sum_{j}m_{j}\left(\mathbf{v}_{i}^{k+1}-\mathbf{v}_{j}^{k+1}\right)\cdot\nabla_{i}W_{ij}=0.$

Then, we can conclude the right-hand side (RHS) of (4.18)

$\text{ RHS
}=-\frac{1}{2}\sum_{i}\sum_{i}m_{i}m_{j}\left(\frac{p_{i}^{k+1}}{\left(\rho_{i}^{k+1}\right)^{2}}+\frac{p_{j}^{k+1}}{\left(\rho_{j}^{k+1}\right)^{2}}\right)\left(\mathbf{v}_{i}^{k+1}-\mathbf{v}_{j}^{k+1}\right)\cdot\nabla_{i}W_{ij}=0.$

Finally, we get

$\frac{E_{k,h}\left(\mathbf{v}^{k+1}\right)-E_{k,h}\left(\tilde{\mathbf{v}}_{i}^{k+1}\right)}{\Delta
t}\leq\sum_{i}\left(m_{i}\mathbf{v}_{i}^{k+1},\frac{\mathbf{v}_{i}^{k+1}-\tilde{\mathbf{v}}_{i}^{k+1}}{\Delta
t}\right)=0.$ (4.19)

(3) movement of particles: Based on the definition of interfacial force:
$\boldsymbol{\mathcal{T}}_{i}:=\frac{\partial
F_{h}^{\nabla}}{\partial\mathbf{x}_{i}}$ and the definition of
$F_{i,h}^{\nabla}$, for $\mathbf{x}_{i}^{k+1}$, $\mathbf{x}_{i}^{k}\in S$ and
the open set $S\subseteq\mathbb{R}^{d}(d=2,3)$, let
$\mathbf{x}_{i}=(1-\alpha_{i})\mathbf{x}_{i}^{k}+\alpha_{i}\mathbf{x}_{i}^{k+1}$,
$\alpha_{i}\in[0,1]$, and $\mathbf{x}$ is the collection of $\mathbf{x}_{i}$.
If we consider the movement as a particle-wise process and based on the mean
value theorem for multivariable, we have

$F_{h}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)=\sum_{i}\boldsymbol{\mathcal{T}}_{i}(\mathbf{x})\cdot\left(\mathbf{x}_{i}^{k+1}-\mathbf{x}_{i}^{k}\right).$

Through the convex-concave splitting, a semi-implicit scheme can make the
following inequality hold:

$\boldsymbol{\mathcal{T}}_{i}(\mathbf{x})\cdot\left(\mathbf{x}_{i}^{k+1}-\mathbf{x}_{i}^{k}\right)\leq\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}\cdot\left(\mathbf{x}_{i}^{k+1}-\mathbf{x}_{i}^{k}\right).$

namely,

$\langle\text{transfer}_{\mathrm{EK-F}}\rangle_{h}^{tem}=\frac{F_{h}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)}{\Delta
t}\leq\sum_{i}\mathbf{v}_{i}^{k+1}\cdot\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}.$
(4.20)

By employing a stabilizing term, the energy can be rewritten in a modified
form: $F_{h}^{\nabla}=F_{h,e}^{\nabla}+F_{h,c}^{\nabla}$, where subscript $e$
means convex and $c$ means concave. Then the original energy $F_{h}^{\nabla}$
can be split into a convex part and a concave part:

$\left\\{\begin{array}[]{l}F_{h,e}^{\nabla}=\frac{\mathcal{M}}{\Delta
t}\sum_{i}\left\|\mathbf{x}_{i}\right\|^{2},\\\
F_{h,c}^{\nabla}=F_{h}^{\nabla}-\frac{\mathcal{M}}{\Delta
t}\sum_{i}\left\|\mathbf{x}_{i}\right\|^{2}.\end{array}\right.$

Their second order derivatives are

$\frac{\partial^{2}F_{h,e}^{\nabla}}{\partial\mathbf{x}_{i}\partial\mathbf{x}_{j}}=\left\\{\begin{array}[]{ll}2\frac{\mathcal{M}}{\Delta
t}&,i=j,\\\ 0&,i\neq j,\end{array}\quad\text{ and
}\quad\frac{\partial^{2}F_{h,c}^{\nabla}}{\partial\mathbf{x}_{i}\partial\mathbf{x}_{j}}=\begin{cases}-2\frac{\mathcal{M}}{\Delta
t}+\frac{\partial\mathcal{T}_{i}}{\partial\mathbf{x}_{i}}&,i=j,\\\
\frac{\partial\mathcal{T}_{i}}{\partial\mathbf{x}_{j}}&,i\neq
j.\end{cases}\right.$

The Hessian matrix of $F_{h,e}^{\nabla}$ is always positive-definite.
Moreover, a sufficiently large modified parameter $\mathcal{M}$ always exists
to make the Hessian matrix of $F_{h,c}^{\nabla}$ negative-definite. Based on
the convex-concave properties, this splitting makes the following inequalities
hold:

$\left\\{\begin{array}[]{l}F_{h,e}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h,e}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)\leq\sum_{i}2\mathcal{M}\mathbf{x}_{i}^{k+1}\cdot\left(\mathbf{x}_{i}^{k+1}-\mathbf{x}_{i}^{k}\right),\\\
F_{h,c}^{\nabla}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h,c}^{\nabla}\left(\mathbf{x}^{k},\phi^{k+1}\right)\leq\sum_{i}\left(\mathcal{T}_{i}\left(\mathbf{x}^{k}\right)-2\mathcal{M}\mathbf{x}_{i}^{k}\right)\cdot\left(\mathbf{x}_{i}^{k+1}-\mathbf{x}_{i}^{k}\right).\end{array}\right.$
(4.21)

Combining these two inequalities (4.21), we find the appropriate semi-implicit
fully discrete scheme for interfacial force

$\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}=\frac{2\mathcal{M}}{\Delta
t}\mathbf{x}_{i}^{k+1}+\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k}\right)-\frac{2\mathcal{M}}{\Delta
t}\mathbf{x}_{i}^{k}=2\mathcal{M}\mathbf{v}_{i}^{k+1}+\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k}\right).$
(4.22)

The scheme for $\boldsymbol{\mathcal{T}}_{i}$ as (4.22) will ensure inequality
(4.20). For the physical consistency of the full discrete scheme, we apply the
same scheme for the term
$\left\langle\boldsymbol{\mathcal{T}}_{i}\left(\tilde{\mathbf{x}}^{k+1},\mathbf{x}^{k}\right)\right\rangle_{h}$
with the intermediate velocity $\tilde{\mathbf{v}}_{i}^{k+1}$ in the momentum
equation as (4.4). Thus,

$\left\langle\text{transfer}_{\mathrm{F-EK}}\right\rangle_{h}^{tem}=-\sum_{i}\tilde{\mathbf{v}}_{i}^{k+1}\cdot\left[2M\tilde{\mathbf{v}}_{i}^{k+1}+\boldsymbol{\mathcal{T}}_{i}\left(\mathbf{x}^{k}\right)\right].$

By recalling the inequality (4.19) derived from the pressure Poisson equation,
we can conclude that there is always a sufficiently large coefficient
$\mathcal{M}$ to make

$\displaystyle\quad\langle\text{transfer}_{\mathrm{EK-F}}\rangle_{h}^{tem}+\langle\text{transfer}_{\mathrm{F-EK}}\rangle_{h}^{tem}$
(4.23)
$\displaystyle\leq\sum_{i}\mathbf{v}_{i}^{k+1}\cdot\left[2\mathcal{M}\mathbf{v}_{i}^{k+1}+\mathcal{T}_{i}\left(\mathbf{x}^{k}\right)\right]-\sum_{i}\tilde{\mathbf{v}}_{i}^{k+1}\cdot\left[2\mathcal{M}\tilde{\mathbf{v}}_{i}^{k+1}+\mathcal{T}_{i}\left(\mathbf{x}^{k}\right)\right]$
$\displaystyle=2\mathcal{M}\sum_{i}\left[\left(\mathbf{v}_{i}^{k+1}\right)^{2}-\left(\tilde{\mathbf{v}}_{i}^{k+1}\right)^{2}\right]+\sum_{i}\left(\mathbf{v}_{i}^{k+1}-\tilde{\mathbf{v}}_{i}^{k+1}\right)\mathcal{T}_{i}\left(\mathbf{x}^{k}\right)\leq
0.$

At last, combining the energy changes at three steps as a (4.17), (4.19), and
(4.23), we reach the conclusion of the lemma as (4.14).

######  Remark 4.2.

We could consider the conservation of energy transfer still maintained at the
fully discrete level. But the energy inequality of (4.23) is equivalent to
adding a numerical dissipative term to help maintain the energy decaying
property. It comes from the numerical approximation of energy transfer.

$\left\langle\operatorname{transfer}_{diss}\right\rangle_{h}^{\mathrm{tem}}=2\mathcal{M}\sum_{i}\left[\left(\mathbf{v}_{i}^{k+1}\right)^{2}-\left(\tilde{\mathbf{v}}_{i}^{k+1}\right)^{2}\right]+\sum_{i}\left(\mathbf{v}_{i}^{k+1}-\tilde{\mathbf{v}}_{i}^{k+1}\right)\mathcal{T}_{i}\left(\mathbf{x}^{k}\right)\leq
0.$ (4.24)

The value of $\mathcal{M}$ is adjusted adaptively posteriori based on the
solved $\mathbf{v}_{i}^{k+1}$. If the above energy inequality (4.24) does not
hold, the value of $\mathcal{M}$ can be appropriately increased to maintain
the energy decaying trend.

Finally, the theorem for the complete NSCH system at the fully discrete level
is proposed here.

######  Theorem 4.1.

Assume the numerical scheme for the NSCH system based on the SPH method
includes (4.1a), (4.1b), (4.1c), (4.1d), (4.1e), (4.1f) and (4.1g). The energy
dissipation law holds at the fully discrete level.

$\frac{E_{\text{total},h}^{k+1}-E_{\text{total},h}^{k}}{\Delta
t}=\frac{E_{k,h}\left(\mathbf{v}^{k+1}\right)-E_{k,h}\left(\mathbf{v}^{k}\right)}{\Delta
t}+\frac{F_{h}\left(\mathbf{x}^{k+1},\phi^{k+1}\right)-F_{h}\left(\mathbf{x}^{k},\phi^{k}\right)}{\Delta
t}\leq 0.$ (4.25)

Proof This energy inequality for the NSCH system can be concluded through
adding these two inequalities, (4.5) and (4.14), from Lemmas 4.1 and 4.2 for
the CH and NS systems, respectively.

##  5 Numerical examples

In this section, we present some numerical results to show the performance of
the proposed energy-stable SPH method for the NSCH two-phase system. We
consider the evolution of a droplets with the interfacial tension. In all of
the numerical cases, the spatial domain is taken as $\Omega=[-1,1]^{2}$. The
particles are placed uniformly in the initial step, and there are $40\times
40$ particles in all examples. The smoothing length is set as $h$.

###  5.1 The dynamics of a square droplet

In the first numerical example, we consider the evolution process of a square
droplet in a domain. The parameters are chosen as: $h=0.05,\Delta
t=0.01,T=10,\mathcal{M}=10^{3},\eta=1,\lambda=1,\epsilon=0.02,\text{Mobility}=0.002,\xi=1.$
The initial phase variables are chosen as $\phi=1$ inside the square droplet
and $\phi=-1$ everywhere else.

(a) Initial

(b) $t=1$

(c) $t=20$

Figure 5.2: The evolution of $\phi$ and spatial particle distribution with
time in example 1

(a) Initial

(b) $t=1$

(c) $t=20$

Figure 5.3: The contour of $\phi$ interpolated by the SPH operator

(a) Energy

(b) Momentum

Figure 5.4: Total energy profile and momentum profile with time

In Fig.5.2, we observe that the droplet in a square shape deforms into a round
one under the interfacial tension. The distribution of particles in the
eventual steady state is totally different from the initial uniform
distribution. But generally, the spacing between particles remains uniform.
So, we conclude that the SPH method can maintain the divergence-free condition
well and does not cause tensile instability. The symmetric velocity field and
the Fig.5.4(b) also imply momentum conservation at the ODE level and the
discrete level. We find the momentum oscillates a little bit numerically
around zero from the zoom window. Particles at the interface are also
distributed one round layer by one round layer, which is caused by the
interfacial force defined by the gradient of $\phi$. While in the bulk phase
part, the particles are located more randomly and messily, which makes sense
because of the homogeneity.

In order to show the diffuse characteristics of the NSCH model, an initial
condition with a sharp interface is set here. The results show the forming of
a diffuse interface with limited width as the numerical calculation process
progresses. Finally, the “color” of particles on the diffuse interface region
also obeys a gradual transition from $-1$ to $1$. This phenomenon coincides
with the physical principles. Furthermore, the contour of $\phi$ at the
continuous domain interpolated by the SPH in Fig.5.3 shows that it agrees well
with the results of the traditional method based on the fine mesh.

Moreover, we observe the consistent decay of the total energy in the
Fig.5.4(a), which validates our design of ODE expression and the fully
discrete scheme. From the zoom window, we can see total energy is dissipating
with time, even at the end of the stage of numerical testing. Gradually, when
this droplet becomes one in a round shape, the system reaches a steady state
with the minimum energy.

###  5.2 The deformation and rotation of a droplet in the force field

We consider the case of a single droplet’s deformation and rotation phenomena
in the square domain. At the initial stage, a force field $g(x,y,t)$ is
enforced in this periodic domain to make this droplet rotate and deform,

$g(x,y,t)=\left\\{\begin{array}[]{l}10\sin(\pi(x+1)),\quad 0\leq t\leq 0.5,\\\
0,\quad\text{otherwise}.\end{array}\right.$

The parameters are chosen as: $h=0.05,\,\Delta
t=0.01,\,T=20,\,\mathcal{M}=10^{3},\,\eta=1,\lambda=1\,,\,\epsilon=0.02,\,\text{Mobility}=0.002,\xi=1.$
Initially, the diffusive interface is smoothed by
$\phi_{0}(x,y)=-\tanh((x^{2}+y^{2}-0.3^{2})/0.02)$

Figure 5.5: The diagram of force field

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=20$

Figure 5.6: The evolution of $\phi$ and spatial particle distribution with
time in example 2

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=20$

Figure 5.7: The evolution of velocity quiver plot

(a) Energy

(b) Momentum

Figure 5.8: Total energy profile and momentum profile with time

The diffuse interface we set up artificially as the initial conditional
gradually evolves into an interface with more physical transition layers. This
example shows the good adaptivity to the two-phase flow problems with high
deformation. The interface can be tracked naturally with the $\phi$
information carried by the particles as Fig.5.6. The symmetric velocity quiver
plots as Fig.5.7 and Fig.5.8(b) indicate the maintenance of divergence-free
conditions and momentum conservation.

In Fig.5.8(a), the total energy of the system is increasing at the beginning,
in the time period $0\textendash 0.5$, due to the external work done by the
force field. As we expected, after cancelling this force field at time =
$0.5$, the total energy of the system decays with time due to the viscous
dissipation and the diffusive effect of the CH system and the momentum
conservation holds well even facing relative high flow rate at the initial
stage as well.

###  5.3 Droplets of different sizes merge into one

In this example, we consider the case of four circular droplets merging with
one large droplet in the computational domain. The parameters are chosen as:
$h=0.05,\Delta
t=0.01,T=10,\mathcal{M}=10^{3},\eta=1,\lambda=1,\epsilon=0.02,\text{Mobility}=0.002,\xi=1.$
The initial phase variable is chosen as:

$\displaystyle\phi_{0}(x,y)$
$\displaystyle=-\tanh((x^{2}+y^{2}-0.3^{2})/0.01)*$
$\displaystyle\tanh(((x-0.4)^{2}+y^{2}-0.1^{2})/0.01)*\tanh(((x-0.2)^{2}+y^{2}-0.1^{2})/0.01)*$
$\displaystyle\tanh(((y-0.2)^{2}+x^{2}-0.1^{2})/0.01)*\tanh(((y-0.4)^{2}+x^{2}-0.1^{2})/0.01).$

The free energy contours in Fig.5.10 show that the maximum free energy
decreases from $1.5$ to $0.3$. At the initial stage, the interfaces of
droplets with different sizes connect with each other with very sharp angles.
Since the free mixing energy characterizes the surface intension, the sharp
connection naturally leads to a higher free energy in this region. The maximum
value of the free energy gradually decreases as the angle becomes round and
smooth. The homogenization of $\phi$ at the interface region also helps
decrease the maximum value of free energy. The energy dissipation and momentum
conservation are preserved well in example 3 as Fig.5.11.

(a) Initial

(b) $t=1$

(c) $t=30$

Figure 5.9: The evolution of $\phi$ and spatial particle distribution with
time in example 3

(a) Initial

(b) $t=1$

(c) $t=30$

Figure 5.10: The evolution of free energy throughout the entire domain

(a) Energy

(b) Momentum

Figure 5.11: Total energy profile and momentum profile with time

##  6 Conclusion

In this work, we exploit the great potential of the Lagrangian particle
treatment in the SPH, which also fulfills the increasing demand for particle
modeling that allows more physical properties to be incorporated. To the best
of our knowledge, this is the first study of energy-stable ODE discretization
and the energy-stable fully discrete scheme by the SPH method for two-phase
problems. Our proposed scheme ensures the inheritance of momentum conservation
and the energy dissipation law from the PDE level to the ODE level, and then
to the fully discrete level. Consequently and desirably, it also helps
increase the stability of the numerical method. The time step size can be much
larger than that of the traditional ISPH methods. This energy-stable SPH
method also alleviates the tensile instability without using any particle-
shifting strategies, which may destroy the rigorous mathematical proof. The
numerical results also demonstrate that our method captures the interface
behavior and the energy variation process well.

## Acknowledgments

This work is supported by King Abdullah University of Science and Technology
(KAUST) through the grants BAS/1/1351-01, URF/1/407401, and URF/1/3769-01. Z.
Qiao’s work is partially supported by the Hong Kong Research Grants Council
(RFS Project No.RFS2021-5S03 and GRF project No.15302919) and the Hong Kong
Polytechnic University internal grant No.4-ZZLS.

## References

  * [1] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Monthly Notices of the Royal Astronomical Society 181 (3) (1977) 375–389.
  * [2] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The Astronomical Journal 82 (1977) 1013–1024.
  * [3] X. Hu, N. A. Adams, A constant-density approach for incompressible multi-phase SPH, Journal of Computational Physics 228 (6) (2009) 2082–2091.
  * [4] A. M. Tartakovsky, N. Trask, K. Pan, B. Jones, W. Pan, J. R. Williams, Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media, Computational Geosciences 20 (4) (2016) 807–834.
  * [5] Z. Chen, L. Shen, A modified smoothed particle hydrodynamics for modelling fluid-fracture interaction at mesoscale, Computational Particle Mechanics (2021) 1–21.
  * [6] V. Zago, G. Bilotta, A. Hérault, R. A. Dalrymple, L. Fortuna, A. Cappello, G. Ganci, C. Del Negro, Semi-implicit 3D SPH on GPU for Lava flows, Journal of Computational Physics 375 (2018) 854–870.
  * [7] D. J. Price, Smoothed particle hydrodynamics and magnetohydrodynamics, Journal of Computational Physics 231 (3) (2012) 759–794.
  * [8] M. Nakajima, D. J. Stevenson, Melting and mixing states of the Earth’s mantle after the Moon-forming impact, Earth and Planetary Science Letters 427 (2015) 286–295.
  * [9] M. Ihmsen, J. Orthmann, B. Solenthaler, A. Kolb, M. Teschner, SPH fluids in computer graphics, in: EUROGRAPHICS 2014/S. LEFEBVRE AND M. SPAGNUOLO, Citeseer, 2014.
  * [10] J. Wang, G. Liu, A point interpolation meshless method based on radial basis functions, International Journal for Numerical Methods in Engineering 54 (11) (2002) 1623–1648.
  * [11] W. G. Hoover, Isomorphism linking smooth particles and embedded atoms, Physica A: Statistical Mechanics and its Applications 260 (3-4) (1998) 244–254.
  * [12] D. Nishiura, M. Furuichi, H. Sakaguchi, Computational performance of a smoothed particle hydrodynamics simulation for shared-memory parallel computing, Computer Physics Communications 194 (2015) 18–32.
  * [13] J. J. Monaghan, Smoothed particle hydrodynamics, Reports on Progress in Physics 68 (8) (2005) 1703.
  * [14] M. Liu, G. Liu, Smoothed particle hydrodynamics (SPH): an overview and recent developments, Archives of Computational Methods in Engineering 17 (1) (2010) 25–76.
  * [15] A. Khayyer, H. Gotoh, Y. Shimizu, Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context, Journal of Computational Physics 332 (2017) 236–256.
  * [16] R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach, Journal of Computational Physics 228 (18) (2009) 6703–6725.
  * [17] X. Zhu, S. Sun, J. Kou, An energy stable sph method for incompressible fluid flow, Advances in Applied Mathematics and Mechanics 14 (5) (2022) 1201–1224.
  * [18] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics 79 (1) (1988) 12–49.
  * [19] C. W. Hirt, B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1) (1981) 201–225.
  * [20] D. Sun, W. Tao, A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows, International Journal of Heat and Mass Transfer 53 (4) (2010) 645–655.
  * [21] J. W. Gibbs, On the equilibrium of heterogeneous substances, American Journal of Science 3 (96) (1878) 441–458.
  * [22] D. J. Eyre, Unconditionally gradient stable time marching the Cahn-Hilliard equation, MRS Online Proceedings Library (OPL) 529 (1998).
  * [23] X. Feng, J. Kou, S. Sun, A novel energy stable numerical scheme for Navier-Stokes-Cahn-Hilliard two-phase flow model with variable densities and viscosities, in: International Conference on Computational Science, Springer, 2018, pp. 113–128.
  * [24] J. Shen, J. Xu, J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, Journal of Computational Physics 353 (2018) 407–416.
  * [25] Q. Du, L. Ju, X. Li, Z. Qiao, Maximum principle preserving exponential time differencing schemes for the nonlocal Allen–Cahn equation, SIAM Journal on Numerical Analysis 57 (2) (2019) 875–898.
  * [26] X. Wang, J. Kou, H. Gao, Linear energy stable and maximum principle preserving semi-implicit scheme for Allen–Cahn equation with double well potential, Communications in Nonlinear Science and Numerical Simulation 98 (2021) 105766\.
  * [27] X. Feng, M.-H. Chen, Y. Wu, S. Sun, A fully explicit and unconditionally energy-stable scheme for peng-robinson vt flash calculation based on dynamic modeling, Journal of Computational Physics (2022) 111275.
  * [28] Z. Qiao, S. Sun, Two-phase fluid simulation using a diffuse interface model with Peng–Robinson equation of state, SIAM Journal on Scientific Computing 36 (4) (2014) B708–B728.
  * [29] J. Kou, S. Sun, X. Wang, A novel energy factorization approach for the diffuse-interface model with Peng–Robinson equation of state, SIAM Journal on Scientific Computing 42 (1) (2020) B30–B56.
  * [30] X. Fan, J. Kou, Z. Qiao, S. Sun, A componentwise convex splitting scheme for diffuse interface models with Van der Waals and Peng–Robinson equations of state, SIAM Journal on Scientific Computing 39 (1) (2017) B1–B28.
  * [31] Q. Yang, J. Yao, Z. Huang, M. Asif, A comprehensive SPH model for three-dimensional multiphase interface simulation, Computers & Fluids 187 (2019) 98–106.
  * [32] M. Hirschler, M. Huber, W. Säckel, P. Kunz, U. Nieken, An application of the Cahn-Hilliard approach to smoothed particle hydrodynamics, Mathematical Problems in Engineering 2014 (2014).
  * [33] D. Li, Z. Qiao, T. Tang, Characterizing the stabilization size for semi-implicit Fourier-spectral method to phase field equations, SIAM Journal on Numerical Analysis 54 (3) (2016) 1653–1681.
  * [34] D. Li, Z. Qiao, On second order semi-implicit fourier spectral methods for 2D Cahn–Hilliard equations, Journal of Scientific Computing 70 (1) (2017) 301–341.
  * [35] X. Li, Z. Qiao, C. Wang, Convergence analysis for a stabilized linear semi-implicit numerical scheme for the nonlocal Cahn–Hilliard equation, Mathematics of Computation 90 (327) (2021) 171–188.