# A finite-difference ghost-point multigrid method for multi-scale modelling
of sorption kinetics of a surfactant past an oscillating bubble

Clarissa Astuto  King Abdullah University of Science and Technology (KAUST),
4700, Thuwal, Saudi Arabia  Armando Coco  Department of Mathematics and
Computer Science, University of Catania, Viale Andrea Doria 6, 95125, Catania,
Italy  Giovanni Russo  Department of Mathematics and Computer Science,
University of Catania, Viale Andrea Doria 6, 95125, Catania, Italy

###### Abstract

We propose a method for the numerical solution of a multiscale model
describing sorption kinetics of a surfactant around an oscillating bubble. The
evolution of the particles is governed by a convection-diffusion equation for
the surfactant concentration $c$, with suitable boundary condition on the
bubble surface, which models the action of the short range attractive-
repulsive potential acting on them when they get sufficiently close to the
surface [1]. In the domain occupied by the fluid, the particles are
transported by the fluid motion generated by the bubble oscillations.

The method adopted to solve the equation for $c$ is based on a finite-
difference scheme on a uniform Cartesian grid and implemented in 2D and 3D
axisymmetric domains. We use a level-set function to define the region
occupied by the bubble, while the boundary conditions are discretized by a
ghost-point technique to guarantee second order accuracy at the curved
boundary. The sparse linear system is finally solved with a geometric
multigrid technique designed ad-hoc for this specific problem. Several
accuracy tests are provided to prove second order accuracy in space and time.

The fluid dynamics generated by the oscillating bubble is governed by the
Stokes equation solved with a second order accurate method based on a
monolithic approach, where the momentum and continuity equations are solved
simultaneously. Since the amplitude of the bubble oscillations are very small,
a simplified model is presented where the computational bubble is actually
steady and its oscillations are represented purely with time-dependent
boundary conditions. A numerical comparison with the moving domain model
confirms that this simplification is perfectly reasonable for the class of
problems investigated in this paper.

##  1 Introduction

Diffusion equations and Stokes problems in a time dependent domain are present
in countless research areas with many applications such as problems of
temperature distribution [2], studies of biological pattern formation and cell
motility on evolving surfaces [3, 4, 5]. The most relevant topic that moves
our research is the modelling of surfactants in two-phase flows using a
diffuse interface, that is part of a long time project [6, 7, 8, 9, 10, 11,
12, 13, 14]. In these works the authors report the experimental trapping
kinetics of a diffusing flux of surfactants sticking at the surface of an
oscillating gas bubble set in the middle of the vessel (see Fig. 1 (a)). The
surfactant concentration past the oscillating bubbles is detected by
conductivity measurements (see Fig. 1 (b)). A different and unexpected
behavior is observed in presence of an empty bubble oscillating at resonance
frequency (black curve in Fig. 1 (b)). The phenomenon is particularly relevant
when the bubbles are exposed to intense forced oscillations near resonance.

Surfactants are important for several industrial applications, such as
processes of emulsification and mixing [15] or the production and
stabilization of 2D nanomaterials [16]. They can be soluble in at least one of
the fluid phases and the exchange of surfactants between the bulk phases and
the fluid interfaces is governed by the process of adsorption and desorption.
In [5] different phase field models are derived for two-phase flow with a
surfactant soluble in possibly both fluids. In [17, 18] they present models
that account for both the diffusive transport inside the solution and the
kinetics taking place at the interface using a free-energy formulation.

Figure 1: Experimental domain and results from [7, 6]. (a) Schematic setup of
the real apparatus. The central sphere represents the oscillating bubble. See
[7, 6] for a detailed description and for experimental values of $H$, $A$,
$L$. (b) Time evolution of the aqueous solution conductance measured above the
bubble (electrodes 2 and 3 of (a)). Red line: no bubble; blue line: saturated
bubble submitted to a flux of surfactants; black line: oscillating bubble
submitted to a flux of surfactants (adapted from [7]). The line thicknesses
represent an estimate of the experimental uncertainty of the conductivity
measurements.

A theoretical work of Ward and Tordai [19] formulated a time-dependent
relation between the surface density of surfactants adsorbed at an interface
and their concentration at the sub-surface layer of the solution, assuming a
diffusive transport from the bulk solution. Consequent theoretical works have
focused on providing a second closure relation between these two variables, as
in [1], while in this paper we focus on the numerical aspects of the problem.

Specifically, we study a diffusion equation in a bulk domain, with a dynamic
time-dependent boundary condition derived by conservation arguments and
stating that at the boundary the flux is proportional to the time derivative
of the solution [4, 20]. We formulate a finite-difference scheme for
advection-diffusion equation with moving curved surfaces/boundaries. Time-
discretization is performed with the Crank-Nicolson method. The bubble region
is implicitly described by a level-set function, while the implementation of
boundary conditions on complex-shaped boundaries/surfaces is based on a ghost-
point method.

The ghost-point approach for domains described by level-set functions has been
successfully proposed in several contexts [21, 22, 23, 24, 25, 26], as it has
the advantage to allow an implicit representation of the boundary and then it
can be employed on meshes that do not necessarily conform with a complex-
shaped boundary. This advantage alleviates the computational burden that is
associated with mesh generation steps, as observed in other approaches based
on fitted-boundary methods. A flexible ghost-point technique, suitable for
different boundary condition types, is presented in [27, 28] and applied to
several contexts [29, 30]. In this paper, we extend the approach to
accommodate time dependant boundary conditions and the presence of second
order tangential derivatives.

A geometric multigrid method is employed to efficiently solve the sparse
linear system arising from the discretization of the problem. The multigrid
approach is extended from [27] in order to account for time-dependent boundary
conditions and the presence of tangential derivatives. A suitable technique is
presented in order to maintain the optimal efficiency of the multigrid method
and contain the boundary effect degradation of the performance. The method is
second-order accurate in space and time, as confirmed by numerical tests.

The convection of particles is driven by the fluid motion around the
oscillating bubble, usually modelled by incompressible Navier-Stokes
equations. We assume the solute does not significantly change the density and
rheology of the fluids, therefore the fluid motion is independent of the
solute concentration (one way coupling). We also assume that the motion of the
bubble is assigned a priori and does not depend on the fluid dynamics. In a
more realistic scenario, the bubble surface deformation would be influenced by
the fluid motion, resulting in a two-way coupling as in fluid/membrane
interaction problems [31, 32].

For the physical parameters range adopted in the experiments, the Reynolds
numbers are very small so that the convective terms of the Navier-Stokes
equations can be neglected, and we can safely model the fluid-dynamics by the
Stokes equations instead. To numerically solve Stokes equations on moving
domains we employ the method proposed in [29] based on a monolithic approach.

The paper is structured as follows. In Sect. 2 we present the mathematical
model for the diffusion of particles. Sect. 3 describes the finite-difference
ghost-point technique to implement the time-dependent boundary conditions. In
Sect. 4 we present the multigrid method to solve the sparse linear system
arising from the ghost-point discretization. In Sect. 5 we perform several
numerical tests to prove the second order accuracy in 2D and 3D axisymmetric
geometries. Sect. 6 deals with the fluid dynamics generated by the bubble
oscillations. Small amplitude of the oscillations suggests that the bubble
motion can be modelled purely from time-dependent boundary condition for the
fluid velocity, while the computational bubble domain remains steady, saving
then meaningful computational efforts. This simplification is justified by
numerical tests. Finally, we couple the Stokes problem with the convection-
diffusion equations to model the particle concentration evolution around an
oscillating attracting bubble, proposing two types of oscillations (harmonic
and ellipsoidal). Conclusions are drawn in Sect. 7.

##  2 Multiscale Model

Modelling the diffusion in presence of a trap is challenging if multiple
scales are involved. In recent papers, such as [6, 7], the range of the
attractive-repulsive core of the trap is of the order of nanometers, a length
that is several orders of magnitude smaller than the size of the domain. In
order to overcome such difficulty, a multiscale model was proposed in [1],
which we briefly recall here, in the single carrier approximation.

The time evolution of a local concentration of ions $c=c(\vec{x},t)$ diffusing
in a steady fluid is governed by the conservation law

$\frac{\partial c}{\partial t}=-\nabla\cdot J.$ (1)

The flux term $J$ contains a diffusion and a drift term:

$J=\ -D\left(\nabla c+\ \frac{1}{k_{B}T}c\nabla V\right)$ (2)

where $k_{B}$ is the Boltzmann’s constant, $T$ is the absolute temperature and
$V=V(\vec{x})$ is a suitable attractive-repulsive potential that models the
particle trap. For simplicity, we describe the 1D model derivation only,
referring the reader to [1] for a detailed derivation of higher order
dimension models. In 1D, equations (1) and (2) read:

$\displaystyle\displaystyle\frac{\partial c}{\partial t}+\frac{\partial
J}{\partial x}$ $\displaystyle=0$ (3) $\displaystyle\displaystyle J$
$\displaystyle=-D\left(\frac{\partial c}{\partial
x}+\frac{1}{k_{B}T}\,c\,V^{\prime}\right)$ (4)

Assume that the trap is located at $x=0$. Initially, particles are located
within the region $x>0$. If they get close to the trap, they are attracted
towards $x=0$. This phenomenon is simulated by a potential $V(x)$ such that
$V^{\prime}(x)>0$ for $x>0$. On the other hand, if particles pass to the
region $x<0$, they are repulsed towards $x=0$. Then, $V^{\prime}(x)<0$ for
$x<0$. The attractive/repulsive mechanism is therefore modeled in the
neighborhood of $x=0$ by a short range potential $V(x)$ that is different from
zero only in a thin region around the trap, say
$\Omega^{\varepsilon}_{b}=[-\varepsilon,\varepsilon L]$ with $\varepsilon\ll
1$ and $L>0$. Therefore, $V(x)=0$ for $x\geq\varepsilon L$ and (assuming that
the first derivative of $V(x)$ is continuous) $V^{\prime}(\varepsilon L)=0$. A
typical shape of the potential is reported in the left panel of Fig. 2.

Figure 2: Representation of $V(x)$, on the left, and $U(\xi)$, on the right,
for $\varepsilon=0.05$ and $L=2$. On the left the dashed line $x=\varepsilon
L$ denotes the right boundary of $\Omega_{b}^{\varepsilon}$.

Assuming that there is a wall at $x=1$, the fluid domain is represented by
$\Omega_{f}^{\varepsilon}=[\varepsilon L,1]$. The problem consists of solving
(1) and (2) in
$\Omega^{\varepsilon}=\Omega_{b}^{\varepsilon}\cup\Omega_{f}^{\varepsilon}=[-\varepsilon,1]$
with boundary conditions $J(-\varepsilon)=J(1)=0$. This problem presents a
multiscale challenge because of the different spatial scales of
$\Omega_{b}^{\varepsilon}=[-\varepsilon,\varepsilon L]$ and
$\Omega_{f}^{\varepsilon}=[\varepsilon L,1]$. To overcome this difficulty, we
aim at approximating the behaviour of the trap in $\Omega^{\varepsilon}_{b}$
with a suitable boundary condition at $x=0$, obtaining then a simplified
problem in $\Omega=[0,1]$ as follows. Using a scaling variable
$\xi=1+x/\varepsilon$, the potential can be written in terms of $U(\xi)$ for
$\xi\in[0,1+L]$ as $V(x)=U(\xi)$. In summary, we first choose a scaling
potential $U(\xi)$ for $\xi\in[0,1+L]$ such that $U^{\prime}(\xi)<0$ in
$[0,1]$, $U^{\prime}(\xi)>0$ in $[1,1+L]$, $U(1+L)=0$, $U^{\prime}(1+L)=0$,
and then we study the behaviour of the trap when the potential is
$V(x)=U(1+x/\varepsilon)$. We assume that the solution
$c_{\varepsilon}(\xi,t)$ of the scaled problem

$\displaystyle\displaystyle\frac{\partial c_{\varepsilon}}{\partial
t}+\frac{1}{\varepsilon}\frac{\partial J_{\varepsilon}}{\partial\xi}$
$\displaystyle=0$ (5) $\displaystyle\displaystyle J_{\varepsilon}$
$\displaystyle=-D\frac{1}{\varepsilon}\left(\frac{\partial
c_{\varepsilon}}{\partial\xi}+\frac{1}{k_{B}T}\,c_{\varepsilon}\,U^{\prime}\right)$
(6)

has the following expansion in $\Omega^{\varepsilon}_{b}$:

$c_{\varepsilon}(\xi,t)=c^{(0)}(\xi,t)+\varepsilon
c^{(1)}(\xi,t)+O(\varepsilon^{2}).$ (7)

Since the flux $J_{\varepsilon}$ must be bounded for $\varepsilon\to 0$, from
$\eqref{eq_flux_eps}$ we have that the coefficient of the term
$\mathcal{O}(\varepsilon^{-1})$ in $J_{\varepsilon}$ has to vanish:

$\frac{\partial
c^{(0)}}{\partial\xi}+\frac{1}{k_{B}T}U^{\prime}(\xi)c^{(0)}=0.$ (8)

This equation can be solved for $c^{(0)}(\xi,t)$, yielding

$c^{(0)}(\xi,t)=c^{(0)}(1+L,t)\exp\left(-\frac{U(\xi)}{k_{B}T}\right)$ (9)

since $U(1+L)=0$. Integrating (3) in $\Omega^{\varepsilon}_{b}$ we have:

$\frac{d}{dt}\int_{-\varepsilon}^{\varepsilon L}c(x,t)\,dx+J(\varepsilon
L)-J(-\varepsilon)=0$

and using the approximation $c(x,t)\approx c^{(0)}(\xi,t)$, the boundary
condition $J(-\varepsilon)=0$ and that $V^{\prime}(\varepsilon L)=0$ we obtain

$\displaystyle\varepsilon\;\frac{\partial c(\varepsilon L,t)}{\partial
t}\;\int_{0}^{1+L}\exp\left(-\frac{U(\xi)}{k_{B}T}\right)d\xi-D\frac{\partial
c(\varepsilon L,t)}{\partial x}$ $\displaystyle=0$

that represents a boundary condition of $c(x,t)$ at $x=\varepsilon L$. Using
this boundary condition at $x=0$ instead of $x=\varepsilon L$, we finally
obtain the following simplified problem related to the multiscale model:

$\displaystyle\frac{\partial c}{\partial t}$
$\displaystyle=\frac{\partial}{\partial x}\left(D\frac{\partial c}{\partial
x}\right)\quad{\rm{for}}\,x\in[0,1]$ (10) $\displaystyle\frac{\partial
c}{\partial x}$ $\displaystyle=0\quad{\rm{at}}\,x=1$ (11) $\displaystyle
M\frac{\partial c}{\partial t}-D\frac{\partial c}{\partial x}$
$\displaystyle=0\quad{\rm{at}}\,x=0$ (12)

where

$M=\varepsilon\int_{0}^{1+L}\exp\left(-\frac{U(\xi)}{k_{B}T}\right)d\xi.$ (13)

We observe that if the potential does not depend on $\varepsilon$,
$M\rightarrow 0$ as $\varepsilon\rightarrow 0$ and then the condition (11)
reduces to a zero Neumann boundary condition, therefore the interesting
multiscale limit is obtained by letting $\varepsilon\to 0$, still maintaining
$M$ finite.111The effective dependence of the small still finite size of
$\varepsilon$ is studied in [1].

Figure 3: Representation of the domain on the left and classification of
inside grid points (green), ghost points (red) and inactive points (blue
circles) on the right.

Let us now describe the problem in higher dimensions. In two space dimensions,
the fluid is contained in a domain $\Omega$ which is the a region external to
a bubble $\mathcal{B}$ and internal to the square box
$\mathcal{S}=(-a,a)^{2}\subset\mathbb{R}^{2}$, with $a>0$ (see Figure (3),
left panel). Eq. (10) reads:

$\frac{\partial c}{\displaystyle\partial t}=\nabla\cdot\left(D\nabla
c\right)\text{ in }\Omega$ (14)

where $D$ is the diffusion coefficient. Imposing zero flux at the wall results
in homogeneous Neumann boundary conditions on
$\Gamma_{\mathcal{S}}=\partial\mathcal{S}$:

$\frac{\partial c}{\partial{n}_{\mathcal{S}}}=0\quad\text{ on
}\Gamma_{\mathcal{S}},$ (15)

where $n_{\mathcal{S}}$ denotes the unit normal vector on
$\Gamma_{\mathcal{S}}$, pointing out of the domain $\Omega$.

In the presence of a steady bubble, the fluid domain is represented by
$\Omega=\mathcal{S}\backslash\mathcal{B}$, where $\mathcal{B}$ is the region
occupied by the bubble and represented by a sphere centred in the origin and
with radius $R_{\mathcal{B}}$ such that $0<R_{\mathcal{B}}<a$. Similarly to
(12), a suitable boundary condition is enforced on the boundary
$\Gamma_{\mathcal{B}}=\partial\mathcal{B}$ to simulate the attractive-
repulsive mechanism of the bubble surface with the particles. In 2D, the
analogue of boundary condition (12) becomes (see [1] for more details):

$M\frac{\partial c}{\partial
t}=MD\frac{\partial^{2}c}{\partial{\tau}^{2}}-D\frac{\partial
c}{\partial{n}_{\mathcal{B}}}\quad\text{ on }\Gamma_{\mathcal{B}},$ (16)

where $M$ is given by the analogue of (13) (the integration of the potential
is performed along the direction normal to the surface of the bubble), $\tau$
denotes the unit vector tangential to $\Gamma_{\mathcal{B}}$,
$n_{\mathcal{B}}$ is the unit normal vector on $\Gamma_{\mathcal{B}}$ pointing
out of the domain $\Omega$ and $\partial^{k}/\partial\tau^{k}$ denotes the
$k$-th derivative along such tangential direction.

In 3D the region $\mathcal{S}$ is the cube $(-a,a)^{3}$, the static bubble
$\mathcal{B}$ is a sphere centered at the origin with radius
$R_{\mathcal{B}}<a$, and the boundary condition on $\mathcal{B}$ becomes:

$M\frac{\partial c}{\partial t}=MD\Delta_{\perp}c-D\frac{\partial c}{\partial
n_{\mathcal{B}}}\quad\text{ in }\Gamma_{\mathcal{B}}$ (17)

where $\Delta_{\perp}$ is the Laplacian-Beltrami operator on the surface of
the bubble $\Gamma_{\mathcal{B}}$.

##  3 Finite-Difference discretization

###  3.1 Discretization in time

Eqs. (14) and (16) can be written in compact form

$\frac{\partial c}{\partial t}=Q\,c$ (18)

where $Q$ is the following (linear) differential operator

$Q\,c=\left\\{\begin{matrix}D\Delta c&\text{ in }\Omega\\\ \\\ \displaystyle
D\frac{\partial^{2}c}{\partial\tau^{2}}-DM^{-1}\frac{\partial c}{\partial
n}\quad&\text{ on }\Gamma_{\mathcal{B}}\end{matrix}\right.$ (19)

with homogeneous Neumann boundary condition (15) on $\Gamma_{\mathcal{S}}$.

Eq. (18) is discretized in time by using the Crank-Nicolson method, which is
second order accurate:

$\displaystyle\frac{c^{n+1}-c^{n}}{k}$
$\displaystyle=\frac{1}{2}\left(Q\,c^{n}+Q\,c^{n+1}\right)$
$\displaystyle\left(I-\frac{k}{2}Q\right)\,c^{n+1}$
$\displaystyle=\left(I+\frac{k}{2}Q\right)\,c^{n}$ (20)

where $k$ is the time step and $I$ is the identity operator.

###  3.2 Discretization in space

The computational domain $\mathcal{S}$ is discretized through a uniform
Cartesian mesh with spatial step $h=2a/N=\Delta x=\Delta y$, where $N^{2}$ is
the number of cells. We choose to use a cell-centered discretization to
facilitate the implementation of homogeneous Neumann boundary conditions on
$\Gamma_{\mathcal{S}}$. However, the accuracy of the method does not rely on
this choice and a vertex-centered discretization would produce similar
results. Therefore, the set of grid points is
$\mathcal{S}_{h}=\\{(x_{i},y_{j})=(-a-h/2+ih,-a-h/2+jh),(i,j)\in\\{1,\cdots,N\\}^{2}\\}$.
Within the set of grid points we define the set of internal points
$\Omega_{h}=\mathcal{S}_{h}\cap\Omega$, the set of bubble points
$\mathcal{B}_{h}=\mathcal{S}_{h}\cap\mathcal{B}$ and the set of ghost points
$\mathcal{G}_{h}$ as grid points inside the bubble with at least one neighbor
point inside $\Omega$:

$(x_{i},y_{j})\in\mathcal{G}_{h}\iff(x_{i},y_{j})\in\mathcal{B}_{h}\text{ and
}\\{(x_{i}\pm h,y_{j}),(x_{i},y_{j}\pm h)\\}\cap\Omega_{h}\neq\emptyset.$ (21)

The remaining grid points that are neither inside nor ghost are called
inactive points. See Fig. 3 (right panel) for a classification of inside,
ghost, and inactive points. Let $N_{I}=|\Omega_{h}|$ and
$N_{G}=|\mathcal{G}_{h}|$ be the number of internal and ghost points,
respectively. We aim at approximating the solution $c$ at grid points of
$\Omega_{h}\cup\mathcal{G}_{h}$, then our numerical solution can be
represented as a column vector
$c_{h}=(\ldots,c_{i,j},\ldots)^{T}\in\mathbb{R}^{N_{I}+N_{G}}$, after choosing
a bijective map between $\left\\{1,\ldots,N_{I}+N_{G}\right\\}$ and the grid
points of $\Omega_{h}\cup\mathcal{G}_{h}$ (the overall numerical method does
not rely on the particular choice of this map). The problem (20) is then
discretized in space, leading to a linear system

$\left(I_{h}-\frac{k}{2}Q_{h}\right)c_{h}^{n+1}=\left(I_{h}+\frac{k}{2}Q_{h}\right)c_{h}^{n},$
(22)

to be solved at each time step, where $I_{h}$ and $Q_{h}$ are the
$(N_{I}+N_{G})\times(N_{I}+N_{G})$ matrices representing the discretization of
the operators $I$ and $Q$ and defined as follows. We denote by
$I_{h}^{(i,j)}=\left(I_{h}^{(i,j),1},\ldots,I_{h}^{(i,j),N_{I}+N_{G}}\right)$
and
$Q_{h}^{(i,j)}=\left(Q_{h}^{(i,j),1},\ldots,Q_{h}^{(i,j),N_{I}+N_{G}}\right)$
the rows of $I_{h}$ and $Q_{h}$, respectively, associated with the grid point
$(x_{i},y_{j})$.

If $P_{ij}=(x_{i},y_{j})\in\Omega_{h}$ in an internal grid point (as in Fig.
4, left panel), then the equation of the linear system is obtained from the
discretization of the internal equation (14) and the standard central
difference on a five-point stencil is used to discretize the Laplace operator
on $(x_{i},y_{j})$, and therefore $I_{h}^{(i,j)}$ and $Q_{h}^{(i,j)}$ are
defined by

$I_{h}^{(i,j)}c_{h}=c_{i,j}$ (23)
$Q_{h}^{(i,j)}c_{h}=D\frac{c_{i+1,j}+c_{i-1,j}+c_{i,j+1}+c_{i,j-1}-4c_{i,j}}{h^{2}}.$

If $(x_{i},y_{j})$ is close the the wall and the five-point stencil contains
grid points outside $\Omega$, we can use the boundary condition (15) to reduce
the five-point stencil and use only internal grid points. For example, looking
at Fig. 4 (right panel), we use the boundary condition

$\frac{c_{i,0}-c_{i,1}}{h}=0\Longrightarrow c_{i,0}=c_{i,1}$

and then

$Q_{h}^{(i,1)}c_{h}=D\frac{c_{i+1,1}+c_{i-1,1}+c_{i,2}-3c_{i,1}}{h^{2}}.$

Figure 4: Representation of the five-point stencil for the discretization of
internal points $P_{ij}=(x_{i},y_{j})$ (left panel) and the reduced stencil
when $P_{ij}$ is close to the wall $\Gamma_{\mathcal{S}}$ (right panel). In
the latter case the stencil is composed by four points.

If $G=(x_{i},y_{j})\in\mathcal{G}_{h}$ is a ghost point, then we discretize
the boundary condition (16), following a ghost-point approach similar to the
one proposed in [27] and summarised as follows. We first compute the closest
boundary point $B\in\Gamma_{\mathcal{B}}$ by

$B=O+R_{\mathcal{B}}\frac{G-O}{|G-O|},$

where $O$ is the centre of the bubble and $R_{\mathcal{B}}$ is the radius.
Then, we identify the $3\times 3$ – point stencil having
$G=(x_{G},y_{G})=(x_{i},y_{j})$ on one corner and whose convex hull contains
$B=(x_{B},y_{B})$ (see Fig. 6, right panel):

$\left\\{(x_{i+s_{x}m_{x}},x_{j+s_{y}m_{y}})\colon
m_{x},m_{y}=0,1,2\right\\},$

where $s_{x}={\rm SGN}(x_{B}-x_{G})$ and $s_{y}={\rm SGN}(y_{B}-y_{G})$, with
${\rm SGN}(\alpha)=-1$ for $\alpha<0$ and ${\rm SGN}(\alpha)=1$ for
$\alpha\geq 0$. The solution $c$ and its first and second derivatives are then
interpolated at the boundary point $B$ using the discrete values $c_{i,j}$ on
the $3\times 3$ – point stencil. The interpolations can be obtained as tensor
products of 1D interpolations in the axis directions. In detail, the 1D
quadratic interpolations using the grid points $x_{i-2},x_{i-1},x_{i}$ to
evaluate the function, its first derivative and the second derivative on
$x_{i}-\vartheta h$ with $0\leq\vartheta<1$ (see Fig. 5) are given by

$\tilde{c}(x_{i}-\vartheta\,h)=\sum_{m=0}^{2}l_{m}(\vartheta)\,c_{i-m},\quad\tilde{c}^{\prime}(x_{i}-\vartheta\,h)=\sum_{m=0}^{2}l^{\prime}_{m}(\vartheta)\,c_{i-m},\quad\tilde{c}^{\prime\prime}(x_{i}-\vartheta\,h)=\sum_{m=0}^{2}l^{\prime\prime}_{m}(\vartheta)\,c_{i-m},$

where

$l(\vartheta)=\left(\frac{(1-\vartheta)(2-\vartheta)}{2},\quad\vartheta(2-\vartheta),\quad\frac{\vartheta(\vartheta-1)}{2}\right)$
$l^{\prime}(\vartheta)=\frac{1}{h}\left(\frac{(2\vartheta-3)}{2},\quad
2(1-\vartheta),\quad\frac{(2\vartheta-1)}{2}\right)$
$l^{\prime\prime}(\vartheta)=\frac{1}{h^{2}}\left(1,\quad-2,\quad 1\right).$

In 2D, we define (see Fig. 6, right panel)

$\vartheta_{x}=s_{x}(x_{B}-x_{G})/h,\qquad\vartheta_{y}=s_{y}(y_{B}-y_{G})/h.$

Figure 5: 1D interpolation on grid points $x_{i-2},x_{i-1},x_{i}$ (circle
markers) to evaluate the function, its first derivative and the second
derivative on $x_{i}-\vartheta h$ (star marker).

Observe that $0\leq\vartheta_{x},\vartheta_{y}<1$. The 2D interpolation
formulas are:

$\tilde{c}(B)=\sum_{m_{x},m_{y}=0}^{2}l_{m_{x}}(\vartheta_{x})l_{m_{y}}(\vartheta_{y})c_{i+s_{x}m_{x},j+s_{y}m_{y}},\\\
\frac{\partial\tilde{c}}{\partial
x}(B)=s_{x}\sum_{m_{x},m_{y}=0}^{2}l^{\prime}_{m_{x}}(\vartheta_{x})l_{m_{y}}(\vartheta_{y})c_{i+s_{x}m_{x},j+s_{y}m_{y}},\\\
\frac{\partial\tilde{c}}{\partial
y}(B)=s_{y}\sum_{m_{x},m_{y}=0}^{2}l_{m_{x}}(\vartheta_{x})l^{\prime}_{m_{y}}(\vartheta_{y})c_{i+s_{x}m_{x},j+s_{y}m_{y}},\\\
\frac{\partial^{2}\tilde{c}}{\partial
x^{2}}(B)=\sum_{m_{x},m_{y}=0}^{2}l^{\prime\prime}_{m_{x}}(\vartheta_{x})l_{m_{y}}(\vartheta_{y})c_{i+s_{x}m_{x},j+s_{y}m_{y}},\\\
\frac{\partial^{2}\tilde{c}}{\partial
y^{2}}(B)=\sum_{m_{x},m_{y}=0}^{2}l_{m_{x}}(\vartheta_{x})l^{\prime\prime}_{m_{y}}(\vartheta_{y})c_{i+s_{x}m_{x},j+s_{y}m_{y}},\\\
\frac{\partial^{2}\tilde{c}}{\partial x\partial
y}(B)=s_{x}\,s_{y}\sum_{m_{x},m_{y}=0}^{2}l^{\prime}_{m_{x}}(\vartheta_{x})l^{\prime}_{m_{y}}(\vartheta_{y})c_{i+s_{x}m_{x},j+s_{y}m_{y}}.$
(24)

Finally, the rows of $I_{h}$ and $Q_{h}$ associated with the ghost point
$G=(x_{G},y_{G})$ are defined by evaluating the boundary condition on $B$,
i.e.

$I_{h}^{(i,j)}c_{h}=\tilde{c}(B)$ (25)
$Q_{h}^{(i,j)}c_{h}=D\left.\frac{\partial^{2}\tilde{c}}{\partial\tau^{2}}\right|_{B}-\frac{D}{M}\left.\frac{\partial\tilde{c}}{\partial\bm{n}}\right|_{B}$
(26)

where

$\displaystyle\frac{\partial}{\partial n}=n_{x}\frac{\partial}{\partial
x}+n_{y}\frac{\partial}{\partial y},$
$\displaystyle\qquad\displaystyle\frac{\partial^{2}}{\partial\tau^{2}}=\displaystyle\tau_{x}^{2}\frac{\partial^{2}}{\partial
x^{2}}+2\tau_{x}\tau_{y}\frac{\partial^{2}}{\partial x\partial
y}+\tau_{y}^{2}\frac{\partial^{2}}{\partial y^{2}},$ (27)
$\displaystyle(n_{x},n_{y})=\frac{O-G}{|O-G|},$
$\displaystyle\qquad(\tau_{x},\tau_{y})=(-n_{y},n_{x}).$ (28)

Observe that, while the row $I_{h}^{(i,j)}$ associated with an internal point
$(x_{i},y_{j})$ is the row of the identity matrix (23), for a ghost point
$G=(x_{i},y_{j})$ this is not true (unless $G=B$, i.e.
$\vartheta_{x}=\vartheta_{y}=0$) since it contains $3^{2}=9$ values
$l_{m_{x}}(\vartheta_{x})l_{m_{y}}(\vartheta_{y})$ for $m_{x},m_{y}=0,1,2$.

###  3.3 Complex-shaped bubbles: a level-set approach

The discretization described in the previous sections for a spherical bubble
can be extended to the case of more complex-shaped bubbles adopting a level-
set approach. In detail, the bubble $\mathcal{B}$ can be implicit defined by a
level set function $\phi(x,y)$ that is positive inside the bubble, negative
outside and zero on the boundary $\Gamma_{\mathcal{B}}$ ([33, 34]):

$\displaystyle\mathcal{B}$ $\displaystyle=\\{(x,y):\phi(x,y)>0\\}$ (29)
$\displaystyle\Gamma_{\mathcal{B}}$ $\displaystyle=\\{(x,y):\phi(x,y)=0\\}.$
(30)

The unit normal vector $n$ in (27) can be computed by:

$n=\frac{\nabla\phi}{|\nabla\phi|}$ (31)

provided that the level-set function is known explicitly. If it is known only
at grid nodes, then the derivatives of $\phi$ in (31) are approximated by
adopting a similar interpolation procedure as the one described in Eq. (24).

We observe that for a given bubble $\mathcal{B}$ there are infinite level-set
functions. For example, $\phi=R_{\mathcal{B}}-\sqrt{x^{2}+y^{2}}$ and
$\phi=R_{\mathcal{B}}^{2}-(x^{2}+y^{2})$ describe the same circular bubble.
For a given bubble, the most convenient level-set function in terms of
numerical stability is the signed distance function $\phi_{d}(x,y)$, i.e. the
distance between $(x,y)$ and $\Gamma_{\mathcal{B}}$ (positive inside the
bubble, negative otherwise). The signed distance function $\phi_{d}$ can be
computed from a generic level-set function $\phi$ by the reinitialization
algorithm [35, 36, 37], consisting of finding the steady-state solution of:

$\frac{\partial\hat{\phi}}{\partial t}={\rm
sgn}(\phi)\left(1-|\nabla\hat{\phi}|\right),\qquad\hat{\phi}=\phi\quad\text{
at time }t=0$ (32)

where $t$ is a fictitious time. A signed distance function is preferred to
avoid numerical instabilities associated with sharp or shallow gradients close
to the boundary (for a signed distance function we have $|\nabla\phi_{d}|=1$).
However, the cases investigated in the present paper involve steady bubbles or
moving bubbles with a pre-determined evolution of the shape, then the
instability issues of a generic level-set function are not observed.

Figure 6: Representation of the five-point stencil when $P_{ij}$ is close to
the boundary $\Gamma_{\mathcal{B}}$ (left panel). In this case the stencil
contains a ghost point $G$. On the right panel we represent the upwind nine-
point stencil associated with the ghost point $G$ and the boundary projection
point $B$.

##  4 Multigrid approach

The linear system (22) can be written as $A_{h}\,c_{h}^{n+1}=b_{h}$, where
$A_{h}=\left(I_{h}-\frac{k}{2}Q_{h}\right)$ and
$b_{h}=\left(I_{h}+\frac{k}{2}Q_{h}\right)c_{h}^{n}$, and it is solved in this
paper using an efficient multigrid approach that is an extension of the method
proposed in [27] for elliptic equations on complex-shaped domains. In brief, a
multigrid method is an iterative solver that starts by performing few steps of
a suitable relaxation scheme to the linear system $A_{h}\,c_{h}^{n+1}=b_{h}$,
obtaining an approximated solution $\bar{c}^{n+1}$. The relaxation scheme is
chosen in such a way that the high frequency Fourier modes of the residual
$r_{h}=b_{h}-A_{h}\,c_{h}^{n+1}$ are dumped away much quicker than the low
frequency Fourier modes. In other words, the relaxation operator smooths the
residual $r_{h}$ after few relaxation steps (say $\nu_{1}$ steps). If so, it
is said to have the smoothing property. Then, the residual $r_{h}$ is
transferred to a coarser grid with spatial step $H=2h$ (without losing much
information, as it is mainly composed of low frequency modes) by a suitable
restriction operator $r_{H}=\mathcal{I}^{h}_{H}r_{h}$ and then the residual
equation $A_{H}e_{H}=r_{H}$ is solved on the coarse grid to obtain an
approximation of the error $e_{H}$. Then, the error is transferred to the fine
grid by an interpolation operator $e_{h}=\mathcal{I}^{H}_{h}e_{H}$ and the
approximation $\bar{c}^{n+1}$ is updated by
$\bar{c}^{n+1}\leftarrow\bar{c}^{n+1}+e_{h}$. Few more steps (say $\nu_{2}$)
of the relaxation operator are then performed on the fine grid to reduce the
errors introduced by the interpolation procedure. The entire scheme is then
performed iteratively until the residual falls below a certain tolerance. In
addition, the residual equation $A_{H}e_{H}=r_{H}$ can be solved recursively
by moving to a coarser grid with spatial step $2H$, and so on. Several types
of multigrid schemes, such as $V-$cycle, $W-$cycle and Full Multigrid, can be
adopted according to the different strategies that can be chosen. In this
paper we use a $W-$cycle approach and describe the main components of the
multigrid method, i.e. relaxation, restriction and interpolation operators,
while we refer the reader to, for example, [38] for a comprehensive treatment
of multigrid methods.

###  4.1 Relaxation scheme

Standard relaxation schemes that show the smoothing property for elliptic
equations in rectangular domains are the Gauss-Seidel scheme and the weighted
Jacobi scheme (with weight $\omega=2/3$ in 1D and $\omega=4/5$ in 2D, see [38]
). While these relaxation schemes converge when the discretization is
performed on a rectangular domain, they might not converge when using a ghost-
point approach for curved boundaries (see [27]). To obtain a convergent scheme
for the problem proposed in this paper, we modify the relaxation on ghost
points as described below (while we keep a Gauss-Seidel scheme on the internal
equations). The proposed relaxation scheme can be written in the Richardson
form (with iterative index $k$)

$c_{h}^{n+1,k+1}=c_{h}^{n+1,k}+P_{h}^{-1}(b_{h}-A_{h}c_{h}^{n+1,k})$ (33)

where $P_{h}$ is a $(N_{I}+N_{G})\times(N_{I}+N_{G})$ matrix called
preconditioner and is chosen as a suitable approximation of $A_{h}$. A
standard Gauss-Seidel scheme (on both internal and ghost points) corresponds
to $P_{h}=(D_{h}+L_{h})$, where $D_{h}$ and $L_{h}$ are the diagonal and lower
part of $A_{h}$, respectively. To modify the scheme for ghost points, we
change the $N_{G}$ diagonal values of $D_{h}$ that corresponds to ghost
points, obtaining a new diagonal matrix $\tilde{D}_{h}$. Finally, the diagonal
entries $\tilde{D}_{h}^{(i,j)}$ of $\tilde{D}_{h}$ are

$\tilde{D}_{h}^{(i,j)}=\left\\{\begin{matrix}D_{h}^{(i,j)}=1+\displaystyle\frac{2kD}{h^{2}}&\text{
if }&(x_{i},y_{j})\in\Omega_{h}\\\ \beta&\text{ if }&(x_{i},y_{j})\in\text{
Ghost }\\\ \end{matrix}\right.$

where $\beta\in\mathbb{R}$ is a suitable value that we determine later.

We observe that in practice the relaxation (33) is performed without storing
the entire matrices $P_{h}$ and $A_{h}$ (matrix-free fashion). This has the
advantage of avoiding the explicit construction of the sparse matrices with
great simplification of implementation aspects and savings in computational
time, especially for moving domains when the discrete operator depends on
time. In fact, the vector $c_{h}^{n+1}$ is computationally stored in a
temporary array $c_{h}$ and its components $c_{i,j}$ are updated (overridden)
in a Gauss-Seidel fashion by iterating all over the grid points. We
distinguish between internal and ghost points. We use the notation
$a\leftarrow b$ to say that the variable $a$ is updated with the value $b$.

##### Internal points

If $(i,j)\in\Omega_{h}$,

$c_{i,j}\leftarrow
c_{i,j}+\frac{h^{2}}{h^{2}+2kD}\left(b_{i,j}-c_{i,j}-\frac{kD}{2h^{2}}\left(c_{i+1,j}+c_{i-1,j}+c_{i,j+1}+c_{i,j-1}-4c_{i,j}\right)\right)$

##### Ghost points

If $(i,j)\in\mathcal{G}_{h}$,

$c_{i,j}\leftarrow
c_{i,j}+\beta^{-1}\left(b_{i,j}-\left(I_{h}^{(i,j)}c_{h}-\displaystyle\frac{k}{2}Q_{h}^{(i,j)}c_{h}\right)\right)$
(34)

where $I_{h}^{(i,j)}$ and $I_{h}^{(i,j)}$ are defined by (25) and (26),
respectively. The iteration (34) can be written as

$c_{i,j}\leftarrow\left(1-\beta^{-1}\left(I_{h}^{(i,j),(i,j)}-\displaystyle\frac{k}{2}Q_{h}^{(i,j),(i,j)}\right)\right)c_{i,j}+\ldots\text{
terms that do not depend on }c_{i,j}\ldots$ (35)

and the value $\beta$ is chosen in such a way that the coefficient of
$c_{i,j}$ on the right-hand side of (35) is not larger than one in absolute
value, i.e.

$\left|1-\beta^{-1}A_{h}^{(i,j),(i,j)}\right|\leq 1,\text{ with
}A_{h}^{(i,j),(i,j)}=\left(I_{h}^{(i,j),(i,j)}-\displaystyle\frac{k}{2}Q_{h}^{(i,j),(i,j)}\right).$
(36)

Using (25) and (26), we have

$I_{h}^{(i,j),(i,j)}=l_{0}(\vartheta_{x})\,l_{0}(\vartheta_{y})$ (37)
$\displaystyle Q_{h}^{(i,j),(i,j)}=$ $\displaystyle
D\left(\tau_{x}^{2}l^{\prime\prime}_{0}(\vartheta_{x})l_{0}(\vartheta_{y})+2\tau_{x}\tau_{y}l^{\prime}_{0}(\vartheta_{x})l^{\prime}_{0}(\vartheta_{y})+\tau_{y}^{2}l_{0}(\vartheta_{x})l^{\prime\prime}_{0}(\vartheta_{y})\right)$
$\displaystyle-\frac{D}{M}\left(n_{x}l^{\prime}_{0}(\vartheta_{x})l_{0}(\vartheta_{y})+n_{y}l_{0}(\vartheta_{x})l^{\prime}_{0}(\vartheta_{y})\right)$
(38)

In order to satisfy condition (36), we require that $\beta$ is chosen in such
a way that

$0\leq\beta^{-1}A_{h}^{(i,j),(i,j)}\leq 2.$ (39)

The left inequality of (39) is satisfied by choosing $\text{sign
}\beta=\text{sign }A_{h}^{(i,j),(i,j)}$ (we conventionally choose $\text{sign
}(0)=1$). To satisfy the right inequality of (39), we have

$\left|\beta\right|\geq\frac{\left|A_{h}^{(i,j),(i,j)}\right|}{2}.$

This condition is always satisfied (regardless of the ghost point) if we
choose

$\left|\beta\right|\geq\tilde{A}/2\quad\text{ with
}\quad\sup_{\vartheta_{x},\vartheta_{y}\in[0,1]}\left|A_{h}^{(i,j),(i,j)}\right|\leq\tilde{A}.$
(40)

The estimate $\tilde{A}$ can be found as follows. Using (37), (4.1) and (24),
we have

$\left|A_{h}^{(i,j),(i,j)}\right|\leq\left|I_{h}^{(i,j),(i,j)}\right|+\frac{k}{2}\left|Q_{h}^{(i,j),(i,j)}\right|\leq\left|l_{0}(\vartheta_{x})\right|\left|l_{0}(\vartheta_{y})\right|\\\
+\frac{Dk}{2}\left(\left|l^{\prime\prime}_{0}(\vartheta_{x})\right|\left|l_{0}(\vartheta_{y})\right|+2\left|l^{\prime}_{0}(\vartheta_{x})\right|\left|l^{\prime}_{0}(\vartheta_{y})\right|+\left|l_{0}(\vartheta_{x})l^{\prime\prime}_{0}(\vartheta_{y})\right|+\frac{1}{\left|M\right|}\left(\left|l^{\prime}_{0}(\vartheta_{x})\right|\left|l_{0}(\vartheta_{y})\right|+\left|l_{0}(\vartheta_{x})\right|\left|l^{\prime}_{0}(\vartheta_{y})\right|\right)\right).$
(41)

Since

$\begin{gathered}\sup_{\vartheta\in[0,1]}\left|l_{0}(\vartheta)\right|=\sup_{\vartheta\in[0,1]}\left|\frac{(1-\vartheta)(2-\vartheta)}{2}\right|=1,\quad\sup_{\vartheta\in[0,1]}\left|l^{\prime}_{0}(\vartheta)\right|=\sup_{\vartheta\in[0,1]}\left|\frac{(2\vartheta-3)}{2h}\right|=\frac{3}{2h},\\\
\sup_{\vartheta\in[0,1]}\left|l^{\prime\prime}_{0}(\vartheta)\right|=\sup_{\vartheta\in[0,1]}\left|\frac{1}{h^{2}}\right|=\frac{1}{h^{2}},\end{gathered}$
(42)

then (from (41))

$\sup_{\vartheta_{x},\vartheta_{y}\in[0,1]}\left|A_{h}^{(i,j),(i,j)}\right|\leq
1+\frac{Dk}{2}\left(\frac{13}{2h^{2}}+\frac{3}{\left|M\right|h}\right)\tilde{A}$
(43)

and finally the condition on $\left|\beta\right|$ is (from (40), taking
$\tilde{A}$ as the right side of (43))

$\left|\beta\right|\geq\frac{1}{2}\left(1+\frac{Dk}{2}\left(\frac{13}{2h^{2}}+\frac{3}{\left|M\right|h}\right)\right).$

###  4.2 Transfer operators

In this section we define the transfer operators $\mathcal{I}^{h}_{H}$
(restriction) and $\mathcal{I}^{H}_{h}$ (interpolation). We adopt a geometric
multigrid, then the Galerkin conditions (required in the Algebraic multigrid)
are not satisfied, meaning that the interpolation and restriction operators
are not one the transpose of the other (multiplied by a suitable constant) and
the coarse grid operator is constructed in the same way as in the fine grid,
without taking into account the transfer operators. In this way we can afford
treating in a simpler manner more complex-shaped geometries, maintaining the
same sparsity pattern of the system at all levels. We use a cell-centered
discretization and denote by $\Omega_{h}$ and $\Omega_{H}$ the fine and coarse
grids respectively, with $\Omega_{H}$ built in the same way as $\Omega_{h}$
(Sect. 3.2) but with spatial step $H=2h$.

####  4.2.1 Restriction operator

The defect $r_{h}$ contains both the defect of the inner relaxations and the
defect of the relaxation of the boundary conditions. Since the discrete
operators of inner equations and boundary conditions scale with different
powers of $h$, the defect may show a sharp gradient crossing the boundary. For
this reason, the restriction operator of the inner equation should involve
only inside grid points and not ghost or inactive nodes to prevent the
degradation of the multigrid performance [27, 28].

In practice, for each inner grid point $(x,y)$ of the coarse grid $\Omega_{H}$
we identify the four surrounding grid nodes of the fine grid $\Omega_{h}$,
namely

$\mathcal{N}_{(x,y)}=\left\\{\left(x\pm\frac{h}{2},y\pm\frac{h}{2}\right)\right\\}$

and we perform the restriction by averaging on the grid nodes of
$\mathcal{N}_{(x,y)}$ that are inside $\Omega$ (see Fig. 7):

$r_{H}(x,y)=\mathcal{I}^{h}_{H}r_{h}(x,y)=\frac{1}{\left|\mathcal{N}_{(x,y)}\cap\Omega_{h}\right|}\sum_{(x^{*},y^{*})\in\mathcal{N}_{(x,y)}\cap\Omega_{h}}r_{h}(x^{*},y^{*}).$
(44)

We observe that the restriction reverts to the classical restriction operator
for cell-centered discretization and rectangular domains when $(x,y)$ is away
from the boundary (see [38]):

$\mathcal{I}^{h}_{H}r_{h}(x,y)=\frac{1}{4}\left[r_{h}\left(x-\frac{h}{2},y-\frac{h}{2}\right)+r_{h}\left(x-\frac{h}{2},y+\frac{h}{2}\right)+r_{h}\left(x+\frac{h}{2},y-\frac{h}{2}\right)+r_{h}\left(x+\frac{h}{2},y+\frac{h}{2}\right)\right].$

A similar approach is adopted to compute the restriction on a ghost point
$(x,y)\in\mathcal{G}_{H}$:

$\mathcal{I}^{h}_{H}r_{h}(x,y)=\frac{1}{\left|\mathcal{N}_{(x,y)}\cap\mathcal{G}_{h}\right|}\sum_{(x^{*},y^{*})\in\mathcal{N}_{(x,y)}\cap\mathcal{G}_{h}}r_{h}(x^{*},y^{*}).$
(45)

Figure 7: Left panel: representation of the restriction operator and
respective weights when the coarse grid point is close to the boundary
$\Gamma_{\mathcal{B}}$. Right panel: restriction operator for ghost points and
respective weights.

####  4.2.2 Interpolation operator

Figure 8: Representation of the interpolation operator for cell-centered
discretization and respective weights.

Once the residual equation $A_{H}e_{H}=r_{H}$ is solved on the coarse grid,
the error $e_{H}$ is interpolated back to the fine grid
$\mathcal{I}_{h}^{H}e_{H}=e_{h}$. Since the error $e_{H}$ is continuous across
the boundary (it is the solution of the residual equation on the entire
domain), then there is no need to separate the cases for inner and ghost
points and the standard bilinear interpolation operator can be adopted across
the entire domain. For example, for the inner grid node $(x,y)$ depicted in
Fig. 8, the interpolation reads:

$e_{h}(x,y)=\mathcal{I}^{H}_{h}e_{H}(x,y)\\\
=\frac{1}{16}\left(9e_{H}\left(x-\frac{h}{2},y+\frac{h}{2}\right)+3e_{H}\left(x+\frac{h}{2},y+\frac{h}{2}\right)+3e_{H}\left(x-\frac{h}{2},y-\frac{h}{2}\right)+e_{H}\left(x+\frac{h}{2},y-\frac{h}{2}\right)\right).$
(46)

##  5 Numerical results

###  5.1 Accuracy test in 2D

In this section we test the accuracy of the method. We choose an exact
solution $c_{exa}$ and augment the system (19) as:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\frac{\partial
c}{\partial t}=D\Delta c+f\quad\text{ in }\Omega\\\
\displaystyle\frac{\partial c}{\partial n}+f_{N}=0\quad\text{ on
}\Gamma_{\mathcal{S}}\\\ \displaystyle M\frac{\partial c}{\partial
t}=MD\frac{\partial^{2}c}{\partial\tau^{2}}-D\frac{\partial c}{\partial
n_{\mathcal{B}}}+f_{B}\quad\text{ on }\Gamma_{\mathcal{B}}\\\
\end{array}\right.$ (50)

choosing $\displaystyle f,f_{N}$ and $f_{B}$ in such a way that $c=c_{exa}$ is
the exact solution. The computational domain is
$\mathcal{S}=[-1,1]\times[-1,1]$, the radius of the bubble is
$R_{\mathcal{B}}=0.4$, while $M=2\times 10^{-4}$ and $D=0.1$. We choose the
following exact solution:

$\displaystyle c_{exa}(x,y,t)$
$\displaystyle=\cos(t)^{2}c_{0}(x,y)+\sin(t)^{2}c_{1}(x,y)$ (51)
$\displaystyle c_{0}(x,y)$
$\displaystyle=\exp\left(-\frac{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}}{\sigma_{0}}\right),\quad
x_{0}=0,\quad y_{0}=-0.6,\quad\sigma_{0}=0.1$ $\displaystyle c_{1}(x,y)$
$\displaystyle=\exp\left(-\frac{\left(x-x_{1}\right)^{2}+\left(y-y_{1}\right)^{2}}{\sigma_{1}}\right),\quad
x_{1}=0,\quad y_{1}=-0.7,\quad\sigma_{1}=0.1.$

We compute the $L^{1},{L}^{2}$ and ${L}^{\infty}$ norms of the relative error
at $t=\pi/8$

$\displaystyle
e_{\gamma}=\frac{||c-c_{exa}||_{\gamma}}{||c_{exa}||_{\gamma}},\quad\gamma=1,2,\infty$
(52)

for different values of $N$ and show the results in Table 1 and in the left
panel of Fig. 10, confirming numerically that the method is second order
accurate.

$N$ | $e_{1}$ | $p_{1}$ | $e_{2}$ | $p_{2}$ | $e_{\infty}$ | ${p_{\infty}}$
---|---|---|---|---|---|---
40 | 5.165E-02 | - | 4.377E-02 | - | 5.418E-02 | -
80 | 1.234E-02 | 2.066 | 1.045E-02 | 2.066 | 1.453E-02 | 1.898
160 | 3.054E-03 | 2.014 | 2.584E-03 | 2.016 | 3.696E-03 | 1.975
320 | 7.766E-04 | 1.975 | 6.472E-04 | 1.997 | 9.323E-04 | 1.987
Table 1: Relative errors $e_{1}$, $e_{2}$ and $e_{\infty}$ at $t=\pi/8$ and
accuracy orders $p_{1}$, $p_{2}$ and $p_{\infty}$ in
$\mathcal{L}^{1},\mathcal{L}^{2}$ and $\mathcal{L}^{\infty}$ norms,
respectively, for $c$ in the 2D test of Sect. 5.1. The exact solution is (51).

###  5.2 Accuracy tests in 3D axisymmetric formulation

Figure 9: Representation of the domain $\Omega$ in 3D axisymmetric:
$\Gamma_{\mathcal{S}}$ is the external wall (top, right and bottom
boundaries); $\mathcal{B}$ is the bubble with boundary $\Gamma_{\mathcal{B}}$
and radius $R_{\mathcal{B}}$; $\Gamma_{c}$ is the axis of symmetry (left
boundary).

In this section we test the accuracy of the method for a 3D axisymmetric model
(see Fig. 9). The computational domain is $\mathcal{S}=[0,2]\times[-1,1]$, the
radius of the bubble $R_{\mathcal{B}}=0.4$ and $M=2\times 10^{-4}$. The
coordinates are the radial distance $\xi$ and the vertical coordinate $z$. and
$D=0.1$. The problem reads (see Fig. 9):

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\frac{\partial
c}{\partial
t}=D\left(\frac{\partial^{2}c}{\partial\xi^{2}}+\frac{1}{\xi}\frac{\partial
c}{\partial\xi}+\frac{\partial^{2}c}{\partial z^{2}}\right)\quad\text{ in
}\Omega\\\ \displaystyle\nabla c\cdot n_{\mathcal{S}}=0\quad\text{ on
}\Gamma_{\mathcal{S}}\cup\Gamma_{c}\\\ \displaystyle M\frac{\partial
c}{\partial t}=MD\frac{\partial^{2}c}{\partial\tau^{2}}-D\frac{\partial
c}{\partial n_{\mathcal{B}}}\quad\text{ on }\Gamma_{\mathcal{B}}\\\
\end{array}\right.$ (56)

We choose the following exact solution:

$\displaystyle c_{exa}(\xi,z,t)$
$\displaystyle=\cos(t)^{2}c_{0}(\xi,z)+\sin(t)^{2}c_{1}(\xi,z)$ (57)
$\displaystyle c_{0}(\xi,z)$
$\displaystyle=\exp\left(-\frac{\left(\xi-\xi_{0}\right)^{2}+\left(z-z_{0}\right)^{2}}{\sigma_{0}}\right),\quad\xi_{0}=0,\quad
z_{0}=-0.6,\quad\sigma_{0}=0.1$ $\displaystyle c_{1}(\xi,z)$
$\displaystyle=\exp\left(-\frac{\left(\xi-\xi_{1}\right)^{2}+\left(z-z_{1}\right)^{2}}{\sigma_{1}}\right),\quad\xi_{1}=0.1,\quad
z_{1}=-0.7,\quad\sigma_{1}=0.1$

and then we calculate the relative errors at $t=\pi/8$ as in Eq. (52). Results
are presented in table 2 and in the right panel of Fig. 10, showing second
order accuracy.

$N$ | $e_{1}$ | $p_{1}$ | $e_{2}$ | $p_{2}$ | $e_{\infty}$ | ${p_{\infty}}$
---|---|---|---|---|---|---
40 | 1.868E-01 | - | 1.233E-01 | - | 1.299E-01 | -
80 | 4.361E-02 | 2.099 | 2.885E-02 | 2.096 | 3.123E-02 | 2.057
160 | 1.075E-03 | 2.021 | 7.103E-02 | 2.022 | 7.728E-02 | 2.015
320 | 2.708E-03 | 1.989 | 1.771E-03 | 2.004 | 1.926E-03 | 2.004
Table 2: Errors $e_{1}$, $e_{2}$ and $e_{\infty}$ at $t=\pi/8$ and accuracy
orders $p_{1}$, $p_{2}$ and $p_{\infty}$ in $\mathcal{L}^{1},\mathcal{L}^{2}$
and $\mathcal{L}^{\infty}$ norms, respectively, for $c$ in the 3D axisymmetric
test of Sect. 5.2. The exact solution is Eq. (57).

Figure 10: Representation of the relative errors at $t=\pi/8$ in logarithmic
scale against the value of $N\in\\{40,80,160,320\\}$ in 2D (left panel, test
of Sect. 5.1, Table 1) and 3D axisymmetric (right panel, test of Sect. 5.2,
Table 2).

##  6 Moving bubbles

In the following sections, unless otherwise specified, we consider the 3D
axisymmetric model.

If the bubble moves over time it generates a fluid motion around it. In that
case, the fluid domain depends on time $\Omega(t)$. Particles are subjected
not only to the diffusion process but they are also transported by the moving
fluids. The motion of a fluid past an oscillating bubble is governed by the
incompressible Navier-Stokes equations. At low Reynolds numbers, the viscous
forces are dominant and the convective term of the Navier-Stokes equations can
be neglected so that the motion can then be described by the Stokes equations:

$\displaystyle\frac{\partial\textbf{u}}{\partial t}+\nabla p$
$\displaystyle=\frac{1}{Re}\nabla^{2}\textbf{u}\quad\text{ in }\Omega(t)$ (58)
$\displaystyle\nabla\cdot\textbf{u}$ $\displaystyle=0\quad\text{ in
}\Omega(t)$

where u is the fluid velocity, $p$ is the pressure, $Re$ is the Reynolds
number.

Driven by the application to Sorption Kinetics, we are interested in modelling
the fluid dynamics generated by an oscillating bubble at extremely small
amplitudes ($\sim 10^{-8}\text{ m}$), resulting in a low Reynolds number
$Re<0.1$ (see [1]). Therefore, Stokes equations provide a reasonable
approximation of the fluid dynamics for the problems investigated in this
paper.

The 3D axisymmetric formulation (the coordinates are the radial distance $\xi$
and the vertical coordinate $z$) is completed by the following boundary
conditions (see Fig. 9):

$\displaystyle\frac{\partial\textbf{u}}{\partial n}$
$\displaystyle=0\quad\text{ on }\Gamma_{c}$ (59) u $\displaystyle=0\quad\text{
on }\Gamma_{\mathcal{S}}$ $\displaystyle\textbf{u}\cdot\textbf{n}$
$\displaystyle=\textbf{u}_{b}\cdot\textbf{n}\quad\text{ on
}\Gamma_{\mathcal{B}}(t)$
$\displaystyle\frac{\partial(\textbf{u}\cdot\tau)}{\partial n}$
$\displaystyle=0\quad\text{ on }\Gamma_{\mathcal{B}}(t).$

The first one is the homogeneous boundary condition dictated by the axis
symmetry, the second one is a no-slip boundary condition on the external wall,
while the third and fourth ones are the free-slip boundary conditions at the
boundary of the bubble (n and $\tau$ are the normal and tangential vectors,
respectively), where $\textbf{u}_{b}$ is the velocity of the bubble surface.

The Stokes equations (58) are discretized in time using Crank-Nicholson:

$\displaystyle\displaystyle\frac{\textbf{u}^{(n+1)}-\textbf{u}^{(n)}}{\Delta
t}+\nabla p^{(n+1/2)}$
$\displaystyle=\frac{1}{2Re}\left(\nabla^{2}\textbf{u}^{(n)}+\nabla^{2}\textbf{u}^{(n+1)}\right)\quad\text{
in }\Omega^{(n+1)}$ (60) $\displaystyle\nabla\cdot\textbf{u}^{(n+1)}$
$\displaystyle=0\quad\text{ in }\Omega^{(n+1)}$

The pressure $p$ and the velocity components $\textbf{u}=(u,v)$ are defined on
a staggered grid (see Fig. 11): $p$ is defined at the centre of each cell,
while $u$ and $v$ are defined at the mid points of vertical and horizontal
side cells, respectively, i.e. the so called Marker-and-cell discretization
introduced by Harlow in the sixties [39].

The differential operators are discretized in space using central difference:

$\displaystyle\displaystyle\frac{\partial
p}{\partial\xi}\Big{|}_{i+1/2,j}=\frac{p_{i+1,j}-p_{i,j}}{h},\quad\frac{\partial
p}{\partial z}\Big{|}_{i,j+1/2}=\frac{p_{i,j+1}-p_{i,j}}{h}$
$\displaystyle\displaystyle\frac{\partial
u}{\partial\xi}\Big{|}_{i+1/2,j}=\frac{u_{i+3/2,j}-u_{i-1/2,j}}{2h},\quad\frac{\partial
u}{\partial z}\Big{|}_{i+1/2,j}=\frac{u_{i+1/2,j+1}-u_{i+1/2,j-1}}{2h}$
$\displaystyle\displaystyle\frac{\partial
v}{\partial\xi}\Big{|}_{i,j+1/2}=\frac{v_{i+1,j+1/2}-v_{i-1,j+1/2}}{2h},\quad\frac{\partial
v}{\partial z}\Big{|}_{i,j+1/2}=\frac{v_{i,j+3/2}-v_{i,j-1/2}}{2h}$
$\displaystyle\displaystyle\nabla^{2}u\Big{|}_{i+1/2,j}=\frac{u_{i+3/2,j}+u_{i-1/2,j}+u_{i+1/2,j+1}+u_{i+1/2,j-1}-4u_{i+1/2,j}}{h^{2}}$
$\displaystyle\displaystyle\nabla^{2}v\Big{|}_{i,j+1/2}=\frac{v_{i+1,j+1/2}+v_{i-1,j+1/2}+v_{i,j+3/2}+v_{i,j-1/2}-4v_{i,j+1/2}}{h^{2}}$
$\displaystyle\displaystyle\nabla\cdot\textbf{u}\Big{|}_{i,j}=\frac{u_{i+1/2,j}-u_{i-1/2,j}+v_{i,j+1/2}-v_{i,j-1/2}}{h}$

Figure 11: Staggered grid for the Stokes problem. Horizontal velocity $u$ is
defined on the middle points of the vertical edges of the cells (circle
markers), vertical velocity $v$ is defined on the middle points of the
horizontal edges of the cells (diamond markers), pressure $p$ is defined on
the centers of the cells (dot markers).

This discretization results in a linear system to be solved for
$(\textbf{u}^{(n+1)},p^{(n+1/2)})$. This linear system is singular, due to the
non uniqueness of $p^{(n+1/2)}$ (it is defined up to an additive constant) and
it has to satisfy a compatibility condition to guarantee the existence of the
solution. Following the approach presented in [29], the issue is circumvented
by augmenting the problem (60) with an additional scalar unknown
$\zeta\in\mathbb{R}$ and an additional equation for $p$ as follows:

$\displaystyle\displaystyle\frac{\textbf{u}^{(n+1)}-\textbf{u}^{(n)}}{\Delta
t}+\nabla p^{(n+1/2)}$
$\displaystyle=\frac{1}{2Re}\left(\nabla^{2}\textbf{u}^{(n)}+\nabla^{2}\textbf{u}^{(n+1)}\right)\quad\text{
in }\Omega^{(n+1)}$ (61) $\displaystyle\nabla\cdot\textbf{u}^{(n+1)}$
$\displaystyle=\zeta\quad\text{ in }\Omega^{(n+1)}$
$\displaystyle\int_{\Omega}p\,d\Omega$ $\displaystyle=0\quad\text{ in
}\Omega^{(n+1)}$

The problem consists then of finding $(\textbf{u}^{(n+1)},p^{(n+1/2)},\zeta)$
that satisfies Eq. (61). We observe that the free divergence condition is not
guaranteed, namely $\zeta$ is usually different from zero. However, the
divergence decades with the same order of the method, i.e. $\zeta=O(h^{2})$,
where $h$ is the spatial step, and then the overall accuracy order is not
degraded (see [29] for more details). We observed numerically that in the
absence of a bubble, namely $\Omega=\mathcal{S}$, then $\xi=0$ within machine
precision. This is usually the case of domains without curved boundaries.

The third equation of (61) is discretized by standard mid-point rule, leading
to the linear equation

$\sum_{(\xi_{i},z_{j})\text{ internal points }}p_{i,j}=0.$ (62)

The curved boundary is treated using the ghost-point technique described in
[29], similarly to the approach presented in Sect. 3.2.

###  6.1 Test 1: pulsating bubble

In this test we want to model the expansion/compression of a (pulsating)
bubble, represented by a sphere $\mathcal{B}(t)$ centred at the origin and
with radius:

$R(t)=R_{\mathcal{B}}(1+A\sin(\omega t))$ (63)

where $R_{\mathcal{B}}$ is the radius of the bubble at time $t=0$. The
velocity of the bubble surface is

$\textbf{u}_{b}(\xi,z)=R^{\prime}(t)\,\textbf{n}=A\,R_{\mathcal{B}}\,\omega\cos(\omega\,t)\,\textbf{n},\text{
where }\textbf{n}=(\xi,z)/\sqrt{\xi^{2}+z^{2}}\text{ and
}\sqrt{\xi^{2}+z^{2}}=R(t).$ (64)

The exact solution for the 3D axisymmetric Stokes problem (58) with free-slip
boundary conditions on the bubble surface (third and fourth equations of (59))
in a semi-infinite domain $\Omega(t)=\left\\{(\xi,z)\colon
0<\xi<+\infty,\right.$ $\left.\xi^{2}+z^{2}>(R(t))^{2}\right\\}$ is:

$\textbf{u}_{\rm
exa}=R^{\prime}(t)\frac{(R(t))^{2}}{(\xi^{2}+z^{2})^{3/2}}\cdot\begin{pmatrix}\xi\\\
z\end{pmatrix},\quad
p=R(t)(R^{\prime\prime}(t)R(t)+2(R^{\prime}(t))^{2})/\sqrt{\xi^{2}+z^{2}}.$
(65)

In a finite domain $\Omega(t)=\mathcal{S}\backslash\mathcal{B}(t)$ we cannot
prescribe the wall boundary conditions $\textbf{u}=0$ on the external boundary
$\Gamma_{\mathcal{S}}$ otherwise the mass conservation is not guaranteed
(since the volume of the bubble is not constant over time). For this specific
test we then prescribe the exact velocity (65) at $\Gamma_{\mathcal{S}}$. We
choose $R_{\mathcal{B}}=0.253$, $A=0.04$ and $\omega=2\pi\,\nu$ with $\nu=50$
and we compute the numerical error at time $t_{fin}=0.1$ as the difference
between the numerical solution and the exact solution (65). The results are
presented in Fig. 12, where the second order accuracy is confirmed for both
velocity components $u$ (radial component) and $v$ (vertical component) in
norms $L^{1}$, $L^{2}$ and $L^{\infty}$.

Figure 12: Representation of the relative error for the oscillating bubble,
Test 1, Sect 6.1. We plot the errors for the two components of the velocity
$u$ and $v$ in $\mathcal{L}^{1}$ and $\mathcal{L}^{\infty}$ norms. In this
test $\Omega=[0,2]\times[-1,1]$, $A$ = 0.04, $nu$ = 1, $R_{\mathcal{B}}$ =
0.253 and final time $t_{fin}$ = 0.1.

###  6.2 A steady computational bubble approach

When the amplitude of the bubble oscillation is sufficiently small compared to
its dimensions, then $R(t)\approx R_{\mathcal{B}}$ and it is reasonable to
simplify the model by assuming that the velocity of the surface bubble is
assigned at a distance of $R_{\mathcal{B}}$ from the origin rather than
$R(t)$:

$\textbf{u}_{b}(\xi,z)=A\,\omega\cos(\omega\,t)\cdot\begin{pmatrix}\xi\\\
z\end{pmatrix}\text{ for }\sqrt{\xi^{2}+z^{2}}=R_{\mathcal{B}}.$ (66)

In this way, the computational domain does not move in time, since
$\mathcal{B}(t)=\mathcal{B}(0)$ for any $t>0$ (steady computational bubble)
and the fluid motion is generated purely from the boundary conditions. The
exact solution for the 3D axisymmetric Stokes problem (58) with free-slip
boundary conditions on the bubble surface (third and fourth equations of
(59)), with surface velocity defined by (66), in a semi-infinite domain
$\Omega=\left\\{0<\xi<+\infty,\xi^{2}+z^{2}>R_{\mathcal{B}}^{2}\right\\}$ is:

$\textbf{u}^{\rm f}_{\rm
exa}=R^{\prime}(t)\frac{(R_{\mathcal{B}})^{2}}{(\xi^{2}+z^{2})^{3/2}}\cdot\begin{pmatrix}\xi\\\
z\end{pmatrix},\quad
p=R^{\prime\prime}(t)R_{\mathcal{B}}^{2}/\sqrt{\xi^{2}+z^{2}}.$ (67)

The difference between the exact solutions (65) and (67) is $\textbf{u}_{\rm
exa}-\textbf{u}^{\rm f}_{\rm exa}=\mathcal{O}(A)$. Therefore, for a fixed
spatial step, the difference between the two approaches decades as
$A\rightarrow 0$. This is confirmed numerically in Fig. 13 (left panel), where
we compute the difference between the numerical solutions of the two
approaches $\textbf{u}_{h}$ and $\textbf{u}_{h}^{f}$ at a fixed value of the
spatial step $h=1/50$ and different values of $A$ over a period $T=1/\nu$ as
follows:

$\displaystyle\frac{1}{T}\int_{0}^{T}\frac{\int_{\Omega(t)}\left|\textbf{u}_{h}-\textbf{u}^{f}_{h}\right|^{p}}{\int_{\Omega(t)}\left|\textbf{u}_{h}^{f}\right|^{p}}dt,\quad
p=1,2,\infty.$ (68)

In the right panel of Fig. 13 we compute the numerical error as the difference
between the numerical and the exact solutions for both approaches at a fixed
spatial step and different values of $A$. We observe that the second approach
generally over performs the first one, although for sufficiently small values
of $A$ the two errors reach a plateau, meaning that the overall error is
dominated by the discretization error at the fixed spatial step. We conclude
that when $A/R_{\mathcal{B}}\ll 1$ (as the cases investigated in this paper)
it is more efficient to keep a steady computational bubble $\mathcal{B}$ and
simulate the bubble motion by assigning a varying velocity at the
computational surface $\partial\mathcal{B}$. We will follow this approach in
the following tests.

Figure 13: Left panel: relative difference between the numerical solutions of
the two approaches described in Sections 6.1 and 6.2 against the value of $A$.
Right panel: relative errors of the two approaches computed as the difference
between the numerical solutions and the respective exact solutions (65) and
(67) against the value of $A$. In both plots we have
$A\in\\{h/64,h/16,h/4,h,4h,16h\\}$, with $h=1/120$, and the time range is
$[t_{\rm in},t_{\rm fin}]=[1,2]$, with $\nu=1$.

###  6.3 Test 2: Oscillating bubble

In the following tests we model the advection-diffusion process of particles
in a moving fluid past an oscillating bubble:

$\frac{\partial c}{\displaystyle\partial t}=\nabla\cdot\left(D\nabla
c-c\textbf{u}\right)$ (69)

where $\textbf{u}=(u,v)$ is the solution of the Stokes problem (58). In
general, we describe the motion of a bubble (that is not necessarily a sphere)
by its parametric equations:

$\displaystyle\xi(\theta,t)$
$\displaystyle=\xi_{c}(t)+\delta_{\xi}(t)\cos(\theta)$ $\displaystyle
z(\theta,t)$ $\displaystyle=z_{c}(t)+\delta_{z}(t)\sin(\theta)$

where $(\xi_{c}(t),z_{c}(t))$ is the centre of the bubble, while
$\delta_{\xi}(t)$ and $\delta_{z}(t)$ regulate the deformation from a
spherical shape. We assume that at time $t=0$ the bubble is a sphere centred
at the origin and with radius $R_{\mathcal{B}}$, then $\xi_{c}(0)=z_{c}(0)=0$
and $\delta_{\xi}(0)=\delta_{z}(0)=R_{\mathcal{B}}$.

Following the same approach as in Sect. 6.2, we keep a steady computational
bubble $\mathcal{B}(t)=\mathcal{B}(0)$ for $t>0$ and we model the velocity of
the surface $\textbf{u}_{b}(\xi,z)=(u_{b}(\xi,z),v_{b}(\xi,z))$ at
$(\xi,z)\colon\sqrt{\xi^{2}+z^{2}}=R_{\mathcal{B}}$ as:

$\displaystyle u_{b}(\xi,z)$
$\displaystyle=\xi^{\prime}_{c}(t)+\delta^{\prime}_{\xi}(t)\cos(\theta)$
$\displaystyle v_{b}(\xi,z)$
$\displaystyle=z^{\prime}_{c}(t)+\delta^{\prime}_{z}(t)\sin(\theta),\quad\text{
where }\quad\theta=\arctan(z/\xi).$

In Test2a, we model an harmonic vertical oscillation of the spherical bubble:

$\xi(t)=0,\quad z(t)=A\sin(2\pi\nu
t),\quad\delta_{\xi}(t)=\delta_{z}(t)=R_{\mathcal{B}},$

while In Test2b, we model an ellipsoidal deformation of the bubble:

$\xi(t)=z(t)=0,\quad\delta_{z}(t)=R_{\mathcal{B}}(1+A\sin(2\pi\nu
t)),\quad\delta_{\xi}(t)=\sqrt{R_{\mathcal{B}}^{3}/\delta_{z}(t)}.$ (70)

In (70), we observe that $\delta_{\xi}(t)$ has been defined to guarantee that
the volume of the ellipsoid, $V(t)=4/3(\pi\delta_{\xi}(t)^{2}\delta_{z}(t))$,
is constant over time. We choose $R_{\mathcal{B}}=0.258$, $A=0.01$ and
$\nu=10$ (Test2a10 and Test2b10) or $\nu=1000$ (Test2a1000 and Test2b1000). In
Fig. 14 we plot the vector fields of the fluid velocity at selected fractions
of the first oscillation period $T$ ($\nu\cdot t=0.25,0.50,0.75,1.00$). The
colormap represents the magnitude of the velocity, while the red dashed line
is the fictitious representation of the bubble (where $A$ has been amplified
by a factor of 20 for graphical purposes).

In Test2a10 we observe that a small vortex is generated next to the bubble,
moving farther to the right after disappearing at around $\xi=1.5$ between
$t\cdot\nu=0.25$ and $t\cdot\nu=0.50$. At that time, a new vortex is generated
next to the bubble, disappearing between $t\cdot\nu=0.75$ and $t\cdot\nu=1.00$
and so on, approaching a periodic behaviour. A similar mechanism is observed
in Test2a1000, except that the vortexes disappear when they are much closer to
the bubble compared to Test2a10, say at around $\xi=0.5$. In Test2b10 two
vortexes are generated at the same time, moving towards the top right and
bottom right corners of the domain, respectively. They disappear in favour of
novel vortexes with the same timeline as in Test2a10. This phenomenon is
observed in experimental results [40]. In Test2b1000 a similar behaviour is
observed, except that the two vortexes disappear when they are much closer to
the bubble than in Test2b10.

Figure 14: Vector fields of the velocity at selected fractions of the first
oscillation period $T$ ($\nu\cdot t=0.25,0.50,0.75,1.00$) for the numerical
tests of Sect. 6.3. Each row of plots corresponds to a numerical test. From
top to bottom: Test2a10, Test2a1000, Test2b10 and Test2b1000. The red dashed
line is the fictitious representation of the bubble (where $A$ has been
amplified by a factor 20 for graphical purposes).

In Table 3 we show that $\nabla\cdot\vec{u}_{h}$ decades with the same order
of accuracy of the numerical error on $\vec{u}_{h}$, i.e.
$\nabla\cdot\vec{u}_{h}=\mathcal{O}(h^{2})$, where $h$ is the spatial step.
For this purpose, we use Test2a1000 and Test2b1000. Since we do not know the
exact solution for these tests, we approximate the order of accuracy $q$ of
$\vec{u}_{h}$ by using the Richardson extrapolation:

$q\approx\log_{2}\left(\frac{\left\|\vec{u}_{h}-\vec{u}_{h/2}\right\|_{\infty}}{\left\|\vec{u}_{h/2}-\vec{u}_{h/4}\right\|_{\infty}}\right).$

We observe numerically that $q\approx 2$. The relative error on $\vec{u}_{h}$
is then approximated by

${e}_{h}=\frac{\left\|\vec{u}_{h}-\vec{u}_{\text{exa}}\right\|_{\infty}}{\left\|\vec{u}_{\text{exa}}\right\|_{\infty}}\approx\frac{4}{3}\frac{\left\|\vec{u}_{h}-\vec{u}_{h/2}\right\|_{\infty}}{\left\|\vec{u}_{h}\right\|_{\infty}}$

The relative error on $\nabla\cdot\vec{u}_{h}$ is computed by normalization
with $\nabla\vec{u}_{h}$:

${e}^{\text{div}}_{h}=\frac{\left\|\nabla\cdot\vec{u}_{h}\right\|_{\infty}}{\left\|\nabla\vec{u}_{h}\right\|_{\infty}}$

with
$\left\|\nabla\vec{u}_{h}\right\|_{\infty}=\max\left\\{\left\|\,\left|\nabla
u_{h}\right|\,\right\|_{\infty},\left\|\,\left|\nabla
v_{h}\right|\,\right\|_{\infty}\right\\}$, where $\left|\nabla u_{h}\right|$
and $\left|\nabla v_{h}\right|$ represent the central finite-difference
approximations of $(u_{x}^{2}+u_{y}^{2})^{1/2}$ and
$(v_{x}^{2}+v_{y}^{2})^{1/2}$, respectively. The order of accuracy
$q^{\text{div}}$ is approximated by

$q^{\text{div}}\approx\log_{2}\left(\frac{e_{h}^{\text{div}}}{e_{h/2}^{\text{div}}}\right)$

We also observe from Table 3 that the relative error on the divergence
${e}^{\text{div}}_{h}$ is about one order of magnitude smaller than the
relative error on velocity ${e}_{h}$.

No. of points | $\left\|\vec{u}_{h}-\vec{u}_{h/2}\right\|_{\infty}$ | order $q$ | $e_{h}$ | $e_{h}^{\text{div}}$ | order $q^{\text{div}}$
---|---|---|---|---|---
32 $\times$ 32 | 2.12 $\cdot 10^{-1}$ | 1.79 | 1.62 $\cdot 10^{-3}$ | 2.02 $\cdot 10^{-4}$ | 1.74
64 $\times$ 64 | 6.15 $\cdot 10^{-2}$ | 2.11 | 4.69 $\cdot 10^{-4}$ | 6.05 $\cdot 10^{-5}$ | 1.96
128 $\times$ 128 | 1.42 $\cdot 10^{-2}$ | 2.02 | 1.08 $\cdot 10^{-4}$ | 1.55 $\cdot 10^{-5}$ | 1.99
256 $\times$ 256 | 3.50 $\cdot 10^{-3}$ | - | 2.67 $\cdot 10^{-5}$ | 3.91 $\cdot 10^{-6}$ | 1.92
512 $\times$ 512 | - | - | - | 1.03 $\cdot 10^{-6}$ | -
No. of points | $\left\|\vec{u}_{h}-\vec{u}_{h/2}\right\|_{\infty}$ | order $q$ | $e_{h}$ | $e_{h}^{\text{div}}$ | order $q^{\text{div}}$
32 $\times$ 32 | 4.51 $\cdot 10^{-1}$ | 1.88 | 3.44 $\cdot 10^{-3}$ | 6.12 $\cdot 10^{-4}$ | 1.79
64 $\times$ 64 | 1.23 $\cdot 10^{-1}$ | 1.89 | 9.36 $\cdot 10^{-4}$ | 1.77 $\cdot 10^{-4}$ | 1.88
128 $\times$ 128 | 3.32 $\cdot 10^{-2}$ | 1.94 | 2.53 $\cdot 10^{-4}$ | 4.84 $\cdot 10^{-5}$ | 1.84
256 $\times$ 256 | 8.68 $\cdot 10^{-3}$ | - | 6.61 $\cdot 10^{-5}$ | 1.35 $\cdot 10^{-5}$ | 1.93
512 $\times$ 512 | - | - | - | 3.55 $\cdot 10^{-6}$ | -
Table 3: Relative numerical errors $e_{h}$ (on $\vec{u}_{h}$) and
$e_{h}^{\text{div}}$ (on $\nabla\cdot\vec{u}_{h}$) and respective orders of
accuracy $q$ and $q^{\text{div}}$ for Test2a1000 (top) and Test2b1000
(bottom).

In Fig. 15 we plot the particle concentration $c$ at a specific detector point
$(\xi_{d}=0.4,z_{d}=0)$ over time $t$ for different numerical tests
(Test2a10,Test2a1000,Test2b10,Test2b1000). At the initial time $t=0$ the
particles follow a Gaussian distribution centered in $(\xi_{0}=0.8,z_{0}=0)$,
as follows (initial condition):

$c(\xi,z,0)=a_{1}e^{-a_{2}((\xi-\xi_{0})^{2}+(z-z_{0})^{2})},\quad
a_{1}=\left(2\pi\sigma^{2}\right)^{-3/2},\quad
a_{2}=\left(2\sigma^{2}\right)^{-1},\sigma=0.4.$ (71)

The black line refers to the case of steady bubble (${\bf u}=0$). All tests
present an oscillating behaviour of the particle concentration in the vicinity
of the black line. The oscillation frequency of the particle concentration is
strictly related with the oscillation frequency of the bubble, while the
amplitude depends on the type of bubble oscillation (harmonic or ellipsoidal).
However, the temporal average of the particle concentration of the proposed
tests, namely
$\overline{c}(x,t)=\frac{1}{T}\int_{t-T/2}^{t+T/2}c(x,\tau)\,d\tau=\nu\int_{t-1/(2\nu)}^{t+1/(2\nu)}c(x,\tau)\,d\tau$,
does not seem to match the green line, suggesting that the bubble oscillation
actually changes the particle distribution in the vicinity of the bubble. The
problem presents a temporal multiscale effect and a rigorous mathematical
explanation of this phenomenon should provide the associated PDEs for
$\overline{c}$, showing that these equations differ from the simple diffusion
equations with ${\bf u}=0$. This is part of our ongoing effort.

Figure 15: Detector values of the particle concentration $c$ of Eq. (69) at
$(\xi_{d}=0.4,z_{d}=0)$. Initial condition is (71). The spatial step is
$h=1/120$. On the left we plot the comparison between Test2a10 (blue line),
Test2b10 (red line) and steady-bubble case with $\textbf{u}=0$ (black line).
Analogously, on the right, we show the comparison between Test2a1000 and
Test2b1000. The dashed lines represent the mean values of the respective
tests.

##  7 Conclusions

We have presented a second order accurate numerical method for the recently
developed multiscale model of sorption kinetic [1], in the single carrier
approximation. The problem consists of a concentration of particles that
diffuse in a fluid agitated by an oscillating bubble. In addition, the bubble
attracts the particles that are in the vicinity of its surface. This
attraction is modelled by a time-dependent boundary condition on the bubble
surface, where normal and (second order) tangential derivatives are involved.
The region occupied by the bubble is implicitly described by a level-set
function and the boundary conditions on the curved boundary (bubble surface)
are discretized by a proper ghost-point technique. A multigrid method is
designed to efficiently solve the linear system arising from the
discretization of the equations. The complexity of the boundary condition on
the bubble surface leads to a specific stability condition that must be
satisfied by the relaxation scheme of the multigrid method.

The fluid motion is modelled by the Stokes equations, solved by a monolithic
approach (continuity and momentum equations are solved simultaneously). A
simplified model is provided for the treatment of the boundary conditions in
the case of moving bubble, in which the bubble is computationally steady and
its oscillations are modeled by a time-dependent fluid velocity imposed on its
surface. This simplification is justified by the small amplitude of the bubble
oscillations, and a comparison with the more realistic moving domain model
confirms that the differences between the two approaches are negligible for
the problems investigated in this paper. Furthermore, this approximation makes
the computation more efficient, because it avoids evolving the domain and its
discretization, so the coefficients of the linear system (22) are time
independent.

Two types of bubble oscillations are implemented: harmonic and ellipsoidal
oscillation, the latter providing a better representation of the experimental
results existing in literature.

The particle concentration at a specific point of the domain is reconstructed
over time for the two types of oscillation. The same test is repeated for a
steady fluid case (actual steady bubble). We observed that the particle
concentration for the oscillating bubble oscillates around an average function
that is close to (but not the same as) the one obtained for the actual steady
bubble. The discrepancy between the two functions can be mathematically
described by numerical approaches developed for temporal multiscale problems.
A more rigorous analysis of this phenomenon is part of our ongoing effort.
Future works will include a saturation effect (already modelled in [1], when
high concentrations is reached near the surface of the bubble) and a full two
carrier model that describes the interaction between the two species of ions
(with a potential obtained from a self-consistent Poisson equation).

## Acknowledgments

G.R. and C.A. acknowledge partial support from ITN-ETN Horizon 2020 Project
ModCompShock, Modeling and Computation on Shocks and Interface (Project
Reference 642768), and from PRIN Project 2017 entitled ”Innovative numerical
methods for evolutionary partial differential equations and applications”
(No.2017KKJP4X), coming from the Italian Ministry of Education, University and
Research (MIUR). All authors acknowledge support from GNCS–INDAM (National
Group for Scientific Computing, Italy).

## References

  * [1] C. Astuto, A. Raudino, and G. Russo.  Multiscale modeling of sorption kinetics.  SIAM Multiscale Modeling and Simulation, submitted arXiv:2202.02552.
  * [2] Mahir Hadžić, Gustavo Navarro, and Steve Shkoller.  Local well-posedness and global stability of the two-phase stefan problem.  SIAM Journal on Mathematical Analysis, 49(6):4942–5006, 2017.
  * [3] Charles M. Elliott, Björn Stinner, and Chandrasekhar Venkataraman.  Modelling cell motility and chemotaxis with evolving surface finite elements.  Journal of The Royal Society Interface, 9(76):3027–3044, 2012.
  * [4] C. Varea, J. L. Aragón, and R. A. Barrio.  Turing patterns on a sphere.  Phys. Rev. E, 60:4588–4592, Oct 1999.
  * [5] Harald Garcke, Kei Fong Lam, and Björn Stinner.  Diffuse interface modelling of soluble surfactants in two-phase flow.  arXiv e-prints, 2013.
  * [6] Antonio Raudino, Antonio Grassi, Giuseppe Lombardo, Giovanni Russo, Clarissa Astuto, and Mario Corti.  Anomalous sorption kinetics of self-interacting particles by a spherical trap.  Communications in Computational Physics, 31(3):707–738, 2022.
  * [7] A. Raudino, D. Raciti, A. Grassi, M. Pannuzzo, and M. Corti.  Oscillations of bubble shape cause anomalous surfactant diffusion: Experiments, theory, and simulations.  Langmuir, 2016.
  * [8] A. Grassi, G.M. Lombardo, M. Pannuzzo, and A. Raudino.  Capture rate and efficiency of an oscillating non-ideal trap interacting with a sea of random diffusing particles. a non-equilibrium fokker-planck picture.  Physics Letters A, 2015.
  * [9] Antonio Raudino, Domenica Raciti, and Antonio Grassi.  Modeling the capture rate by a radially oscillating spherical bubble. a bio-mimetic model for studying the mechanically-mediated uptake by cells.  Physica A: Statistical Mechanics and its Applications, 2016.
  * [10] M. Corti, M. Bonomo, and A. Raudino.  New interferometric technique to evaluate the electric charge of gas bubbles in liquids.  Langmuir, 2012.
  * [11] M. Corti, M. Pannuzzo, and A. Raudino.  Out of equilibrium divergence of dissipation in an oscillating bubble coated by surfactants.  Langmuir, 2014.
  * [12] M. Corti, M. Pannuzzo, and A. Raudino.  Trapping of sodium dodecyl sulfate at the air-water interface of oscillating bubbles.  Langmuir, 2015.
  * [13] A. Raudino, D. Raciti, and M. Corti.  Anomalous behavior of ultra-low-amplitude capillary waves. A glimpse of the viscoelastic properties of interfacial water?  Langmuir, 2017.
  * [14] P. Brocca, A. Saponaro, B. Introini, V. Rondelli, M. Pannuzzo, D. Raciti, M. Corti, and A. Raudino.  Protein adsorption at the air-water interface by a charge sensing interferometric technique.  Langmuir, 2019.
  * [15] WH Wu, DG Eskin, A Priyadarshi, T Subroto, I Tzanakis, and W Zhai.  New insights into the mechanisms of ultrasonic emulsification in the oil–water system and the role of gas bubbles.  Ultrasonics Sonochemistry, 73:105501, 2021.
  * [16] Anastasia V Tyurnina, Justin A Morton, Tungky Subroto, Mohammad Khavari, Barbara Maciejewska, Jiawei Mi, Nicole Grobert, Kyriakos Porfyrakis, Iakovos Tzanakis, and Dmitry G Eskin.  Environment friendly dual-frequency ultrasonic exfoliation of few-layer graphene.  Carbon, 185:536–545, 2021.
  * [17] A. F. H. Ward and L. Tordai.  Time‐dependence of boundary tensions of solutions i. the role of diffusion in time‐effects.  The Journal of Chemical Physics, 14(7):453–461, 1946.
  * [18] H Diamant and D Andelman.  Kinetics of surfactant adsorption at fluid/fluid interfaces: non-ionic surfactants.  Europhysics Letters (EPL), 34(8):575–580, jun 1996.
  * [19] AFH Ward and L Tordai.  Time-dependence of boundary tensions of solutions i. the role of diffusion in time-effects.  The Journal of Chemical Physics, 1946.
  * [20] Padilla P. et al. Plaza R., Sánchez-Garduño F.  The effect of growth and curvature on pattern formation.  Journal of Dynamics and Differential Equations, 2004.
  * [21] R.P. Fedkiw, T. Aslam, B. Merriman, and S. Osher.  A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method).  Journal of Computational Physics, 152:457–492, 1999.
  * [22] F. Gibou and R. Fedkiw.  A second-order-accurate symmetric discretization of the Poisson equation on irregular domains.  Journal of Computational Physics, 176:205–227, 2002.
  * [23] F. Gibou and R. Fedkiw.  A fourth order accurate discretization for the Laplace and heat equations on arbitary domains, with applications to the Stefan problem.  Journal of Computational Physics, 202:577–601, 2005.
  * [24] Yen Ting Ng, Chohong Min, and Frédéric Gibou.  An efficient fluid-solid coupling algorithm for single-phase flows.  J. Comput. Phys., 228:8807–8829, 2009.
  * [25] Javier Fernández-Fidalgo, Stéphane Clain, Luis Ramírez, Ignasi Colominas, and Xesús Nogueira.  Very high-order method on immersed curved domains for finite difference schemes with regular cartesian grids.  Computer Methods in Applied Mechanics and Engineering, 360:112782, 2020.
  * [26] Stéphane Clain, Diogo Lopes, and Rui MS Pereira.  Very high-order cartesian-grid finite difference method on arbitrary geometries.  Journal of Computational Physics, 434:110217, 2021.
  * [27] Armando Coco and Giovanni Russo.  Finite-difference ghost-point multigrid methods on cartesian grids for elliptic problems in arbitrary domains.  Journal of Computational Physics, 241:464–501, 2013.
  * [28] Armando Coco and Giovanni Russo.  Second order finite-difference ghost-point multigrid methods for elliptic problems with discontinuous coefficients on an arbitrary interface.  Journal of Computational Physics, 361:299 – 330, 2018.
  * [29] Armando Coco.  A multigrid ghost-point level-set method for incompressible navier-stokes equations on moving domains with curved boundaries.  Journal of Computational Physics, 418:109623, 2020.
  * [30] Alina Chertock, Armando Coco, Alexander Kurganov, and Giovanni Russo.  A second-order finite-difference method for compressible fluids in domains with moving boundaries.  Commun Comput Phys, 23:230–263, 2018.
  * [31] Christiana Mavroyiakoumou and Silas Alben.  Large-amplitude membrane flutter in inviscid flow.  Journal of Fluid Mechanics, 891, 2020.
  * [32] Christiana Mavroyiakoumou and Silas Alben.  Dynamics of tethered membranes in inviscid flow.  Journal of Fluids and Structures, 107:103384, 2021.
  * [33] Stanley Osher and Ronald Fedkiw.  Level set methods and dynamic implicit surfaces, volume 153 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2003.
  * [34] J. A. Sethian.  Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science.  Cambridge monographs on applied and computational mathematics 3. Cambridge University Press, 2nd ed edition, 1999.
  * [35] Mark Sussman, Peter Smereka, and Stanley Osher.  A level set approach for computing solutions to incompressible two-phase flow.  Journal of Computational physics, 114(1):146–159, 1994.
  * [36] Giovanni Russo and Peter Smereka.  A remark on computing distance functions.  Journal of computational physics, 163(1):51–67, 2000.
  * [37] Antoine Du Chéné, Chohong Min, and Frédéric Gibou.  Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes.  Journal of Scientific Computing, 35(2):114–131, 2008.
  * [38] U. Trottemberg, C.W. Oosterlee, and A. Schuller.  Multigrid.  Academic Press, 2000.
  * [39] Francis H Harlow and J Eddie Welch.  Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface.  The physics of fluids, 8(12):2182–2189, 1965.
  * [40] Paul Tho, Richard Manasseh, and Andrew Ooi.  Cavitation microstreaming patterns in single and multiple bubble systems.  Journal of fluid mechanics, 576:191–233, 2007.