# A 2D hybrid method for interfacial transport of passive scalars

Yu Fan  Yujie Zhu  Xiaoliang Li  Xiangyu Hu  Nikolaus A. Adams  Technical
University of Munich Boltzmannstr. 15 D-85748 Garching  Department of
Mathematics ”Federigo Enriques”, Università degli Studi di Milano, Via Cesare
Saldini, 50, 20133, Milan, Italy

###### Abstract

A hybrid Eulerian-Lagrangian method is proposed to simulate passive scalar
transport on arbitrary shape interface. In this method, interface deformation
is tracked by an Eulerian method while the transport of the passive scalar on
the material interface is solved by a single-layer Lagrangian particle method.
To avoid particle clustering, a novel remeshing approach is proposed. This
remeshing method can resample particles, adjust the position of particles by a
relaxation process, and transfer mass from pre-existing particles to resampled
particles via a redistribution process, which preserves mass both globally and
locally. Computational costs are controlled by an adaptive remeshing strategy.
Accuracy is assessed by a series of test cases.

###### keywords:

(regional) level set method , smoothed particle hydrodynamics , remeshing
method

††journal: Journal of LaTeX Templatesmytitlenotemytitlenotefootnotetext: Fully
documented templates are available in the elsarticle package on CTAN.

##  1 Introduction

The presence of surfactants on interfaces in multi-phase flow can alter the
surface tension, interface shapes and topology [1], and affect the relevant
engineering processes [2, 3] . Numerical simulation is one of the most
important ways to understand surfactant dynamics. In simulations, usually,
surfactants are modeled as passive scalars. The computation of the transport
of passive scalars is a prerequisite for further study of surfactant dynamics.
Currently, there are several numerical approaches to compute passive scalars:
Eulerian, Lagrangian, and hybrid. However, these approaches encounter
difficulties with respect to the accuracy, conservation and efficiency.

In the Eulerian approach, the passive scalar on the interface is usually
represented on a grid [4, 5, 6, 7, 8, 9], which causes a mismatch between the
grid data and the off-grid interface, and eventually violates mass
conservation. To cope with this difficulty, a global correction was introduced
[7], but still failed to satisfy mass conservation locally. Another approach
proposed for insoluble surfactants based on a volume-of-fluid method tracks
mass and volume separately, which enforces the conservation of former [10].
However, for large deformation or long-time simulation, the consistency error
between mass and volume can increase considerably.

Contrary to Eulerian methods, Lagrangian methods maintain mass conservation
inherently. A typical example is Smoothed Particle Hydrodynamics (SPH), which
has been used to study surfactant dynamics on bubble surfaces [11]. The
interface was represented by multiple layers of particles so that a 2D problem
is transformed into 3D. Since a 3D problem involves more neighbor connections
when computing physical variables, it leads to extra computational costs in
run-time search of the neighbor connections. To save costs, a codimension-1
SPH method, based on differential geometry, was proposed for surfactant
dynamics on thin film and bubbles [12]. However, as the particles are added by
reseeding when the surface is enlarged, mass is not conserved.

In contrast to either pure Eulerian or Lagrangian approaches, many researchers
propose hybrid approaches, which take advantage of both properties. While the
interface is tracked by an Eulerian method, Lagrangian surface meshes are used
to solve the advection-diffusion equation of the passive scalar [13, 14, 15,
16]. The difficulty is that, when the surface mesh is deformed by the flow,
its mesh quality can not be guaranteed. To avoid degeneration of mesh quality,
an artificial tangential velocity [17] is introduced in Refs. [13, 15]. Other
similar methods proposed to solve this problem used either a map between the
initial parameter of the curve and the arc-length parameter of the current
interface [14], or a surface mesh optimization method [18]. Although these
methods can alleviate the problem of mesh quality, this geometric correction
is quite complicated. An additional problem of these hybrid methods is that,
to solve the transport equation on the surface, connectivity is necessary for
all these methods.

In this work, we investigate an alternative hybrid approach, by using the
mesh-free particle Lagrangian formulation to discretize the transport equation
of the passive scalar on a co-dimension 1 interface while tracking the
interface motion by the level-set method. This discretization is on a single
layer of particles, by which the mass conservation is satisfied automatically.
Since the particles are only distributed on the co-dimension 1 interface, 2D
problems are reduced to 1D, every particle has less neighboring connections,
and computational costs are saved significantly. The method does not require
explicit geometry/topology information and higher-order derivatives, which
simplifies implementation and improves robustness. When the interface
experiences strong stretching or compression, a generalized remeshing method
is proposed based on classical remeshing [19, 20]. This generalized remeshing
method can achieve arbitrary order consistency and the mass conservation holds
automatically as the zeroth order consistency. Besides, this method is
generalized for multi-region problems based on the regional level-set method
[21, 22]. The availability of particles can also be used to improve the
accuracy of the regional level-set method without taking too much extra
memory.

The remainder of this paper is organized as follows: Hybrid method and its
generalized version for multi-region problems are detailed in Section 2.
Section 3 comprehensively validates the accuracy, efficiency and conservation
through several benchmark tests.

##  2 The hybrid method

In this section, details of the hybrid Eulerian-Lagrangian method are
presented. We adopt the conservative interface method [23] to track the
interface deformation, while a Lagrangian particle method is used to compute
the evolution of passive scalar on the interface, see Fig. 1. Particles are
generated on the interface based on the level-set field and regularized by a
relaxation process. At each time step, the position of each particle is
updated by the velocity interpolated from the neighboring cell centers. If the
particles are unevenly distributed after advection, a generalized remeshing
method will modify both the position of particles and the value of passive
scalars. To reduce the expense, an adaptive strategy is introduced to decrease
the frequency of remeshing.

Figure 1: Flow chart for the hybrid method.

###  2.1 The Eulerian component

On the Eulerian side, we use a level set function $\phi(x,y,t)$ to represent
the interface [24]. The level set function has the property $|\nabla\phi|=1$
as signed distance, and identifies the inside region from the outside by the
signs: $\phi<0$ corresponds to the inside, $\phi>0$ the outside, and $\phi=0$
the interface. The level set function can be updated by solving the advection
equation $\frac{\partial\phi}{\partial t}+u\phi_{x}+v\phi_{y}=0$ with
reinitialization steps [25, 26]. The velocity field is either given with
analytical function for simple cases or obtained from the solution of two-
phase N-S equations [23].

###  2.2 The Lagrangian component

We use Lagrangian particles to solve the transport equation of the passive
scalar $\frac{\partial C}{\partial t}+\vec{u}\cdot\nabla_{s}C=0$, where
$\nabla_{s}$ is the gradient on the surface. The mass of every particle is
initialized as $m_{p}=C_{p}V_{p}$, where the concentration $C$ is assigned by
a given initial condition. The concentration field can be approximated by a
weighted average of the mass of the particles within the cut-off radius.

$C_{p}=\Sigma_{q}{W_{pq}m_{q}}.$ (1)

There are three main steps to solve the transport equation: particle
initialization, particle advection and remeshing. The first step is to
initialize the particles. We generate the particles at the centers of cut
cells, which can be identified by the signs of level set values at 9 points
($8$ points on cell boundary and one at cell center). Then, the particles are
pushed to the interface with a distance $|\phi|$ in the direction
$-sgn(\phi)\vec{n}$. However, this distribution of particles is non-uniform,
which will cause a zeroth-order inconsistency.

Figure 2: A sketch for the relaxation algorithm. The red and yellow particles
are the neighboring particles and the considered particle respectively. The
position is adjusted (purple arrow) by the resultant force of the repulsive
force (green dashed arrow) between every pair of particles, and then confine
the particles to the interface (black arrow).

To improve the consistency, a relaxation method proposed by Fu et al. [27] and
simplified by Zhu et.al [28] for SPH particle generation, is adopted here to
assure a uniform distribution of particles along the interface, see Fig. 2.
The second step is to advect the particles with the velocity interpolated from
the level set field by Eq. (2).

$\frac{d\vec{x_{p}}}{dt}=\vec{u}.$ (2)

The final step is to remesh particles for avoiding particle clustering during
advection. The remeshing process consists three sub-steps: resampling
(optional), relaxation, and redistribution.

The resampling enables adjusting the number of particles, which can be very
helpful to avoid void regions when the interface is extremely stretched.
However, after resampling a uniform distribution of particles often requires
many relaxation iterations. To reduce the frequency of remeshing, we introduce
a criterion based on a measure of the particle distribution,

$Q=\frac{\max\\{|v_{i}-\bar{v}|\\}}{\bar{v}}$ (3)

where $\bar{v}=\sum_{p}v_{p}/n$ is the average volume of particles and $n$ is
the number of particles, $v_{p}=\frac{1}{\Sigma_{q}{W_{pq}}}$ is the smoothed
volume [29] of the particle $p$, and $W_{pq}$ is the 5th-order Wendland kernel
function $W(|x_{p}-x_{q}|,h)$ with smoothed length $h$. Uniform distributions
are indicated by $Q\ll 1$. We can measure the distribution quality every $N$
time steps and in the following we set $N=30$. Accordingly, we use the
following remeshing criterion to control $Q$:

  * 1.

$Q<Q_{0}$, we remesh particles. In the following we use $Q_{0}=0.01$.

  * 2.

$Q_{1}<Q$, we resample particles followed by relaxation iterations until
$Q<Q_{0}$. Finally, we redistribute the mass from old particles to the new
ones. In the following we set $Q_{1}=0.03$.

  * 3.

$Q_{0}\leq Q\leq Q_{1}$, we do not resample particles but do $M$ times
relaxation iterations. To save costs, $M$ should be much smaller than the
number of iterations when $Q>Q_{1}$. In the following we use $M=30$. Finally,
we redistribute the mass from pre-existing particles to the resampled ones.

The redistribution process transfers the passive scalar from the old particles
to the new ones, see Fig. 3.

Figure 3: Redistribution of the mass from old cloud (green points) to new
cloud (red points). The yellow point and black curve is the particle and the
interface at current time step, while the green point and blue curve is the
particle and interface at next time step.

The redistribution method is generalized from the traditional particle
remeshing method [19]. The key idea is to use a new set of particles to
represent the old field of passive scalar represented by the old ones. We can
define the error of concentration between the two set of particles as

$\epsilon_{c}(x)=C^{old}(x)-C^{new}(x),$ (4)

where $C^{old}$ and $C^{new}$ are the concentrations represented by the old
and new particle sets, respectively. The redistribution process can be
expressed as

$m_{q}=\Sigma_{p}{m_{p}\beta_{pq}},$ (5)

where $\left\\{x_{p}\right\\}$ is the old set, $\left\\{x_{q}\right\\}$ is the
new one, and the subscript $m$ means the particle mass. In order to preserve
the accuracy, we can look for a proper set of coefficients $\beta_{pq}$ in Eq.
(5) to minimize the error defined in Eq. (4) to a given order $O(h^{p})$. Here
we let $p=4$ to balance the cost and accuracy.

The $C^{old}(x)$ and $C^{new}(x)$ in Eq. (4) can be further written as

$C^{old}=\Sigma_{p}{W(|x-x_{p}|,h)m_{p}},\quad
C^{new}=\Sigma_{q}{W(|x-x_{q}|,h)m_{q}}.$ (6)

Thus the error becomes

$\displaystyle\epsilon_{c}(x)$ $\displaystyle=C^{old}(x)-C^{new}(x)$ (7)
$\displaystyle=\Sigma_{p}{W(|x-x_{p}|,h)m_{p}}-\Sigma_{q}{W(|x-x_{q}|,h)m_{q}}$
(8)
$\displaystyle=\Sigma_{p}{W(|x-x_{p}|,h)m_{p}}-\Sigma_{q}{W(|x-x_{q}|,h)\Sigma_{p}{m_{p}\beta_{pq}}}$
(9)
$\displaystyle=\Sigma_{p}{m_{p}[W(|x-x_{p}|,h)-\Sigma_{q}{W(|x-x_{q}|,h)\beta_{pq}}]}.$
(10)

The function $W(|x-x_{q}|,h)$ can be approximated by a Taylor expansion at
$(x-x_{p})$ to third order as

$W(|x-x_{q}|,h)=W((|x-x_{p}|,h)+\Sigma_{k=1}^{3}{W_{y}^{(k)}|_{y=x_{p}}\frac{(x_{q}-x_{p})^{k}}{k!}}+{\cal
O}((x_{p}-x_{q})^{4}).$ (11)

We can minimize the error by eliminating the first 4 terms of
$(x_{q}-x_{p})^{k}$ with suitable weights $\beta_{pq}$ from

$\left(\begin{array}[]{cccc}1&1&1&1\\\
x_{p}-x_{1}&x_{p}-x_{2}&x_{p}-x_{3}&x_{p}-x_{4}\\\
(x_{p}-x_{1})^{2}&(x_{p}-x_{2})^{2}&(x_{p}-x_{3})^{2}&(x_{p}-x_{4})^{2}\\\
(x_{p}-x_{1})^{3}&(x_{p}-x_{2})^{3}&(x_{p}-x_{3})^{3}&(x_{p}-x_{4})^{3}\end{array}\right)\left(\begin{array}[]{c}\beta_{1}\\\
\beta_{2}\\\ \beta_{3}\\\ \beta_{4}\\\
\end{array}\right)=\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\
\end{array}\right)$ (12)

This coefficient matrix is a transformed Vandermonde matrix, such that the
solution exists and is unique:

$\displaystyle\beta_{0}$
$\displaystyle=\frac{x_{4}}{x_{4}-x_{1}}\cdot\frac{x_{3}}{x_{3}-x_{1}}\cdot\frac{x_{2}}{x_{2}-x_{1}}$
(13) $\displaystyle\beta_{1}$
$\displaystyle=\frac{x_{4}}{x_{4}-x_{2}}\cdot\frac{x_{3}}{x_{3}-x_{2}}\cdot\frac{x_{1}}{x_{1}-x_{2}}$
(14) $\displaystyle\beta_{2}$
$\displaystyle=\frac{x_{4}}{x_{4}-x_{3}}\cdot\frac{x_{2}}{x_{2}-x_{3}}\cdot\frac{x_{1}}{x_{1}-x_{3}}$
(15) $\displaystyle\beta_{3}$
$\displaystyle=\frac{x_{3}}{x_{3}-x_{4}}\cdot\frac{x_{2}}{x_{2}-x_{4}}\cdot\frac{x_{1}}{x_{1}-x_{4}}$
(16)

Apparently, this method can be generalized to arbitrary order if there are
enough neighboring particles in the compact support of each particle.

###  2.3 Extension to multi-region problems

The multi-region problem can be solved with a multiple level-set method [30],
and combined with the present method in a straightforward way. However, with
increasing number of regions, the multiple level-set method demands larger
memory size because each region is represented by a separated level-set field.
This problem is avoided by a regional level-set method [21, 22], in which only
a pair of data is stored in the memory, i.e. an unsigned distance function and
a region index (also called as color index). However, to implement the
regional level-set method, some modifications are necessary. Before
introducing the modifications, we review the mechanism and some related basic
concepts of the regional level-set method.

Similar to classical level-set method, each region is still tracked by its own
level-set field, but this level-set field, called local level-set field, is
constructed from the regional level-set field. After advancing the local
level-set field, the regional level-set field can be updated by the new local
level-set field. With assignment of all indicator colors, the updating of the
regional level-set field is done. This process can be abstracted as a
reconstruction mapping, which maps all local level-set fields to a regional
level-set field.

We denote the local level-set function of color index $r$ as $\phi^{r}(x,y)$,
with positive value inside and negative outside, and the regional level-set
function as $\psi=(\phi(x,y),r(x,y))$, where the function $\phi$ here is the
unsigned distance function from the cell center to the nearest interface, and
$r$ is its indicator color. Because the local level-set fields will not be
stored in the memory, it needs to be constructed by the construction mapping
based on the regional level-set field. The construction mapping from regional
level-set field to local level-set field with color $\chi$ can be represented
as

$\phi^{\chi}=\mathcal{C}^{\chi}(\phi,r)=\bigg{\\{}\begin{aligned} \phi,\quad
r=\chi\\\ -\phi,\quad r\neq\chi\end{aligned},$ (17)

where $\mathcal{C}^{\chi}$ is the construction mapping from regional level-set
field to local level-set field indexed $\chi$. If $\chi$ is removed,
$\mathcal{C}$ is the reconstruction mapping from regional level-set field to
all local level sets.

Because only the regional level-set field is stored in the memory, the
regional level-set field needs to be updated from all local level-set fields.
For cell $C$ whose center is $(x,y)$, if there exists at least one positive
local level-set field, the reconstruction mapping can be defined as

$\psi=(\phi(x,y),\chi(x,y))=\mathcal{R}(\\{\phi^{r},r\in I^{+}\\})=(\min_{r\in
I^{+}}\\{\phi^{r}\\},\operatorname*{arg\,min}_{r\in I^{+}}(\phi^{r})),$ (18)

where $I^{+}$ is the set of local indices with positive local level-set value.
If all the local level-set values are negative, we need to choose the one with
smallest absolute value, and the reconstruction mapping becomes

$\psi(\phi(x,y),\chi(x,y))=\mathcal{R}(\\{\phi^{r},r\in I\\})=(\min_{r\in
I}\\{|\phi^{r}|\\},\operatorname*{arg\,min}_{r\in I}(\phi^{r})),$ (19)

where $I$ is the set of all local indices. The constructing operator is not
one-to-one, which causes an inconsistency
$\mathcal{RC}((\phi,\chi))\neq(\phi,\chi)$, see Fig. 7. This will introduce
extra error and lead to an inaccurate concentration field of the passive
scalar. Thus, we introduce some extra treatments to reduce the error. Here, we
consider the application of present method with the $2D$ regional level-set
method.

####  2.3.1 Correction near triple points

Figure 4: An example of the error introduced by regional level-set method near
the triple points. For the level-set field of region $1$, the green line is
the exact distance while the red line is the inaccurate distance reproduced by
the regional level set method.

We introduce two techniques to fix this problem: a particle-based correction
for the local level-set field and a local mass exchange to smoothing the
passive scalar. Since particles on the interface are readily available, we can
use them to improve the local level-set field without requiring additional
memory. Specifically, we utilize particles to correct the distance in $3\times
3$ set of neighboring cells around the triple point. There are two important
parts for correcting the local level-set field , the correction of the sign
and of the distance.

Figure 5: Measuring the distance. Looking for two nearest particles to the
cell center $C$, denote them as $P_{1}$ and $P_{2}$. Calculate the height of
$\triangle CP_{1}P_{2}$ to the edge $P_{1}P_{2}$ as the distance.

The simplest method to estimate the distance is looking for the 2 nearest
particles of a cell, see Fig. 5.

Figure 6:  Determination of the sign of the level-set function near a triple
point. Yellow and blue arrows are the normal vectors estimated by the
disturbed interface (yellow curve) and particles (red points) respectively.

On the other hand, a wrong sign is more critical because it may introduce a
very large error, which can be about two times as the distance. However, this
can be modified by using particles as well. Firstly, we estimate the
interfacial normal vector from particles, see Fig. 6. Its direction points
from the cell center to the interface segment. We also have a normal vector
estimated by the the local level-set field
$\vec{n}=\frac{\nabla\phi}{|\nabla\phi|}$. If the angle between these two
normal vectors is larger than 90 degrees, the sign of local level-set is wrong
and should be changed.

The second technique is to introduce mass exchange near the triple point to
reduce the oscillations of concentration field. Because of the limited
resolution provided by the particles, the error introduced by the regional
level-set method can not be resolved completely by the above correction. In a
$3\Delta x\times 3\Delta y$ neighboring region of the triple point, we
introduce a mass exchange between each pair of distinct particles $p$ and $q$

$f_{pq}=-(M_{p}-M_{q})W_{pq}\frac{V_{p}+V_{q}}{2},$ (20)

where $f_{pq}$ is the net amount of mass exchanged between particle $p$ and
particle $q$. $M_{p}$ and $M_{q}$ are the mass of the two particles
respectively, $W_{pq}$ is the weight function of the distance between $p$ and
$q$, and $V_{p}$ and $V_{q}$ is the volume of the particles respectively.

The corrected mass of particle $p$ is $M_{p}^{*}$

$\displaystyle M_{p}^{*}$ $\displaystyle=M_{p}-\sum_{q\neq
p}{(M_{p}-M_{q})W_{pq}\frac{V_{p}+V_{q}}{2}}$ (21)
$\displaystyle=M_{p}-\sum_{q}{(M_{p}-M_{q})W_{pq}\frac{V_{p}+V_{q}}{2}}$ (22)
$\displaystyle\approx\sum_{q}{M_{q}W_{pq}V},$ (23)

in which the last approximation results from assuming that the volume of the
particles is nearly equal to the average volume $V$. This mass exchange makes
the mass of every particle in this region be a weighted average of all the
other particles in this region. Obviously, this process maintains the mass
conservation. Moreover, the measure of influenced neighboring cells (namely
the total area) tends to zero as the resolution tends to infinite, thus this
mass exchange is consistent.

##  3 Test cases

Different test cases are considered in this section. The simplest one is
translational movement. The uniform case tests the error from the remeshing
process, and the non-uniform case measures the numerical diffusion.
Furthermore, a droplet in the given simple shear flow without surface tension
and in two-phase shear flow with surface tension is considered. Besides the
multi-region case tests the compatibility with regional level set method.

###  3.1 Translation of passive scalar distributed on a circle

The transport of passive scalar with uniform and non-uniform initial
distribution on a unit circle in a uniform and constant flow field is tested.
Theoretically, the concentration of the passive scalar on the interface should
remain unchanged.

The computational domain is $[0,1]\times[0,1]$ on the $(x,z)$ plane. The
circle with radius $r=0.2$ is centered at $(0.3,0.3)$ in a constant flow
$u=0.4/\sqrt{2}$ and $w=0.4/\sqrt{2}$. The final time is $t_{final}=1.0$. The
initial concentrations are $1.0$ and $c(\theta)=2+\cos(\theta)$, where
$\theta$ is the angle, for the uniform and non-uniform cases, respectively.

To study the accuracy of remeshing and control for the impact of particle
distribution, we conduct particle relaxation at each time step with a
sufficiently large number of iterations to achieve a uniform particle
distribution. However, in practical use, remeshing is typically only required
when particle distribution becomes degenerate, and therefore a much smaller
number of iterations is usually sufficient.

In the first part, the uniform case is tested to verify the zeroth-order
consistency of the redistribution method. Fig. 7 and Tab. 1 show that our
redistribution method has good zeroth-order consistency and convergence
properties when a sufficiently uniform particle distribution is achieved.

(a) Initial Condition

(b) $t=1.0$

Figure 7: Translation case with a uniform initial concentration, $256\times
256$ cells

Nnumerical diffusion is tested by a non-uniform case. From the results in Fig.
8 and Tab. 2, we can see that the maximum and minimum concentration is
preserved in the advection process, which implies the numerical diffusion is
quite small.

(a) Initial Condition

(b) $t=1.0$

Figure 8: Translation case with non-uniform initial concentration,$256\times 256$ cells. Table 1: The relative error of translation case with uniform concentration Resolution | $64\times 64$ | $96\times 96$ | $128\times 128$ | $256\times 256$
---|---|---|---|---
Relative Error | 1.932% | 1.005% |  0.879% |  0.994%
Table 2: The range of concentration in non-uniform case Resolution |  $64\times 64$ |  $128\times 128$ |  $256\times 256$
---|---|---|---
Range | [0.999, 2.992] | [1.0018, 3.00174] | [0.999, 3.002]

###  3.2 Droplet in the simple shear flow

A circle in the simple shear flow can have a large deformation [5], if there
is no surface tension. We compare the results with the analytical solution to
verify the accuracy of our method.

####  3.2.1 Derivation of the analytical solution

The governing equation is the mass conservation:

$\frac{d}{dt}\int_{L(t)}{C(l,t)dl}=0,$ (24)

where $L(t)$ is an arbitrary element on the interface at time $t$, and $l$ is
the arc length. In other words, the total mass $M_{L}$ is constant on an
element

$\int_{L(t)}{C(l,t)dl}=M_{L}\bigg{|}_{(t=0)}.$ (25)

The right hand side can be obtained by the initial condition. The initial
condition is that the interface is a unit circle centered at origin, and the
passive scalar has a uniform distribution $C=1$, thus
$M_{L}=\theta_{1}-\theta_{0}$, where $\theta_{0}$ and $\theta_{1}$ are the two
ends of the considered curve $L(0)$. $L(t)$ is continuously deformed by the
given velocity field, thus it can be determined by velocity field and initial
element $L(0)$. We assume that the flow is: $u=0$, $v=Mx$. L(0) can be
parameterized as

$\displaystyle x(\theta)$ $\displaystyle=\cos\theta,$ (26) $\displaystyle
y(\theta)$ $\displaystyle=\sin\theta,$ (27)

where $\theta\in[\theta_{0},\theta_{1}]$. Points on $L(t)$ are given as

$\displaystyle x(\theta,t)$ $\displaystyle=\cos\theta,$ (28) $\displaystyle
y(\theta,t)$ $\displaystyle=\sin\theta+Mxt=\sin\theta+Mt\cos\theta,$ (29)

and $dl$ can be represented by

$dl=\left|\frac{\partial}{\partial\theta}\vec{r}(\theta,t)\right|d\theta,$
(30)

where $\vec{r}(\theta,t)=(x(\theta,t),y(\theta,t))$. Substituting this into
Eq. (25), we have

$\int_{\theta_{0}}^{\theta_{1}}{C(\theta,t)\left|\frac{\partial}{\partial\theta}\vec{r}(\theta,t)\right|d\theta}=\theta_{1}-\theta_{0},$
(31)

which holds for arbitrary $\theta_{0}$ and $\theta_{1}$. Consider the partial
derivative with respect to $\theta_{1}$, we have

$C(\theta_{1},t)\left|\frac{\partial}{\partial\theta}\vec{r}(\theta_{1},t)\right|=1.$
(32)

The second term in the left hand side can be obtained from Eq. (28-29). For
this simple shear flow case, we have

$\left|\frac{\partial}{\partial\theta}\vec{r}(\theta,t)\right|=\sqrt{1+M^{2}t^{2}\sin^{2}\theta-2Mt\sin\theta\cos\theta}.$
(33)

Therefore, the solution for arbitrary $\theta$ is

$C(\theta,t)=\frac{1}{\sqrt{1+M^{2}t^{2}\sin^{2}\theta-2Mt\sin\theta\cos\theta}}.$
(34)

Here we use $\theta$ to replace $\theta_{1}$ because $\theta_{1}$ can be
arbitrary in Eq. (32).

####  3.2.2 Numerical Results

(a) Initial Condition

(b) $M=0.5$, time$=1.0$

Figure 9: A circle in simple shear flow

In this test case, $M=0.5$, and the final time $t=1.0$. The particles are
remeshed after every step, with a target distribution quality
$Q_{0}=Q_{1}=1\%$. The maximum value of the passive scalar exact solution is
about $1.2806$. It can be shown from the Fig. 9 that the maximum value of the
numerical result is about $1.28$. The $L_{\infty}$ norm of relative error of
concentration field is shown in Tab. 3 .

Table 3: The relative error of concentration field on droplet in the given shear flow. Resolution | $128\times 128$ | $256\times 256$ | $512\times 512$
---|---|---|---
Relative Error | 0.504% |  0.397% |  0.289%

We emphasize that, in addition to remeshing, errors in the simulation results
may be introduced by other factors, such as the particle distribution and the
consistency of the SPH method. Although sufficient iterations of relaxation
have been done, the particle distribution is not perfectly uniform. Even with
a perfectly uniform one, the SPH method still can not achieve the zero-th
order exactly. Thus, these errors cannot be remedied solely by increasing the
resolution of the simulation.

###  3.3 Drop deformation Case

We consider the deformation of a droplet immersed in another shearing fluid
[31] to verify our method in two-phase flow. There are two walls at top and
bottom of the bulk fluid moving in opposite directions with constant velocity
$\pm u_{\infty}$, which induces a shear flow in the bulk phase, and finally
leads to the deformation of the immersed droplet. The deformation was measured
experimentally [32], theoretically [31] and by numerical simulation [11, 33,
29]. The geometry is $8R_{0}\times 8R_{0}$ in $XoZ$ plane, and the center of
the drop is $(4R_{0},4R_{0})$. The ratio of the viscosity
$\lambda=\eta_{d}/\eta_{w}=100$. The shear rate is $G=2u_{\infty}/8R_{0}=1$.
The non-dimensionless parameters are

$\displaystyle Ca=\frac{G\eta_{w}R_{0}}{\alpha}=0.15,$ (35) $\displaystyle
Re=\frac{\rho GR_{0}^{2}}{\eta_{w}}=1.$ (36)

We use the fluid solver in [33], which has been verified before. The different
thing is that we consider the variations of the concentration field on the
interface with an uniform initial distribution $C=1$.

In the past study, the convergence of the deformation is studied, and the
shape almost converges even in a low resolution, like $128\times 128$. If only
the interface advection is considered, the surface concentration is only
dependent on the interface velocity. We sample the interfacial velocity on the
particles and find that the convergence for the velocity requires much higher
resolution than the shape, see Fig. 11.

Figure 10: The convergence study of interfacial velocity magnitude Figure 11:
The convergence study of interfacial velocity components

Finally we observe from Fig. 10 that the velocity almost converges above the
$1024\times 1024$ resolution, but usually $128\times 128$ or $256\times 256$
(or equivalent number of particles in particle methods) is enough to track the
interface deformation [23, 29]. It is worth stressing that the convergence of
velocity is just the prerequisite for the convergence of the concentration for
pure advection. We can see that the magnitude of the interfacial velocity is
highly non-uniform, which results in some points being highly compressed while
others are stretched. This non-uniformity makes convergence more difficult. If
there is no diffusion, the concentration near the compressed region would
become infinitely high, which is unphysical and violates the mono-layer
assumption. Therefore, a mechanism is needed to constrain the concentration
from approaching infinity. Therefore, we introduce a diffusion term into the
passive scalar transport and the transport equation becomes

$\frac{\partial C}{\partial t}+C\nabla_{s}\cdot\vec{u}=D\nabla^{2}_{s}C.$ (37)

The discretized scheme is

$\frac{Dm}{Dt}\bigg{|}_{i}=2D\sum_{k\neq
i}{\frac{(C_{i}-C_{k})}{r_{ik}}}W^{\prime}_{ik}V_{i}V_{k},$ (38)

where the $D=1.0$ is the diffusion coefficient and the corresponding Peclet
number is $Pe=GR_{0}^{2}/D=1$. Eq. (38) is based on [34]. We modify it from
the concentration to the mass by multiplying the volume of the particle.

After introducing the diffusion term, the convergence of the concentration
field is improved since the diffusion term helps to smooth out large
concentrations caused by high compression, see Fig. 12. Moreover, the
deformation can be measured by a deformation parameter which is defined as

$D=Ca\frac{16+19\lambda}{16+16\lambda}.$ (39)

The theoretical and numerical deformation parameters for a clean droplet are
0.1778 and 0.1718 (with resolution $512\times 512$). From Fig. 12, we can see
that the deformation parameter is almost convergent at resolution of
$1024\times 1024$, and the parameter is much larger than the clean one, which
implies the Marangoni force plays an important role.

(a) Initial Condition, $1024\times 1024$

(b) $T=8,\ 256\times 256$, $D=0.1886$

(c) $T=8,\ 512\times 512$, $D=0.2069$

(d) $T=8,\ 1024\times 1024$, $D=0.2102$

Figure 12: A droplet in shear flow with different resolution, $Q_{0}=1\%$ and
$Q_{1}=3\%$.

We would like to stress that relaxation iterations are usually very expensive.
However, with our method, relaxation is only invoked when the particle
distribution is not sufficiently uniform. Moreover, with the adaptive
strategy, overheads are drastically reduced. We measured the costs (wall clock
time) in a $128\times 128$ droplet deformation case on a workstation with an
I9-10980XE 18-core CPU. The relaxation process only costs a total of $0.959$
seconds, which is a very small fraction compared to the full cost of $2362.11$
seconds.

###  3.4 Multi-region case: Two tracks case

This case is constructed by the following steps. First, on the $(x,z)$ plane,
we draw two circles with $r=0.1$ close to each other, centered at $x_{1}=0.3$,
$x_{2}=0.5$. We slice the two circles at $z=0.5$. Then we insert straight
lines in between the $4$ half circles with length 0.2, see Fig. (13(a)).

The computational domain is $1.0\times 1.0$. The background velocity is given
as $u=w=0.3/\sqrt{2}$. The resolution in this case is $128\times 128$. For the
uniform case, the initial condition is $C=1$. For the non-uniform case, the
initial concentration depends linearly on the height from $2$ to $3$. Such a
multi-region case will not produce the singularity, namely the interface near
the triple point is $C^{1}$ smooth.

(a) Uniform initial condition

(b) Uniform case, time$=1.0$

(c) Non-uniform initial condition

(d) non-uniform, time$=1.0$

Figure 13: Translation of two tracks

From Fig. 13, we can observe that the relative error of concentration near the
triple point is approximately $3.5\%$ in the uniform case. This error is
relatively higher compared to other cases due to the less accurate regional
level set near the triple point, despite our efforts to enhance its accuracy
through various techniques in Sec. 2.3.1. However, these techniques have
significantly improved the accuracy, resulting in a four to five times
improvement as illustrated in Fig. 14. If the local level-set fields near
triple points can be further improved, both the accuracy and convergence
reported in Section 3.1 and Section 3.2 can be achieved.

(a) Interface for left region, time$=1.0$

(b) Concentration distribution, about $16\%$ error, time$=1.0$

Figure 14: The results for the left region, without modifications for uniform
initial concentration.

##  4 Conclusion

To simulate the concentration evolution of a passive scalar on a moving
interface, we propose a hybrid Eulerian-Lagrangian method inherits the
advantages of either side. Conservation of mass is preserved by the particle
method, and the level-set method captures the evolution of the interface. For
validation, we use a circle with a translational motion to prove the accuracy.
The error introduced by SPH method and generalized remeshing method is
reasonably small and can be smaller as resolution increases. Subsequently, we
test a droplet in a given shear flow with equally good results. Furthermore,
the droplet in the two-phase shear flow is validated with converging results.
Finally, the two-track case is devised and tested for multi-region problems.
Although the relaxation algorithm is expensive, the overhead is well
controlled by using interfacial particles with co-dimension one and reducing
the frequency of remeshing adaptively. The Lagrangian interpolation, due to
its unavailability in higher dimensions, presents a challenge in extending it
to three-dimensional space, and requires further investigation in the future.
Furthermore, the enhancement of regional level-set method accuracy near triple
points remains an unresolved issue that warrants scholarly attention.

##  5 References

## References

  * [1] P. Chu, J. Finch, G. Bournival, S. Ata, C. Hamlett, R. J. Pugh, A review of bubble break-up, Advances in Colloid and Interface Science 270 (2019) 108–122.
  * [2] R. Chowdhury, N. P. Holmes, N. Cooling, W. J. Belcher, P. C. Dastoor, X. Zhou, Surfactant engineering and its role in determining the performance of nanoparticulate organic photovoltaic devices, ACS omega 7 (11) (2022) 9212–9220.
  * [3] S. Chowdhury, S. Shrivastava, A. Kakati, J. S. Sangwai, Comprehensive review on the role of surfactants in the chemical enhanced oil recovery process, Industrial & Engineering Chemistry Research 61 (1) (2022) 21–64.
  * [4] D. Adalsteinsson, J. A. Sethian, Transport and diffusion of material quantities on propagating interfaces via level set methods, Journal of Computational Physics 185 (1) (2003) 271–288.
  * [5] J.-J. Xu, H.-K. Zhao, An eulerian formulation for solving partial differential equations along a moving interface, Journal of Scientific Computing 19 (1) (2003) 573–594.
  * [6] J.-J. Xu, Z. Li, J. Lowengrub, H. Zhao, A level-set method for interfacial flows with surfactant, Journal of Computational Physics 212 (2) (2006) 590–616.
  * [7] F. S. Schranner, N. A. Adams, A conservative interface-interaction model with insoluble surfactant, Journal of Computational Physics 327 (2016) 653–677.
  * [8] K. E. Teigen, P. Song, J. Lowengrub, A. Voigt, A diffuse-interface method for two-phase flows with soluble surfactants, Journal of computational physics 230 (2) (2011) 375–393.
  * [9] M. A. Olshanskii, A. Reusken, X. Xu, A stabilized finite element method for advection–diffusion equations on surfaces, IMA journal of numerical analysis 34 (2) (2014) 732–758.
  * [10] A. J. James, J. Lowengrub, A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant, Journal of computational physics 201 (2) (2004) 685–722.
  * [11] S. Adami, X. Hu, N. A. Adams, A conservative sph method for surfactant dynamics, Journal of Computational Physics 229 (5) (2010) 1909–1926.
  * [12] M. Wang, Y. Deng, X. Kong, A. H. Prasad, S. Xiong, B. Zhu, Thin-film smoothed particle hydrodynamics fluid, ACM Transactions on Graphics (TOG) 40 (4) (2021) 1–16.
  * [13] H. D. Ceniceros, The effects of surfactants on the formation and evolution of capillary waves, Physics of Fluids 15 (1) (2003) 245–256.
  * [14] M.-C. Lai, Y.-H. Tseng, H. Huang, An immersed boundary method for interfacial flows with insoluble surfactant, Journal of Computational Physics 227 (15) (2008) 7279–7293.
  * [15] M.-C. Lai, Y.-H. Tseng, H. Huang, et al., Numerical simulation of moving contact lines with surfactant by immersed boundary method, Communications in Computational Physics 8 (4) (2010) 735.
  * [16] W. C. de Jesus, A. M. Roma, M. R. Pivello, M. M. Villar, A. da Silveira-Neto, A 3d front-tracking approach for simulation of a two-phase fluid with insoluble surfactant, journal of Computational Physics 281 (2015) 403–420.
  * [17] T. Y. Hou, J. S. Lowengrub, M. J. Shelley, Removing the stiffness from interfacial flows with surface tension, Journal of Computational Physics 114 (2) (1994) 312–338.
  * [18] M. Botsch, L. Kobbelt, A remeshing approach to multiresolution modeling, in: Proceedings of the 2004 Eurographics/ACM SIGGRAPH symposium on Geometry processing, 2004, pp. 185–192.
  * [19] G.-H. Cottet, P. D. Koumoutsakos, et al., Vortex methods: theory and practice, Vol. 8, Cambridge university press Cambridge, 2000.
  * [20] D. Rossinelli, P. Koumoutsakos, Vortex methods for incompressible flow simulations on the gpu, The visual computer 24 (2008) 699–708.
  * [21] W. Zheng, J.-H. Yong, J.-C. Paul, Simulation of bubbles, Graphical Models 71 (6) (2009) 229–239.
  * [22] S. Pan, X. Y. Hu, N. A. Adams, High-resolution method for evolving complex interface networks, Computer Physics Communications 225 (2018) 10–27.
  * [23] X. Y. Hu, B. Khoo, N. A. Adams, F. Huang, A conservative interface method for compressible flows, Journal of Computational Physics 219 (2) (2006) 553–578.
  * [24] S. Osher, R. P. Fedkiw, Level set methods: an overview and some recent results, Journal of Computational physics 169 (2) (2001) 463–502.
  * [25] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, Journal of Computational physics 114 (1) (1994) 146–159.
  * [26] L. Fu, X. Y. Hu, N. A. Adams, Single-step reinitialization and extending algorithms for level-set based multi-phase flow simulations, Computer Physics Communications 221 (2017) 63–80.
  * [27] L. Fu, L. Han, X. Y. Hu, N. A. Adams, An isotropic unstructured mesh generation method based on a fluid relaxation analogy, Computer Methods in Applied Mechanics and Engineering 350 (2019) 396–431.
  * [28] Y. Zhu, C. Zhang, Y. Yu, X. Hu, A cad-compatible body-fitted particle generator for arbitrarily complex geometry and its application to wave-structure interaction, Journal of Hydrodynamics 33 (2) (2021) 195–206.
  * [29] X. Hu, N. A. Adams, An incompressible multi-phase sph method, Journal of computational physics 227 (1) (2007) 264–278.
  * [30] B. Merriman, J. K. Bence, S. J. Osher, Motion of multiple junctions: A level set approach, Journal of computational physics 112 (2) (1994) 334–363.
  * [31] G. I. Taylor, The viscosity of a fluid containing small drops of another fluid, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 138 (834) (1932) 41–48.
  * [32] G. I. Taylor, The formation of emulsions in definable fields of flow, Proceedings of the Royal Society of London. Series A, containing papers of a mathematical and physical character 146 (858) (1934) 501–523.
  * [33] J. Luo, X. Hu, N. A. Adams, A conservative sharp interface method for incompressible multiphase flows, Journal of Computational Physics 284 (2015) 547–565.
  * [34] L. Brookshaw, A method of calculating radiative heat diffusion in particle simulations, Publications of the Astronomical Society of Australia 6 (2) (1985) 207–210.