Low Order Finite Element Methods for the Navier-Stokes-Cahn-Hilliard Equations

Aleksander Lovrićcor1
[cor1]Corresponding author.

Wulf G. Dettmer
Djordje Perić

Zienkiewicz Centre for Computational Engineering, College of Engineering, Swansea University, Fabian Way, Swansea SA1 8EN, Wales, UK

A computationally efficient, low order finite element formulation is developed for modelling the Navier-Stokes-Cahn-Hilliard equations, which have been established as a promising phase field modelling approach for simulation of immiscible multiphase flows. The present study suggests that traditional Navier-Stokes-Cahn-Hilliard models do not allow for surface tension effects to be neglected due to the presence of the surface tension parameter in the Cahn-Hilliard equation. This motivates the proposed formulation, which allows surface tension effects to be changed without affecting the behaviour of the phase transport. Two methods are proposed: The first uses stabilised SUPG/PSPG linear elements, while the second is based on mixed Taylor-Hood elements. The proposed models are applied to a number of benchmark and example problems, including both capillary regime in which surface tension effects are dominant, and inertial regime in which surface tension effects are negligibly small. All results obtained agree very well with reference solutions.

Phase-field modelling, Navier-Stokes-Cahn-Hilliard equations, Surface tension.

§ INTRODUCTION

Modelling multiphase fluid flow in the presence of surface tension has for a long period been the subject of scientific investigation. These types of flows can be observed in numerous industrial and natural applications, for instance in emulsification, fluidized beds, combustion reactors, gas-liquid pipeline flows, etc. Typically multiphase problems are subject to a high degree of topological change, making interface tracking difficult, if not impossible. Strategies for simulating such flows include volume-of-fluid [1], level-set [2] and diffuse interface methods [3, 4, 5], as well as arbitrary Lagrangian-Eulerian [6, 7, 8, 9] or purely Lagrangian finite element formulations [10, 11, 12, 13, 14].

Diffuse-interface methods have become a popular choice for the simulation of multiphase problems involving complex topological changes. The fundamental concepts behind diffuse-interface methods for immiscible multiphase fluids were pioneered by Van der Waals [15] and Korteweg [16]. For an in depth explanation of the types of phase-field models, refer to [17, 18]. The fundamental principle is that instead of having a sharp interface partitioning two phases, the interface is modelled as a region within which physical quantities (density, viscosity etc.) and interfacial forces vary smoothly. In order to achieve this an auxiliary field variable, typically called the phase-field variable or order parameter, is introduced to smoothly transition quantities across the interface between two distinct values in the bulk phases. The dynamics of the phase-field variable are typically governed by the Cahn-Hilliard equation [19], or its non-conserving counterpart the Allen-Cahn equation [20]. Various numerical strategies for solving the Cahn-Hilliard and Allen-Cahn equations are explored in for instance [21, 22, 23, 24].

Coupling the Cahn-Hilliard equation with hydrodynamic Navier-Stokes equations results in a set of equations known as the Navier-Stokes-Cahn-Hilliard equations. These equations have been extensively used in the simulation of immiscible two-component fluid flows, and have been widely researched from a mathematical and computational perspective, see for instance [25, 26, 27, 28, 29]. These types of equations were utilised as early as the 1970's [30], with the so-called Model H. This model combines the phase dynamics of the Cahn-Hilliard equation with the hydrodynamics of the Navier-Stokes equation for the purpose of modelling immiscible, incompressible two-component flows. The drawback of the Model H is that it assumes that the mixture is comprised of constant matching densities, and is hence severely restricted in its application.
The difficulty with considering non-matching densities is that the macroscopic description of the density differs from the direct average of the microscopic descriptions, meaning that even if the individual components are incompressible, the mixture may not be. Typically models proposed for mixtures of non-matching densities are characterised as either: incompressible - where the volume averaged velocity field is divergence-free (i.e. satisfies incompressibility), and  quasi-incompressible - where the mass-averaged velocity field satisfies mass conservation only, resulting in a slight compressibility in the interface region.
A quasi-incompressible modification of the Model H by Lowengrub and Truskinovsky [28] allows for the the use of non-matching densities in the bulk regions. However, as mentioned the definition of the average velocity results in the loss of the divergence-free velocity field. Furthermore, a stronger coupling is present since the pressure from the momentum equation of the Navier-Stokes equations is present in the Cahn-Hilliard part of the equations.
More recently, another quasi-incompressible form of the equations is proposed in [31], where notably a linear method is used which satisfies discrete energy dissipation unconditionally.
On the other hand, models by Boyer [5] and Ding et al. [32] present appealing formulations with divergence-free velocity fields; however the energy inequalities are unknown for the model.
Shen et al. [29] propose a physically consistent model for incompressible mixtures with non-matching densities, implementing energy stable, accurate time integration schemes. More recently, Abels et al. [4] derived a variable density variation of the Model H which is thermodynamically consistent. The model is shown to recover sharp interfaces when the interface thickness tends to zero. The Abels et al. model has been effectively used in numerous numerical studies, and is chosen as the base model of this paper. The Allen-Cahn equation has also been successfully coupled with the Navier-Stokes equations in [33, 34, 35].

Numerous numerical strategies have been developed to solve the types of Navier-Stokes-Cahn-Hilliard equations discussed above.
The use of a spectral-element type discretisation has been shown effective in [36, 37, 38].
A finite difference discretisation is effectively demonstrated by Kim et al. [39].
Approaches based on finite element methods, including those with adaptive re-meshing, are also frequently employed (see, for instance [40, 41, 42, 43, 44]).
Guo et al. [45] present an adaptive mesh strategy, using an energy preserving $C_0$ finite element method to solve problems of high topological complexity. A discontinuous Galerkin finite element approach is used by Giesselmann and Pryer [46] to model the quasi-incompressible Navier-Stokes-Cahn-Hilliard equations.
One of the more recent approaches is based on using isogeometric analysis, which is successfully demonstrated in [18, 47, 48].
There have also been recent advancements in fractional step type schemes for the mentioned incompressible phase field models. Shen et al. [29] present a three stage decoupling where the computations of the Allen-Cahn equation, i.e. the phase-field variable and chemical potential, are decoupled from the Navier-Stokes equations, in which the pressure and velocity fields are also decoupled. Similar strategies for the Navier-Stokes-Chan-Hilliard equations is proposed in [49, 50, 51], and [52], where the latter only decouples the Cahn-Hilliard equation from the Navier-Stokes equations.

Similar numerical strategies have been demonstrated using the Navier-Stokes-Allen-Cahn equations. For instance Yang et al. [53] and more recently Chiu [54] successfully implemented a Navier-Stokes-Allen-Cahn type model to simulate surface tension dominated drop dynamics problems.

The primary objectives of this work are:

* To present a new Navier-Stokes-Cahn-Hilliard model capable of simulating problems  where the choice of the surface tension coefficient is based on the physics of the problem at hand.
This addresses the drawback of traditional forms of the Navier-Stokes-Cahn-Hilliard equations, wherein the surface tension coefficient cannot be set to zero without removing all stabilising effects of the surface and bulk energies.

* To present computationally efficient, low order finite element formulations of the Navier-Stokes-Cahn-Hilliard equations. As mentioned previously, the simulation of these equations generally involves the use of higher order b-splines or adaptive re-meshing,

which are generally associated with high computational cost. In this work emphasis is placed on computational efficiency, hence an equal order SUPG/PSPG stabilisation strategy as well as a mixed Taylor-Hood type methodology based on standard finite element discretisation with Lagrange polynomials are presented.

* To verify the presented methods with benchmark problems and realistic examples of two-phase flows. A focus is placed on reproducing a number of surface tension dominated benchmark problems. Additionally, sloshing and breaking dam problems are considered where surface tension effects are negligible and have therefore been deactivated during simulation.

The remainder of this work is structured as follows: In Section <ref>, the governing equations associated with the Navier-Stokes-Cahn-Hilliard equations are presented. In Section <ref> these equations are formulated for mixed and stabilised finite elements. In Section <ref> various numerical examples are demonstrated using the proposed formulations. In Section <ref> conclusions are drawn.

§ GOVERNING EQUATIONS

§.§ Cahn-Hilliard equation

Phase-field models are characterised by the introduction of an auxiliary function $\varphi$, a phase-field variable which localises the individual phases. $\varphi$ is represented by distinct values outside the interface region (i.e. the bulk region), for example,

\begin{align}
\varphi(x,t) =
\begin{cases}
1,  &\quad \textrm{phase $a$}\\
-1, &\quad \textrm{phase $b$} .
\end{cases}
\label{eq:order_paramater_def}
\end{align}

The Cahn-Hilliard equation governs the transport and decomposition of $\varphi$ through dissipation of the energy functional,

\begin{equation}
\mathcal{F}(\varphi) = \intom \Psi \dom  ,
\label{eq:GL_functional}
\end{equation}

\begin{align}
\Psi = \left( F(\varphi) + \frac{\epsilon^2}{2}\lvert \nabla \varphi\rvert^2 \right).
\end{align}

$\Psi$ is the Ginzburg-Landau free energy composed of a hydrophobic bulk component $F(\varphi)$, and a hydrophilic surface component $\frac{\epsilon^2}{2}\lvert \nabla \varphi\rvert^2$. Here $\epsilon$ is a length scale related to the interface thickness. In the context of this work we consider the bulk free energy component as a double well potential,

\begin{align}
F(\varphi) = W \left( 1 - \varphi^2 \right)^2,
\end{align}

where $W$ is the height of the well.

An in depth mathematical analysis and derivation of the Cahn-Hilliard equation can be found for instance in [55, 56, 57, 22].
The Cahn-Hilliard equation is presented in a two equation form, similarly to [23]: Consider a domain $\Omega \subset \mathbb{R} ^d$ ($d \leq 3$) with boundary $\Gamma$ which is separable into Dirichlet and Neumann subsets, $\Gamma_g$ and $\Gamma_q$. The solution variables $\varphi$ and $\eta$ are described by the following governing equations:

+ a ·∇φ- ∇·( M(φ) ∇η)
= 0           in  Ω_N

η- f(φ) + ϵ^2 Δφ   = 0           in  Ω_N

∇η·   = 0           on  Γ_q

∇φ·+ ‖∇φ‖cos(α)    = 0           on  Γ_q.

Here $\eta$ is the chemical potential, which is the variational derivative of the free energy given by Equation (<ref>), with $f(\varphi)=4W(\varphi^3-\varphi)$. $M(\varphi)$ is the mobility function, which is described in Section <ref>. The angle $\alpha$ is the three-phase contact angle at the boundary.

§.§ Navier-Stokes-Cahn-Hilliard equations

Consider again the domain $\Omega \subset \mathbb{R} ^d$ ($d \leq 3$) bounded by $\Gamma$ and containing a mixture of two immiscible incompressible fluids with different densities $\rho_a$ and $\rho_b$. In this work, the following NSCH model is used:

ρ(φ) (∂/∂t
+ (·∇)- b )
+ (J(φ,η)·∇)+ ∇p
- ∇·( 2μ(φ) ∇^s )
- κη∇φ   = 0           in  Ω

∇·   = 0           in  Ω

+ u ·∇φ- ∇·( M(φ) ∇η)
= 0           in  Ω

η- f(φ) + ϵ^2 Δφ   = 0           in  Ω

∇η·   = 0           on  Γ_q

∇φ·+ ‖∇φ‖cos(α)    = 0           on  Γ_q

where $\vu$ and $p$ represent the velocity and pressure, respectively, and $\bm{b}$ is the body force. The symmetric gradient operator is defined as $\nabla^s := \frac{1}{2}\left(\nabla + \nabla^T\right)$.

The model presented in Equation (<ref>) is based on the Abels et al. [4] and is thermodynamically consistent. It allows for variable densities, and agrees with sharp interface models for $\epsilon \to 0$.

Cahn and Hilliard [19] describe surface tension as the excess free energy per unit surface area. It follows that in the case of a plane phase-field interface at equilibrium, the surface tension coefficient, $\gamma$, is related to the phase-field variable $\varphi$ by

\begin{align}
\gamma = \kappa \int^{\infty}_{-\infty} \left( \sdr{}{\varphi}{x} \right) \mathrm{d} x.
\label{eq:gamma}
\end{align}

In order to achieve consistency of the surface tension term in (<ref>) with the Young-Laplace sharp interface surface tension model, the factor $\kappa$ must be chosen as

\begin{align}
\kappa = \frac{3}{4\sqrt{2W}} \frac{\gamma}{\epsilon},
\label{eq:kappa}
\end{align}

Equation (<ref>) is successfully implemented in numerous studies [27, 58, 32, 48]. An NSCH model using Equation (<ref>) is shown to approach the sharp interface model as $\epsilon \to 0$ in [59].

In the present work, the surface tension stress term $\eta\nabla\varphi$ (see also Boyer [5]) is used in Equation (<ref>), instead of $\gamma\epsilon\nabla\cdot(\nabla\varphi\otimes \nabla\varphi)$ (see Abels et al. [4]). The former term is easily derived from the latter, and the resulting pressure term in Equation (<ref>) is now a modified pressure term $\nabla \hat{p}$, where $\hat{p}=(p+\Psi)$. This modified pressure is similar to the original pressure, in fact, it is identical in the bulk regions where $\Psi$ vanishes. Advantageously the modified pressure varies generally more smoothly across interface regions than the original pressure.

In Equation (<ref>), the relative diffusive flux is expressed as

\begin{align}
\bm{J} = -\rho_{dif} M(\varphi) \nabla \eta .
\end{align}

The density and viscosity are described by the linear approximations

\begin{align}
\begin{alignedat}{2}
\rho(\varphi) &= \rho_{dif}\varphi &&+ \rho_{avg}
\\\mu(\varphi) &= \mu_{dif}\varphi &&+  \mu_{avg},
\end{alignedat}
\label{eq:density viscosity interp}
\end{align}

\begin{align}
\rho_{dif} = \frac{\rho_a-\rho_b}{2}, \quad \rho_{avg} = \frac{\rho_a+\rho_b}{2} \quad \textrm{and} \quad
\mu_{dif} = \frac{\mu_a-\mu_b}{2}, \quad \mu_{avg} = \frac{\mu_a+\mu_b}{2}.
\label{eq:density viscosity interp2}
\end{align}

If necessary, a cut-off function $\bar{\varphi}$ can be incorporated in Equation (<ref>), where

\begin{align}
\bar{\varphi} =
\begin{cases}
\varphi \quad&\text{if} \lvert\varphi\rvert \leq 1 \\
\mathrm{sign}(\varphi) \quad&\text{if} \lvert\varphi\rvert > 1.
\end{cases}
\label{eq:cut-off}
\end{align}

This would ensure that $\rho$ and $\mu$ remain within the physical bounds of the specified bulk phases, and has been effectively used in [60].

Remark 1.

In the model given by Equation (<ref>) the physical quantities and the Young-Laplace surface tension term occur only in the momentum equation. The CH Equation (<ref>) merely governs the transport of the order parameter $\varphi$. Traditional representations of the NSCH equations such as those used in [28, 32, 4] can easily be recovered by introducing appropriate scalar factors in Equations (<ref>) or (<ref>).

§.§ Mobility function

In this work three options are considered for the mobility, i.e.

\begin{align}
M_0 = D,
\end{align}

\begin{align}
M_2(\varphi) =
\begin{cases}
D(1 - \varphi^2), \quad &\text{if} \ \lvert \varphi \rvert \leq 1 \\
0, \quad &\text{elsewhere},
\end{cases}
\end{align}

\begin{align}
M_3(\varphi) =
\begin{cases}
D(-2\varphi^3 - 3\varphi^2 + 1) \quad&\text{if} \ \varphi \geq -1 \\
D(2\varphi^3 - 3\varphi^2 + 1) \quad&\text{if} \ \varphi \leq  1 \\
0 \quad&\text{elsewhere},
\end{cases}
\end{align}

where $D$ is a constant. The latter two mobility functions are degenerate, i.e. they are non-zero only in the interface region. This is illustrated in Figure <ref> for $D=1$.

Degenerate mobility functions, $M_{2,3}$ (left) and the respective first derivatives (right), with $D=1$ in both cases.

Notably, $M_3(\varphi)$, transitions smoother in the derivative.

On coarse meshes, the choice of mobility can drastically change the behaviour of the model. Choosing a large constant mobility can result in the acceleration of the Ostwald ripening or coarsening effect. In essence this means the total interfacial area will reduce with time in an effort to reach the lowest energy state and hence thermodynamic equilibrium. For this reason $D$ must be chosen small enough to ensure that the associated timescale is far larger than the time domain of interest. The reader is referred to [61] for a detailed explanation of Ostwald ripening.

§ NUMERICAL FORMULATION

§.§ Navier-Stokes-Cahn-Hilliard mixed Taylor-Hood formulation

In order to satisfy the LBB stability condition, a Taylor-Hood element is chosen for the spatial discretisation (see Figure <ref>). Thus, the velocity, phase-field variable and chemical potential interpolations are piecewise quadratic, while the pressure interpolation is piecewise linear.

[] (0,0) – (1,0) – (1,1) – (0,1) – cycle;
[black,fill] (0,0) circle(0.04);
[black,fill] (0.5,0) circle(0.04);
[black,fill] (1,0) circle(0.04);
[black,fill] (1,0.5) circle(0.04);
[black,fill] (1,1) circle(0.04);
[black,fill] (0.5,1) circle(0.04);
[black,fill] (0,1) circle(0.04);
[black,fill] (0,0.5) circle(0.04);
[black,fill] (0.5,0.5) circle(0.04);
[black] (0,0) circle(0.07);
[black] (1,0) circle(0.07);
[black] (1,1) circle(0.07);
[black] (0,1) circle(0.07);

[white] (0,0) – (1,0) – (1,1) – (0,1) – cycle;
[white,fill] (0,0) circle(0.05);
[white,fill] (0.5,0) circle(0.05);
[white,fill] (1,0) circle(0.05);
[white,fill] (1,0.5) circle(0.05);
[white,fill] (1,1) circle(0.05);
[white,fill] (0.5,1) circle(0.05);
[white,fill] (0,1) circle(0.05);
[white,fill] (0,0.5) circle(0.05);
[white,fill] (0.5,0.5) circle(0.05);
[white] (0,0) circle(0.1);
[white] (1,0) circle(0.1);
[white] (1,1) circle(0.1);
[white] (0,1) circle(0.1);
[] at (0.05,0.6) (A) Nodes:;
[black,fill] (0.45,0.6) circle(0.04);
[] at (1.11,0.6) (A) $\vu$, $\varphi$, $\eta$ (quadratic);
[black] (0.45,0.4) circle(0.05);
[] at (0.90,0.4) (A) $p$ (linear);

Two-dimensional Taylor-Hood interpolations.

The weak form of the NSCH problem described in Equation (<ref>) reads: Find $(\vu^h$, $p^h$, $\varphi^h$, $\eta^h) \in (\mathcal{S}^h,\mathcal{P}^h,\mathcal{X}^h,\mathcal{Z}^h)$, such that, for all $(\vw^h$, $q^h$, $s^h$, $v^h) \in (\mathcal{V}^h,\mathcal{P}^h,\mathcal{Y}^h,\mathcal{W}^h)$,

\begin{align}
% momentum
\begin{split}
\int_{\Omega}
\vw^h &\cdot \left( \rho(\varphi^h) \left(  \frac{\partial\vu^h}{\partial t}
+ \vu^h \cdot \nabla \vu^h - \bm{b} \right) - \bm{J}(\varphi^h,\eta^h) \cdot \nabla \vu^h - \kappa\eta^h \nabla \varphi^h \right)
-\left(\nabla \cdot \vw^h \right) p^h
\\&+\nabla \vw^h : \left( 2 \mu(\varphi^h) \nabla^s \vu^h \right) \dom = 0
\end{split}
\label{eq:NSCHweakA}
\end{align}

\begin{equation}
% continuity
\int_{\Omega} \nabla q^h \cdot \vu^h \dom = 0
\label{eq:NSCHweakB}
\end{equation}

\begin{equation}
% Cahn-Hilliard
\intom v^h \left( \frac{\partial \varphi^h}{\partial t} + \vu^h\cdot \nabla \varphi^h \right)
+ \nabla v^h \cdot \left( M(\varphi^h) \nabla \eta^h \right) \dom = 0
\label{eq:NSCHweakC}
\end{equation}

\begin{equation}
% auxiliary
\intom s^h \ \left( \eta^h - f(\varphi^h) \right) - \nabla s^h \cdot \epsilon^2 \nabla \varphi^h \dom
- \intgam s^h \ \epsilon^2 \lVert\nabla\varphi^h\rVert \cos{(\alpha)} \dgam = 0,
\label{eq:NSCHweakD}
\end{equation}

where $\mathcal{S}^h$, $\mathcal{P}^h$, $\mathcal{X}^h$, $\mathcal{Z}^h$, $\mathcal{V}^h$, $\mathcal{Y}^h$ and $\mathcal{W}^h$ are the appropriate finite element spaces of piecewise continuous quadratic and linear basis functions.

§.§ Navier-Stokes-Cahn-Hilliard stabilised formulation

The standard SUPG/PSPG stabilisation strategy (see, for instance [62, 63, 64, 65, 66] ) is applied to the Navier-Stokes momentum equation. This allows for the use of piecewise linear equal order interpolations for $u$, $p$, $\varphi$ and $\eta$.

The SUPG/PSPG stabilised weak form of (<ref>) reads: Find $(\vu^h$, $p^h$, $\varphi^h$, $\eta^h) \in (\mathcal{S}^h,\mathcal{P}^h,\mathcal{X}^h,\mathcal{Z}^h)$, such that, for all $(\vw^h$, $q^h$, $s^h$, $v^h) \in (\mathcal{V}^h,\mathcal{P}^h,\mathcal{Y}^h,\mathcal{W}^h)$,

\begin{align}
% momentum
\begin{split}
\vw^h \cdot \left( \rho(\varphi^h)\left(  \frac{\partial\vu^h}{\partial t}
+ \vu^h \cdot \nabla \vu^h - \bm{b} \right)
- \bm{J}(\varphi^h,\eta^h) \cdot \nabla \vu^h - \kappa\eta^h \nabla \varphi^h \right)
-\left(\nabla \cdot \vw^h\right) p^h
\\&\quad\qquad\qquad+\nabla \vw^h : \left( 2 \mu(\varphi^h) \nabla^s \vu^h \right) \dom
\\&+\sum_{e=1}^{n_{el}} \intome \Big[ \tau_{\vu} \rho(\varphi^h) \left(\vu^h \cdot \nabla \vw^h \right) + \tau_p \nabla q^h \Big]
\\&\quad\qquad\qquad\cdot \Big[ \rho(\varphi^h)\left( \frac{\partial\vu^h}{\partial t} + \vu^h \cdot \nabla \vu^h - \bm{b} \right) - \bm{J}(\varphi^h,\eta^h) \cdot \nabla \vu^h - \kappa\eta^h \nabla \varphi^h + \nabla p^h \Big] \dom = 0
\end{split}
\label{eq:NSCHstabWeakA}
\end{align}

\begin{equation}
% continuity
\int_{\Omega} \nabla q^h \cdot \vu^h \dom = 0
\label{eq:NSCHstabWeakB}
\end{equation}

\begin{equation}
% Cahn-Hilliard
\intom v^h \left( \frac{\partial \varphi^h}{\partial t} + \vu^h\cdot \nabla \varphi^h \right)
+ \nabla v^h \cdot \left(M(\varphi^h)\nabla\eta^h\right) \dom = 0
\label{eq:NSCHstabWeakC}
\end{equation}

\begin{equation}
% Chemical potential
\intom s^h \left( \eta^h - f \right) - \nabla s^h \cdot \epsilon^2 \nabla \varphi^h \dom
- \intgam s^h \ \epsilon^2 \lVert\nabla\varphi^h\rVert \cos{(\alpha)} \dgam= 0.
\label{eq:NSCHstabWeakD}
\end{equation}

where $\mathcal{S}^h$, $\mathcal{P}^h$, $\mathcal{X}^h$, $\mathcal{Z}^h$, $\mathcal{V}^h$, $\mathcal{Y}^h$ and $\mathcal{W}^h$ are the appropriate finite element spaces of piecewise continuous linear basis functions.

Following [67], the stabilisation parameters, $\tau_{\vu}$ and $\tau_p$ are defined as

\begin{align}
\tau_{\vu} = \left( \frac{1}{\tau_p^2} + \left(\frac{2 \rho \lVert \vu_e \rVert}{h_e}\right)^2 \right)^{-\frac{1}{2}},
\quad
\tau_{p} = \frac{h_e^2}{4\mu},
\label{eq:stab parameters}
\end{align}

where $h_e$ is the characteristic size of the element, evaluated as $h_e=V_e^{1/d}$, with $d$ representing the number of spatial dimensions and $V_e$, the element volume or area. The vector $\vu^e$ is the velocity in the element centroid.

§.§ Temporal discretisation

The generalised-$\alpha$ method is employed for the temporal discretisation.
This method is an unconditionally stable, implicit single-step time integration scheme (refer to [68, 69]). The scheme allows for high frequency damping to be controlled by without compromising the second order accuracy (see [64, 70]). Applying the generalised-$\alpha$ method to a generic first order problem gives

\begin{equation}
\dvu^{n + \alpha_m} = \bm{f} \left(t^{n + \alpha_f},\bm{u}^{n + \alpha_f} \right),
\label{eq:AMgeneral}
\end{equation}

\begin{align}
\dvu^{n+\alpha_m} &= \left(1 - \alpha_m \right) \dvu^n + \alpha_m \dvu^{n+1}
\label{eq:AMacc}
\\t^{n+\alpha_f} &= \left(1 - \alpha_f \right)t^n + \alpha_f t^{n+1}
\label{eq:AMtime}
\\\vu^{n+\alpha_f} &= \left(1 - \alpha_f \right)\vu^n + \alpha_f \vu^{n+1}
\label{eq:AMvelf}
\\\frac{\vu^{n+1} - \vu^{n}}{\dt} &= \left(1 - \gamma \right) \dvu^n + \gamma \dvu^{n+1},
\label{eq:AMvel}
\end{align}

\begin{equation}
\alpha_m = \frac{1}{2} \ \frac{3 - \rho_\infty}{1 + \rho_\infty}
\text{,}\qquad
\alpha_f = \frac{1}{1 + \rho_\infty}
\text{,}\qquad
\gamma = \frac{1}{2} + \alpha_m - \alpha_f,
\label{eq:AMpar}
\end{equation}

where $\gamma$ must not be confused with the surface tension parameter. Note that $\alpha_m$, $\alpha_f$ and $\gamma$ are expressed in terms of the spectral radius $\rhoi$ for an infinitely large time step size.
This scheme has been successfully employed in solution of a number of challenging coupled field problems [9, 67, 71, 72].

§.§ Solver

The problems described by Equations (<ref>) and (<ref>) are highly nonlinear and thus a Newton-Raphson procedure is employed based on the consistent linearisation of all nonlinear terms. For smaller scale problems, a direct linear solver (PARDISO [73, 74, 75]) is used, whereas for larger scale problems, an iterative parallel solver (PETSC [76]) with a block Jacobi pre-conditioner is employed.

§.§ Adaptive time stepping

Several of the numerical examples in Section <ref> implement an adaptive time stepping procedure defined by

\begin{align}
\Delta t = \left( t^{n+1}-t^{n} \right)\theta^{\left(n_{\mathrm{iter}}-n_{\mathrm{opt}}\right)},
\end{align}

where $n_{\mathrm{iter}}$ and $n_{\mathrm{opt}}$ are respectively the number of Newton iteration steps required in the previous time step and the desired number of iteration steps. The constant $\theta$ is typically chosen between $0.5$ and $1$.

§ NUMERICAL EXAMPLES

In the following examples the proposed methodology is applied to a number of two and three-dimensional problems. Sections <ref>-<ref> present surface tension dominated problems, while the simulations described in Sections <ref>-<ref> do not feature any surface tension effects.

The properties of water and air used in the following examples are:

water:      ρ= 0.998        [per-mode=symbol]^3

μ= 0.0101        [per-mode=symbol]

air:      ρ= 0.0012        [per-mode=symbol]^3

μ= 0.000182        [per-mode=symbol]

water/air:      γ= 73.0        [per-mode=symbol]^2

gravity:      g = 980.0        [per-mode=symbol]^2

Unless otherwise stated, the computations are based on the stabilised formulation given by Equation (<ref>) and the Cahn-Hilliard parameters are set as follows:

mobility function:      M_0       = 10^-3  [per-mode=symbol]^2

well height:      W       = 0.25

interface thickness:      ϵ      = 2h

where $h$ corresponds to the characteristic size of the element used, which is taken as the largest element size of the respective mesh.
The relative diffusive flux term $\left(\bm{J}\cdot\nabla\right) \vu$ is neglected in the following examples as its effect is observably negligible.

§.§ Static bubble

A bubble of radius $R$ is placed at coordinate $[0.5,0.5]$ in the centre of a square domain of dimensions $[0,1]\times[0,1]$. The densities and viscosities, as well as the surface tension coefficient, are set to $1$. The initial condition is taken as

\begin{align}
\varphi(\bm{x}) = \tanh\left( \frac{R - d(\bm{x})}{\sqrt{2}\epsilon} \right).
\end{align}

where $d(\bm{x})$ is the Euclidean distance between the bubble centre and $\bm{x}$.

The pressure difference between the centre of the bubble and a point outside the bubble at coordinate $[1.0,0.5]$, are compared to the Young-Laplace equation,

\begin{align}
\Delta p = \frac{\gamma}{R}.
\label{eq: Young-Laplace}
\end{align}

Two meshes are considered: a mesh with $256\times256$ mixed Taylor-Hood elements, and a mesh with $512\times512$ stabilised linear quadrilateral elements. For both meshes, the interface thickness is set as $\epsilon = 1/128$. All simulations are run until $T=10$.
Figure <ref> shows the excellent agreement between the analytical expression and the numerical experiment, for a range of values of $R$.

Static bubble: Comparison between numerical and analytical (Young-Laplace) pressure drop.

§.§ Small amplitude oscillations of a two-dimensional drop

A small two-dimensional drop of water with radius $R=0.0125 \ \si{cm}$ is placed in a square domain filled with air. The initial configuration of the drop is described by the equation

\begin{align}
r_n(\theta) = R + A \cos(n\theta),
\end{align}

where $A$ is set to $0.02 R$ for all simulations, and $n$ is the mode order (see [77, 9, 13, 78]). The effects of gravity are neglected, and the mobility is set as $M_3(\varphi)$, with $D=10^{-2}$. The mesh considered has $512\times512$ linear quadrilateral elements, and a fixed time step size of $\Delta t =1\times10^{-6} \ \si{\second}$ is used.
Figure <ref> shows the evolution of the oscillation amplitude for modes $n=2,3,4$. The observed period matches very well with the analytical period,

\begin{align}
\tau = \frac{2\pi}{\omega_n},
\label{eq: oscillations analytical period}
\end{align}

where following [13],

\begin{align}
\omega_n^2 = \left( n^3 - n \right) \frac{\gamma}{\rho R^3},
\end{align}

Small amplitude oscillations of a two-dimensional drop: Amplitude evolution for $n=2,3,4$, with $512\times512$ linear elements.

A comparison of the numerical and analytical periods is given in Table <ref>.

Small amplitude oscillations of a two-dimensional drop: Numerical and analytical periods.

$n$     $\tau$ [numerical]     $\tau$ [analytical (<ref>)]

$2$     $4.190 \cdot 10^{-4}$     $4.192 \cdot 10^{-4}$

$3$     $2.110 \cdot 10^{-4}$     $2.096 \cdot 10^{-4}$

$4$     $1.350 \cdot 10^{-4}$     $1.326 \cdot 10^{-4}$

§.§ Large amplitude oscillation of a two-dimensional drop

A two-dimensional drop is set up similarly to Section <ref>, with a larger initial amplitude $A=0.2R$.
Figure <ref> shows the evolution of the amplitude for the cases $n=2,3,4$. It is clear that the periods observed in the present study slightly lag behind the analytical period. This is due to the larger amplitudes (compared to Section <ref>), which exceed the linear range and cause strongly nonlinear behaviour.

Large amplitude oscillation of a two-dimensional drop: Amplitude evolution for $n=2,3,4$, with $512\times512$ linear elements.

Figure <ref> shows the interface in terms of the isolines of $\varphi$ at $\varphi=0$.

$t=0 \ \si{\second}$

$t=1.08\cdot 10^{-4} \ \si{\second}$

$t=2.2\cdot 10^{-4} \ \si{\second}$

$t=3.332\cdot 10^{-4} \ \si{\second}$

$t=4.4\cdot 10^{-4} \ \si{\second}$

$t=0 \ \si{\second}$

$t=6\cdot 10^{-5} \ \si{\second}$

$t=1.12\cdot 10^{-4} \ \si{\second}$

$t=1.72\cdot 10^{-4} \ \si{\second}$

$t=2.2\cdot 10^{-4} \ \si{\second}$

$t=0 \ \si{\second}$

$t=4\cdot 10^{-5} \ \si{\second}$

$t=7.6\cdot 10^{-5} \ \si{\second}$

$t=1.16\cdot 10^{-4} \ \si{\second}$

$t=1.52\cdot 10^{-4} \ \si{\second}$

Large amplitude oscillation of a two-dimensional drop: Interface isolines for $n=2,3,4$.

Figure <ref> shows the velocity streamlines at different time instances. Notably the axes of symmetry for all three modes are clearly defined.

$t=7.2\cdot 10^{-5} \ \si{\second}$

$t=4\cdot 10^{-5} \ \si{\second}$

$t=3.2\cdot 10^{-5} \ \si{\second}$

Large amplitude oscillations of a two-dimensional drop: Velocity streamlines at different time instances for $n=2,3,4$

Figure <ref> shows the pressure distributions at different times within a single period of oscillation.
Clearly the pressure variation in the air is minimal, and higher pressure concentrations are present in regions with more curvature, as expected from the Young-Laplace equation.

$t=4\cdot 10^{-6} \ \si{\second}$

$t=4\cdot 10^{-6} \ \si{\second}$

$t=4\cdot 10^{-6} \ \si{\second}$

$t=1.08\cdot 10^{-4} \ \si{\second}$

$t=6\cdot 10^{-5} \ \si{\second}$

$t=4\cdot 10^{-5} \ \si{\second}$

$t=2.2\cdot 10^{-4} \ \si{\second}$

$t=1.12\cdot 10^{-4} \ \si{\second}$

$t=7.6\cdot 10^{-5} \ \si{\second}$

$t=3.332\cdot 10^{-4} \ \si{\second}$

$t=1.72\cdot 10^{-4} \ \si{\second}$

$t=1.16\cdot 10^{-4} \ \si{\second}$

$t=4.4\cdot 10^{-4} \ \si{\second}$

$t=2.2\cdot 10^{-4} \ \si{\second}$

$t=1.52\cdot 10^{-4} \ \si{\second}$

Large amplitude oscillation of a two-dimensional drop: Pressure contour plots.

The accurate conservation of the volumes of water and air are demonstrated in Figure <ref> for $n=4$, where the relative volume error is defined as $\epsilon(t_n) = (V_n - V_0)/V_0$. It is observed that for both air and water, the relative error does not exceed $0.015\%$.

Large amplitude oscillation of a two-dimensional drop: Relative volume changes for air and water for $n=4$.

§.§ Capillary rise

The capillary rise of fluid between two parallel plates is simulated. Figure <ref> illustrates typical initial and final configurations of this problem.

Capillary rise: Initial (left) and final (right) configurations.

The column height can accurately be approximated using Jurin's law [79]:

\begin{align}
\Delta h &= \frac{\gamma \cos{\alpha}}{\rho g R} \quad \textrm{(2D)}
\label{eq:JurinsLawa}
\\\Delta h &= \frac{2\gamma \cos{\alpha}}{\rho g R} \quad \textrm{(axisymmetric)},
\label{eq:JurinsLawb}
\end{align}

where $R$ is the pipe radius and $\alpha$ is the contact angle.
The initial configuration of the two-dimensional problem is as follows: The initial water level in a $[0,10]\times[0,8] \ \si{\milli\metre^2}$ container filled with air is set to $4 \ \si{\milli\metre}$. Two parallel plates at a distance $R=1 \ \si{\milli\metre}$, are placed in the centre of the container. Two meshes are considered: a coarse mesh with $64\times24$ linear quadrilateral elements and a fine mesh with $256\times96$ linear quadrilateral elements. Slip boundary conditions are applied on the side and top boundaries, as well as on the tube walls. The contact angle $\alpha$ is applied on the inner tube walls. The external air pressure is set to zero. Adaptive time stepping is used, with the initial time step size being $\Delta t = 0.001$.

The fluid levels $\Delta h$ are shown with respect to various hydrophilic and hydrophobic contact angles in Figure <ref>. The numerical fluid levels are obtained as

\begin{align}
\Delta h = \frac{1}{3}\left[ \Delta h(0) + \Delta h(R/2) + \Delta h(R) \right].
\end{align}

Capillary rise between two parallel plates: Fluid level evolution with respect to time (left) and contact angle (right) compared to Jurin's law, with $256\times96$ linear elements.

The fluid level heights are tabulated in Table <ref>.

Capillary rise between two parallel plates: Fluid level heights.

$\alpha$     $\Delta h$ [coarse]     $\Delta h$ [fine]     $\Delta h$ [Jurin's law (<ref>)]

$70 ^{\circ}$     $0.2587$      $0.2567$      $0.2553$

$80 ^{\circ}$     $0.1295$      $0.1292$     $0.1296$

$100^{\circ}$     $-0.1294$     $-0.1292$     $-0.1296$

$110^{\circ}$     $-0.2583$     $-0.2567$     $-0.2553$

The equilibrium configurations are shown in Figure <ref>.

Capillary rise between two parallel plates: One half of the equilibrium configurations of $\varphi$ for various contact angles, with $256\times96$ linear elements.

Capillary rise in a circular pipe is also considered. Following Equation (<ref>) it is expected that twice the fluid level height is reached. The geometry and mesh are shown in Figure <ref>.

Capillary rise in a circular pipe: Slice through horizontal axis of finite element mesh with 1,559,048 linear hexahedron elements.

The steady state solutions to $\varphi$ are shown in Figure <ref> and

Capillary rise in a circular pipe: Equilibrium configurations of $\varphi$, the part of the domain filled with air is not shown.

Table <ref> shows corresponding fluid levels for the present numerical solution and the analytical solution.

Capillary rise in a circular pipe: Fluid level heights.

$\alpha$     $\Delta h$     $\Delta h$ [Jurin's law (<ref>)]

$80 ^{\circ}$     $0.2574$      $0.2595$

$100^{\circ}$     $-0.2557 $     $-0.2595$

Clearly, excellent correspondence between numerical and analytical results is observed in all comparisons.

§.§ Sessile drop

A small water bubble in the presence of gravity is placed on a solid surface at the bottom of a rectangular domain filled with air, as shown in Figure <ref>.

Sessile drops: Geometry.

The rectangular domain is of dimensions $[0,1.5]\times[0,0.5] \ \si{cm}^2$ and the drop radius is $R=0.25 \ \si{\centi\metre}$. Both two and three-dimensional sessile drops are considered with the same radius. In order to reduce dynamic effects and accelerate the computation of the equilibrium configuration, the viscosities are multiplied by a factor of $100$. Slip boundary conditions are applied at the lower boundary.
A mesh with $96\times64$ linear quadrilateral elements is considered. Adaptive time stepping is used, with an initial time step size of $\Delta t = 0.001$.

The final configurations of the two-dimensional drops subject to various contact angles are shown in Figure <ref>.

Sessile drop in two dimensions: Equilibrium configuration of $\varphi$ using $96\times64$ linear elements.

In Figure <ref>, the results are compared to the exact equilibrium configuration obtained by numerically solving a set of ordinary differential equations [80].

Sessile drop in two dimensions: Comparison of numerical results with the analytical solution given by Pozrikidis [80], using $96\times64$ linear elements.

Figures <ref> and <ref> show the results obtained for the three-dimensional case. A very good match is observable for both two and three-dimensional cases.

Sessile drop: Steady state $\varphi$ solution for a small three-dimensional sessile drops with various contact angles, using $48\times32\times48$ linear elements.

Sessile drop: Comparison of a small three-dimensional sessile drops with various contact angles to the analytical solution given by Pozrikidis [80], using $48\times32\times48$ linear elements.

§.§ Rising bubble

In this example a bubble of radius $R=0.25$ is placed near the bottom of a rectangular domain of size $[0,1]\times[0,2]$ filled with a heavier fluid. The parameters are set as given in [48], with the heavier fluid properties $\rho = 1000$, $\mu = 10$, and the lighter bubble properties $\rho=1$, $\mu=1$. The surface tension parameter and the acceleration vector are set, respectively, to $\gamma = 1.96$ and $b = [0,-0.98]^T$. The mesh consists of $128\times256$ linear elements, and a fixed time step size of $\Delta t = 0.01$ is used. The boundaries are set up such that the fluids are allowed to slip along the left and right boundaries, but not along the top and bottom boundaries.
It is expected that due to the large density difference, there would be an overpowering body force contribution acting on either side of the interface in comparison to the rather negligible surface tension forces. Thus the bubble would experience large deformations characterised by long and narrow filaments at the edges. This behaviour is observable in Figure <ref>, where the solution to $\varphi$ is shown at different time instances. The observations match well with those made in [48].

Rising bubble: $\varphi$ evolution using $128\times256$ linear elements.

§.§ Faucet leak in two dimensions

In this example we consider a leaky faucet problem, which consists of water dripping from a faucet into a pool of water at the bottom of a rectangular box filled with air, as shown in Figure <ref>. The dimensions of the box are $[0,2] \times [0,4.5] \ \si{\centi\metre^2}$, the height of the pool is $H=0.4 \ \si{\centi\metre}$, and the radius of the faucet is $R=0.26 \ \si{\centi\metre}$. The simulations are run with a filling velocity of $\vu_{\mathrm{fill}} = 10 \ \si{\centi\metre\per\second}$, which is set uniformly at the inlet.

At the beginning of the simulation the faucet is completely closed, after which the inlet velocity is ramped smoothly up to the the filling velocity with the function $\frac{1}{2} \left( 1 - \cos(\pi t/T) \right)$, where $T=0.01 \ \si{\second}$. A mesh comprising of 35,712 linear elements is considered in all simulations. Adaptive time stepping is used, with an initial time step size of $\Delta t = 0.0001$. The mobility is set as $M=0.5 \ \si{\centi\metre\per\second\squared}$. Slip boundary conditions are applied on all surfaces.

Filling drop in two dimensions: Geometry.

The evolution of $\varphi$ is shown in Figure <ref>, illustrating the ability of the described formulation to deal with highly complex physical problems resulting in large geometry modifications including significant topology changes.

$t=0 \ \si{\second}$

$t=0.0354 \ \si{\second}$

$t=0.0816 \ \si{\second}$

$t=0.1198 \ \si{\second}$

$t=0.1444 \ \si{\second}$

$t=0.1641 \ \si{\second}$

$t=0.1708 \ \si{\second}$

$t=0.1771 \ \si{\second}$

$t=0.1940 \ \si{\second}$

$t=0.2118 \ \si{\second}$

$t=0.2204 \ \si{\second}$

$t=0.2442 \ \si{\second}$

$t=0.2549 \ \si{\second}$

$t=0.2837 \ \si{\second}$

$t=0.2974 \ \si{\second}$

$t=0.3242 \ \si{\second}$

$t=0.3400 \ \si{\second}$

$t=0.3850 \ \si{\second}$

Faucet leak in two dimensions: Evolution of $\varphi$, using $35,712$ linear elements.

It is noted that small filling velocities lead to the development of non-physical flows in the air domain that originate from the tip of the drop. This phenomenon is likely to be triggered by the combination of the massively different values of density and viscosity of water and air with an imbalance of the selected mobility, interface thickness, mesh density and relatively large time domain of interest. It is currently further investigated by the authors.

§.§ Faucet leak in three dimensions

A three-dimensional faucet leak is now considered with a similar setup to Section <ref>. In this analysis the drop break-off is particularly of interest, thus the domain size is reduced to $[0,2 \ \si{\centi\metre}] \times [0,3 \ \si{\centi\metre}]$, and the pool at the bottom has been removed. A $15^\circ$ wedge is considered due to the symmetry of the problem, as shown in Figure <ref>, which consists of 1,004,665 linear tetrahedron elements. Adaptive time stepping is used, with an initial time step size of $\Delta t = 0.0001$. The filling velocity is set as $\vu_{\mathrm{fill}} = 10 \ \si{\centi\metre\per\second}$.

The evolution of $\varphi$ is shown in Figure <ref>.

Faucet leak in three dimensions: Geometry with $15^{\circ}$ wedge.

$0 \ \si{\second}$

$0.0363 \ \si{\second}$

$0.0550 \ \si{\second}$

$0.0869 \ \si{\second}$

$0.1098 \ \si{\second}$

$0.1241 \ \si{\second}$

$0.1359 \ \si{\second}$

$0.1453 \ \si{\second}$

$0.1467 \ \si{\second}$

$0.1478 \ \si{\second}$

$0.1499 \ \si{\second}$

Faucet leak in three dimensions: Evolution of $\varphi$, using $1,004,665$ linear elements.

A visual comparison of the drop configuration at different time instances is made with Dettmer and Perić [9] in Figure <ref>.

Faucet leak in three dimensions: Comparison of interface isolines with [9].

The vertical velocity component during break-off is visible in Figure <ref>.

$0.1438 \ \si{\second}$

$0.1463 \ \si{\second}$

$0.1468 \ \si{\second}$

$0.1486 \ \si{\second}$

Faucet leak in three dimensions: Evolution of vertical velocity, using $1,004,665$ linear elements.

In Figure <ref>, the tip position evolution compared with [9] is shown along with the water the volume evolution. Unlike in the present study, the initial configuration in [9] is set as an equilibrium pendant drop, thus these results are shifted in time to the appropriate position. It is observable that the tip position evolution matches well with [9].

Faucet leak in three dimensions: Tip position comparison with [9] (left), and water volume evolution comparison with analytical solution (right).

§.§ Broken dam

A fluid column ($\rho = 1$, $\mu=0.01$) of width $b = 3.5$ and height $h = 7$ is placed in a rectangular domain ($\rho = 0.001$, $\mu=0.0001$) of dimensions $[0,15] \times [0,10]$. The gravitational acceleration is set by using $\bm{b} = [0,-1]^T$, and the surface tension coefficient is set to zero, i.e. $\gamma=0$. The fluid is allowed to slip on the horizontal and vertical surfaces, and is allowed to adopt any contact angle by replacing the boundary terms in Equations (<ref>) and (<ref>), such that the terms become

\begin{align}
% 	\intgam  s \ \epsilon^2 \nabla \varphi \cdot \vn \dgam
\intom s^h \left( \eta^h - f \right) - \nabla s^h \cdot \epsilon^2 \nabla \varphi^h \dom
- \intgam s^h \ \epsilon^2 \nabla \varphi^h \cdot \vn \dgam = 0.
\label{eq:free boundary}
\end{align}

The mesh consists of $192\times128$ linear elements, and the time step size is set to $\Delta t = 0.01$. The geometry of the problem and the tip displacement evolution are shown in Figure <ref>.

Broken dam: Geometry (left) and evolution of tip displacement (right).

The tip displacement evolutions agrees well with experimental results taken from [11].

The evolution of $\varphi$ is shown at several time instances in Figure <ref>.

broken dam: $\varphi$ evolution using $192\times128$ elements.

§.§ Sloshing tank

Following [9], a sloshing tank of size $[0,1] \times [0,1.2]$ is set up with the following parameters $\rho_1 = 1$, $\rho_2 = 0.001$, $\mu_1 = 0.01$, $\mu_2 = 0.0001$ , $\bm{b} = [0,-1]^T$. The surface tension effects are neglected. A uniform mesh with 48,000 linear quadrilateral elements is selected, and a time step size of $\Delta t = 0.1$ is chosen. The fluid at the edges is allowed to adopt any contact angle by using (<ref>). Figure <ref> shows that the frequency at the left and right edges agrees excellently with [9].

Sloshing tank: Geometry (left) and amplitude frequency (right), with 48,000 linear quadrilateral elements.

§.§ Rayleigh-Taylor instability

The Rayleigh-Taylor instability problem described in [48] is considered. The conventional Rayleigh-Taylor problem consists of a fluid $A$ sitting on top of a less dense fluid $B$. Any perturbation between the fluid layers in combination with gravitational force, will cause fluid $A$ to drive into fluid $B$ causing the well know mushroom cloud effect synonymous with Rayleigh-Taylor instability. A domain of size $[0,1]\times[0,4]$ is considered with the higher density fluid having $\rho_1=3$ and $\mu_1 = 0.0031316$, and the lower density fluid having $\rho_2=1$, $\mu_2 = 0.0031316$. The gravitational acceleration is taken as $g=9.80665$. The initial perturbation is given by

\begin{align}
\varphi(x,y) = \tanh\left( \frac{y - 2 - 0.1 \cos(2\pi x)}{\sqrt{2}\epsilon} \right).
\end{align}

The upper and lower boundaries are set to no-slip, and the left and right boundaries are set to slip conditions.

In concurrence with [48], the other parameters are chosen as follows: $\gamma = 0.01$, $\epsilon = 0.005$, $\alpha = \pi/2$. The mobility function is chosen as $M_2(\varphi)$ with $D=4\cdot 10^{-5}$. Two meshes are compared, a mesh with $256\times1024$ linear stabilised elements and a mesh with $128\times512$ mixed Taylor-Hood elements. For comparison we consider the results from Guermond and Quartapelle [81] where surface tension effects are ignored.
It should be mentioned that in [48], the author states that the surface tension coefficient is set as small but not zero, in order to avoid the CH equation becoming a pure transport equation, since setting $\gamma=0$ would render $\eta=0$ in the conventional formulation. In this work, we consider $\gamma=0.01$ with the conventional Abels et al. formulation [4], as well as $\gamma=0$ with the formulation presented in Section <ref>. For the latter formulation the mobility coefficient is increased to $D=10^{-3}$, since it is no longer necessary to set it so small.

The evolution of $\varphi$ is shown in Figure <ref>.

Rayleigh-Taylor instability: Evolution of phase field variable $\varphi$, with $256\times1024$ linear elements.

Figure <ref> shows the position along the vertical axis of the rising and falling interface tips for the stabilised and mixed formulations.

Rayleigh-Taylor instability: $y$ position of interface at the left wall and at the centre. Comparison with [81].

The time scaling of the reference solution [81] required mapping according to $t = \sqrt{2/g} \ \hat{t}$ to account for the non-dimensionality of the variables. The results for the stabilised and mixed formulations are in agreement with the reference, and nearly indistinguishable from each other. Observing the plot corresponding to $\gamma=0$, it appears that neglecting surface tension effects does not really alter the displacements observed, although it does confirm that $\gamma$ can be set to zero without encountering numerical problems.

§ CONCLUSIONS

In this work two novel finite element formulations are presented for modelling the Navier-Stokes-Cahn-Hilliard equations; The first uses mixed Taylor-Hood elements, while the second uses linear equal order stabilised SUPG/PSPG elements. The former formulation is primarily considered for comparative purposes. The models are formulated with two key aspects in mind:

* The ability to deactivate surface tension effects. This is done by removing the surface tension coefficient from the Cahn-Hilliard equations, and hence any instabilities which would result from setting it to zero. Thus, it is ensured that the Cahn-Hilliard equation exclusively deals with phase dynamics, while the Navier-Stokes equations alone control the physical phenomena.

* Computational efficiency, which relates specifically to the stabilised formulation. Here the standard SUPG/PSPG stabilisation is introduced in the momentum equation, which allows for the use of efficient equal order linear elements.
The employment of linear elements greatly improves the overall computational efficiency, which is crucial for the simulation of realistic three dimensional problems.

A number of benchmark/example problems are solved with the proposed methodologies. The examples in Sections <ref>-<ref> demonstrate problems dominated by surface tension, while problems without surface tension are investigated in Sections <ref>-<ref>.
In all cases the available reference solutions are reproduced accurately. Particularly noteworthy is the agreement of the results with those obtained from an arbitrary Lagrangian-Eulerian (ALE) based strategy in [9]. It is demonstrated in Sections <ref> and <ref> that the volumes of the two fluid phases are conserved with good accuracy.

§ REFERENCES

[1]
C. W. Hirt and B. D. Nichols.
Volume of fluid (VOF) method for the dynamics of free boundaries.
Journal of Computational Physics, 39(1):201–225, 1981.

[2]
S. Osher and J. A. Sethian.
Fronts propagating with curvature-dependent speed: Algorithms based
on Hamilton-Jacobi formulations.
Journal of Computational Physics, 79(1):12–49, 1988.

[3]
G. Tryggvason.
Numerical simulations of the Rayleigh-Taylor instability.
Journal of Computational Physics, 75(2):253–282, 1988.

[4]
H. Abels, H. Garcke, and G. Grün.
Thermodynamically Consistent, Frame Indifferent Diffuse Interface
Models for Incompressible Two-Phase Flows with Different Densities.
Mathematical Models and Methods in Applied Sciences,
22(03):1150013, 2011.

[5]
F. Boyer.
A theoretical and numerical model for the study of incompressible
mixture flows.
Computers and Fluids, 31(1):41–68, 2002.

[6]
T. J. R. Hughes, W. K. Liu, and T. K. Zimmermann.
Lagrangian-Eulerian finite element formulation for incompressible
viscous flows.
Computer Methods in Applied Mechanics and Engineering,
29(3):329–349, 1981.

[7]
D. J. Benson.
Computational methods in Lagrangian and Eulerian hydrocodes.
Computer Methods in Applied Mechanics and Engineering,
99(2-3):235–394, 1992.

[8]
J. Donea, A. Huerta, J. Ponthot, and A. Rodriguez-Ferran.
Arbitrary Lagrangian – Eulerian Methods.
Encyclopedia of Computational Mechanics, pages 413–437, 2004.

[9]
W. G. Dettmer and D. Perić.
A computational framework for free surface fluid flows accounting
for surface tension.
Computer Methods in Applied Mechanics and Engineering,
195(23-24):3038–3071, 2006.

[10]
P. Bach and O. Hassager.
An algorithm for the use of the Lagrangian specification in
Newtonian fluid mechanics and applications to free-surface flow.
Journal of Fluid Mechanics, 152:173–190, 1985.

[11]
B. Ramaswamy and M. Kawahara.
Lagrangian finite element analysis applied to viscous free surface
fluid flow.
International Journal for Numerical Methods in Fluids,
7(9):953–984, 1987.

[12]
P. H. Saksono and D. Perić.
On finite element modelling of surface tension Variational
formulation and applications - Part I: Quasistatic problems.
Computational Mechanics, 38(3):265–281, 2006.

[13]
P. H. Saksono and D. Perić.
On finite element modelling of surface tension: Variational
formulation and applications – Part II: Dynamic problems.
Computational Mechanics, 38(3):251–263, 2006.

[14]
A. Javili and P. Steinmann.
A finite element framework for continua with boundary energies. Part
I: The two-dimensional case.
Computer Methods in Applied Mechanics and Engineering,
198(27-29):2198–2208, 2009.

[15]
J. D. van der Waals.
The thermodynamic theory of capillarity under the hypothesis of a
continuous variation of density.
Journal of Statistical Physics, 20(2):197–200, 1979.

[16]
D. J. Korteweg.
Sur la forme que prennent les équations des mouvements des
fluides si l'on tient compte des forces capillaires par des variations de
Arch. N´eer. Sci. Exactes S´er. II, (6):1–24, 1901.

[17]
D. M. Anderson, G. B. McFAdden, and A. A. Wheeler.
Diffuse-Interface Methods in Fluid Mechanics.
Annu. Rev. Fluid Mech, (30):139–65, 1998.

[18]
H. Gomez and K. G. van der Zee.
Computational Phase-Field Modeling.
In Encyclopedia of Computational Mechanics Second Edition,
chapter 11, pages 1–35. John Wiley & Sons, Ltd, Chichester, UK, second
edition, 2017.

[19]
J. W. Cahn and J. E. Hilliard.
Free energy of a nonuniform system. III. Nucleation in a
two-component incompressible fluid.
The Journal of Chemical Physics, 31(3):688–699, 1959.

[20]
S. M. Allen and J. W. Cahn.
Ground state structures in ordered binary neighbor interactions.
Acta Metallurgica, 20:423–433, 1972.

[21]
G. N. Wells, E. Kuhl, and K. Garikipati.
A discontinuous Galerkin method for the Cahn-Hilliard equation.
Journal of Computational Physics, 218(2):860–877, 2006.

[22]
H. Gómez, V. M. Calo, Y. Bazilevs, and T. J. R. Hughes.
Isogeometric analysis of the Cahn-Hilliard phase-field model.
Computer Methods in Applied Mechanics and Engineering,
197(49-50):4333–4352, 2008.

[23]
H. Gómez and T. J. R. Hughes.
Provably unconditionally stable, second-order time-accurate, mixed
variational methods for phase-field models.
Journal of Computational Physics, 230(13):5310–5327, 2011.

[24]
F. Guillén-González and G. Tierra.
On linear schemes for a cahn-hilliard diffuse interface model.
Journal of Computational Physics, 234(1):140–171, 2013.

[25]
L. K. Antanovskii.
A phase field model of capillarity.
Physics of Fluids, 7(4):747–753, 1995.

[26]
M. E. Gurtin.
Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a
microforce balance.
Physica D: Nonlinear Phenomena, 92(3-4):178–192, 1996.

[27]
D. Jacqmin.
Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field
Journal of Computational Physics, 155(1):96–127, 1999.

[28]
J. S. Lowengrub and L. Truskinovsky.
Quasi-incompressible Cahn-Hilliard fluids and topological
Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 454(1978):2617–2654, 1998.

[29]
J. Shen and X. Yang.
A Phase-Field Model and Its Numerical Approximation for Two-Phase
Incompressible Flows with Different Densities and Viscosities.
SIAM Journal on Scientific Computing, 32(3):1159–1179, 2010.

[30]
P. Hohenberg and B. Halperin.
Theory of dynamic critical phenomena.
Reviews of Modern Physics, 49(3):435–479, 1977.

[31]
M. Shokrpour Roudbari, G. Simsek, E. H. van Brummelen, and K. G. van der Zee.
Diffuse-interface two-phase flow models with different densities: A
new quasi-incompressible form and a linear energy-stable method.
Mathematical Models and Methods in Applied Sciences,
28(04):733–770, 2018.

[32]
H. Ding, P. D. M. Spelt, and C. Shu.
Diffuse interface model for incompressible two-phase flows with
large density ratios.
Journal of Computational Physics, 226(2):2078–2095, 2007.

[33]
X. Xu, L. Zhao, and C. Liu.
Axisymmetric Solutions to Coupled Navier-Stokes/Allen-Cahn
SIAM Journal on Mathematical Analysis, 41(6):2246–2282, 2010.

[34]
V. Joshi and R. K. Jaiman.
A positivity preserving and conservative variational scheme for
phase-field modeling of two-phase flows.
Journal of Computational Physics, 360:137–166, 2018.

[35]
G. Favre and G. Schimperna.
On a Navier-Stokes-Allen-Cahn model with inertial effects.
Journal of Mathematical Analysis and Applications,
475(1):811–838, 2019.

[36]
S. Dong and J. Shen.
A time-stepping scheme involving constant coefficient matrices for
phase-field simulations of two-phase incompressible flows with large density
Journal of Computational Physics, 231(17):5788–5804, 2012.

[37]
H. G. Lee, J. S. Lowengrub, and J. Goodman.
Modeling pinchoff and reconnection in a Hele-Shaw cell. II. Analysis
and simulation in the nonlinear regime.
Physics of Fluids, 14(2):514–545, 2002.

[38]
C. Liu and J. Shen.
A phase field model for the mixture of two incompressible fluids and
its approximation by a Fourier-spectral method.
Physica D: Nonlinear Phenomena, 179(3-4):211–228, 2003.

[39]
J. Kim, K. Kang, and J. Lowengrub.
Conservative multigrid methods for Cahn-Hilliard fluids.
Journal of Computational Physics, 193(2):511–543, 2004.

[40]
X. Feng.
Fully discrete finite element approximations of the
Navier-Stokes-Cahn- Hilliard diffuse interface model for two-phase fluid
SIAM Journal on Numerical Analysis, 44(3):1049–1072, 2006.

[41]
Q. Du and J. Zhang.
Adaptive Finite Element Method for a Phase Field Bending Elasticity
Model of Vesicle Membrane Deformations.
SIAM Journal on Scientific Computing, 30(3):1634–1657, 2008.

[42]
D. Kay, V. Styles, and R. Welford.
Finite element approximation of a Cahn-Hilliard-Navier-Stokes
Interfaces and Free Boundaries, 10(1):15–43, 2008.

[43]
Y. Zhang, H. Wang, and T. Tang.
Simulating two-phase viscoelastic flows using moving finite element
Communications in Computational Physics, 7(2):333–349, 2010.

[44]
K. Bao, Y. Shi, S. Sun, and X. Wang.
A finite element method for the numerical solution of the coupled
Cahn–Hilliard and Navier–Stokes system for moving contact line problems.
Journal of Computational Physics, 231(24):8083–8099, 2012.

[45]
Z. Guo, P. Lin, and J. S. Lowengrub.
A numerical method for the quasi-incompressible
Cahn-Hilliard-Navier-Stokes equations for variable density flows with a
discrete energy law.
Journal of Computational Physics, 276:486–507, 2014.

[46]
J. Giesselmann and T. Pryer.
Energy consistent discontinuous Galerkin methods for a
quasi-incompressible diffuse two phase flow model.
ESAIM: Mathematical Modelling and Numerical Analysis,
49(1):275–301, 2015.

[47]
L. F. R. Espath, A. F. Sarmiento, P. Vignal, B. O. N. Varga, A. M. A. Cortes,
L. Dalcin, and V. M. Calo.
Energy exchange analysis in droplet dynamics via the
Navier-Stokes-Cahn-Hilliard model.
Journal of Fluid Mechanics, 797:389–430, 2016.

[48]
B. S. Hosseini, S. Turek, M. Möller, and C. Palmes.
Isogeometric Analysis of the Navier-Stokes-Cahn-Hilliard equations
with application to incompressible two-phase flows.
Journal of Computational Physics, 348:171–194, 2017.

[49]
S. Minjeaud.
An unconditionally stable uncoupled scheme for a triphasic
Cahn-Hilliard/Navier-Stokes model.
Numerical Methods for Partial Differential Equations,
29(2):584–618, 2013.

[50]
C. Liu, J. Shen, and X. Yang.
Decoupled Energy Stable Schemes for a Phase-Field Model of Two-Phase
Incompressible Flows with Variable Density.
Journal of Scientific Computing, 62(2):601–622, 2015.

[51]
Y. Chen and J. Shen.
Efficient, adaptive energy stable schemes for the incompressible
Cahn-Hilliard Navier-Stokes phase-field models.
Journal of Computational Physics, 308:40–56, 2016.

[52]
F. Guillén-González and G. Tierra.
Splitting schemes for a Navier-Stokes-Cahn-Hilliard model for two
fluids with different densities.
Journal of Computational Mathematics, 32(6):643–664, 2014.

[53]
X. Yang, M. G. Forest, H. Li, C. Liu, J. Shen, Q. Wang, and F. Chen.
Modeling and simulations of drop pinch-off from liquid crystal
filaments and the leaky liquid crystal faucet immersed in viscous fluids.
Journal of Computational Physics, 236(1):1–14, 2013.

[54]
P. H. Chiu.
A coupled phase field framework for solving incompressible two-phase
Journal of Computational Physics, 392:115–140, 2019.

[55]
J. W. Cahn and J. E. Hilliard.
Free energy of a nonuniform system. I. Interfacial free energy.
The Journal of Chemical Physics, 28(2):258–267, 1958.

[56]
C. M. Elliott and S. Zheng.
On the Cahn-Hilliard equation.
Nonlinear Analysis, 96:339–357, 1986.

[57]
J. F. Blowey and C. M. Elliott.
The Cahn-Hilliard gradient theory for phase separation with
non-smooth free energy Part II: Numerical analysis.
European Journal of Applied Mathematics, 3(2):147–179, 1992.

[58]
W. Villanueva and G. Amberg.
Some generic capillary-driven flows.
International Journal of Multiphase Flow, 32(9):1072–1086,

[59]
X. Xu, Y. Di, and H. Yu.
Sharp-interface limits of a phase-field model with a generalized
Navier slip boundary condition for moving contact lines.
Journal of Fluid Mechanics, 849:805–833, 2018.

[60]
H. Yu and X. Yang.
Numerical approximations for a phase-field moving contact line model
with variable densities and viscosities.
Journal of Computational Physics, 334:665–686, 2017.

[61]
P W Voorhees.
Ostwald Ripening of Two-Phase Mixtures.
Annual Review of Materials Science, 22(1):197–215, 1992.

[62]
T. E. Tezduyar, M. Behr, S. Mittal, and J. Liou.
A new strategy for finite element computations involving moving
boundaries and interfaces-The deforming-spatial-domain/space-time procedure:
II. Computation of free-surface flows, two-liquid flows, and flows with
drifting cylinders.
Computer Methods in Applied Mechanics and Engineering,
94(3):353–371, 1992.

[63]
T. E. Tezduyar and S. Sathe.
Stabilization Parameters in SUPG and PSPG Formulations.
Journal of Computational and Applied Mechanics, 4(1):71–88,

[64]
W. G. Dettmer and D. Perić.
An analysis of the time integration algorithms for the finite
element solutions of incompressible Navier - Stokes equations based on a
stabilised formulation.
Computer Methods in Applied Mechanics and Engineering,
192(9-10):1177–1226, 2003.

[65]
G. Scovazzi.
A discourse on Galilean invariance, SUPG stabilization, and the
variational multiscale framework.
Computer Methods in Applied Mechanics and Engineering,
196(4-6):1108–1132, 2007.

[66]
T. J. R. Hughes, G. Scovazzi, and L. P. Franca.
Multiscale and Stabilized Methods.
In Encyclopedia of Computational Mechanics Second Edition,
Volume 5: Fluids - Part 1, pages 5–68. 2017.

[67]
W. G. Dettmer, C. Kadapa, and D. Perić.
A stabilised immersed boundary method on hierarchical b-spline
Computer Methods in Applied Mechanics and Engineering,
311:415–437, 2016.

[68]
J. Chung and G. M. Hulbert.
A Time Integration Algorithm for Structural Dynamics With Improved
Numerical Dissipation: The Generalized-$\alpha$ Method.
Journal of Applied Mechanics, 60(2):371, 2008.

[69]
K. E. Jansen, C. H. Whiting, and G. M. Hulbert.
Generalized-$\alpha$ method for integrating the filtered
Navier-Stokes equations with a stabilized finite element method.
Computer Methods in Applied Mechanics and Engineering,
190(3-4):305–319, 2000.

[70]
A. Lovrić, Wulf G. Dettmer, Chennakesava Kadapa, and Djordje Perić.
A new family of projection schemes for the incompressible
Navier–Stokes equations with control of high-frequency damping.
Computer Methods in Applied Mechanics and Engineering,
339:160–183, 2018.

[71]
C. Kadapa, W. G. Dettmer, and D. Perić.
A stabilised immersed boundary method on hierarchical b-spline grids
for fluid–rigid body interaction with solid–solid contact.
Computer Methods in Applied Mechanics and Engineering,
318:242–269, 2017.

[72]
C. Kadapa, W. G. Dettmer, and D. Perić.
A stabilised immersed framework on hierarchical b-spline grids for
fluid-flexible structure interaction with solid–solid contact.
Computer Methods in Applied Mechanics and Engineering,
335:472–489, 2018.

[73]
A. De Coninck, B. De Baets, D. Kourounis, F. Verbosio, O. Schenk,
S. Maenhout, and J. Fostier.
Needles: Toward Large-Scale Genomic Prediction with
Marker-by-Environment Interaction.
Genetics, 203(1):543–555, 2016.

[74]
F. Verbosio, A. De Coninck, D. Kourounis, and O. Schenk.
Enhancing the scalability of selected inversion factorization
algorithms in genomic prediction.
Journal of Computational Science, 22(Supplement C):99–108,

[75]
D. Kourounis, A. Fuchs, and O. Schenk.
Towards the Next Generation of Multiperiod Optimal Power Flow
IEEE Transactions on Power Systems, PP(99):1–10, 2018.

[76]
S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith.
Efficient Management of Parallelism in Object Oriented Numerical
Software Libraries.
In E Arge, A M Bruaset, and H P Langtangen, editors, Modern
Software Tools in Scientific Computing, pages 163–202. Birkhäuser
Press, 1997.

[77]
W. G. Dettmer, P. H. Saksono, and D. Perić.
On a finite element formulation for incompressible Newtonian fluid
flows on moving domains in the presence of surface tension.
Communications in Numerical Methods in Engineering,
19(9):659–668, 2003.

[78]
G. B. Foote.
A numerical method for studying liquid drop behavior: Simple
Journal of Computational Physics, 11(4):507–530, 1973.

[79]
S. Liu, S. Li, and J. Liu.
Jurin's law revisited: Exact meniscus shape and column height.
European Physical Journal E, 41(3):1–7, 2018.

[80]
C. Pozrikidis and J. H. Ferziger.
Introduction to Theoretical and Computational Fluid Dynamics.
Physics Today, 50(9):72–74, 1997.

[81]
J. L. Guermond and L. Quartapelle.
A Projection FEM for Variable Density Incompressible Flows.
Journal of Computational Physics, 165(1):167–188, 2000.


11institutetext:  GHN Data Science, Ho Chi Minh, Vietnam
11email<EMAIL_ADDRESS>
# Hub and Spoke Logistics Network Design for Urban Region with Clustering-
Based Approach

Quan Duong  0000-0002-4493-2331 Dang Nguyen  0000-0002-3664-4364 Quoc Nguyen
0000-0002-1692-1918

###### Abstract

This study aims to propose effective modeling and approach for designing a
logistics network in the urban area in order to offer an efficient flow
distribution network as a competitive strategy in the logistics industry where
demand is sensitive to both price and time [3]. A multi-stage approach is
introduced to select the number of hubs and allocate spokes to the hubs for
flow distribution and hubs’ location detection. Specifically, a fuzzy
clustering model with the objective function is to minimize the approximate
transportation cost is employed, in the next phase is to focus on balancing
the demand capacity among the hubs with the help of domain experts, afterward,
the facility location vehicle routing problems within the network is
introduced. To demonstrate the approach’s advantages, an experiment was
performed on the designed network and its actual transportation cost for the
real operational data in which specific to the Ho Chi Minh city infrastructure
conditions. Additionally, we show the flexibility of the designed network in
the flow distribution and its computational experiments to develop the
managerial insights which contribute to the network design decision-making
process.

††Proceedings of the 34th International Conference on Industrial, Engineering
& Other Applications of Applied Intelligent Systems (IEA/AIE 2021), Kuala
Lumpur, Malaysia. Copyright 2021 by the author(s)

###### Keywords:

Logistics network design Hub-and-spoke Clustering.

##  1 Introduction

The volume of express delivery gets increasing with the rapid development of
e-commerce. At the same time, intense competition has forced express delivery
service companies to compete in different categories to get ahead in the race.
Transportation cost and delivery lead-time are considered to be the most
important metrics for success. On the other hand, logistics network design is
concerned with the number of locations of stations that allocate for
customer’s demand points, and the assignment of the stations to the
distribution centers within their coverage. The optimal setting must
distribute the goods to the customers with the least cost and satisfy the
service level agreement. To illustrate these strategies, considering the
current configuration of the delivery network of logistics providers in
Vietnam, an attempt to design a network that distributes more efficiently and
reduces the transportation cost. Thus, the primary goal of this work is to,
firstly, examine the limitation of the single distribution center. Secondly,
design a delivery network that optimizes the transportation cost by
demonstrating for an urban region. Moreover, we develop an analytical
framework that offers multiple distribution scenarios within the designed
network and its impact on the operational performance of the firm.

The rest of the paper is organized as follows: In Section 2, the relevant
studies on logistics network design for the urban region are reviewed.
Afterward, Section 3 presents our practical approach and framework and
establishes the flow distribution scenarios within the network. The experiment
results take into account the real operational data from a delivery network
with a case study for Ho Chi Minh city will be conducted in Section 4.
Finally, we conclude the main findings and future research directions in
Section 5.

##  2 Literature Review

The logistics network design problem consisting of multiple facilities
location problem [6], each facility has a limit of serving demand points
contribute the difficulty of problem where need to solve the single facility
location problems, but also searching for the optimal number of facilities and
assigning the customer demand points for each cluster. Additionally, Cheong et
al [3] explored the benefit of segmenting demand points into different
classes, they assumed that each demand point can be divided into two classes
of the sensitivity to delivery lead time, the segmentation results in more
effective allocation in demand point and reduce the network cost. A typical
use of flow distribution, utilizes the distribution centers (DC) to serve the
customer’s demand points within their coverage. In a single DC scenario,
vehicles depart from DC’s location and deliver to the customer’s demand
points, this process makes the logistics resource difficult to be shared in
high-density areas. Moreover, the logistics service provider manages to
accommodate customer’s demanding at time window and price make the operations
more challenging. In 2016, Roca-Riu et al [9] took a further step to evaluate
the urban good distribution for several European cities with a single
independent DC and proposed continuous models to improve the efficiency of
urban distribution with the use of urban consolidation centers (UCC), the
quantitative metric indicated that by deploying the UCC can reduce the
transportation cost up to 14%. The collaborative and multi-depots in urban
areas got a lot of attention in recent years [7, 10, 12]. Later in 2020, Wang
et al [11] solve collaborative multi-depots vehicle routing problems to deal
with the impact of changing time window. In order to group the demand points
into a single, weighted K-mean clustering and Fuzzy C-means clustering
algorithms stand out for the facility location problem [4, 6, 5]. However,
Meng et al [13] figured out that weighted K-mean always traps into local
optima and is sensitive to the initial setting, and very limited research
incorporate the domain knowledge for the target region with fuzzy c-means
clustering capacity.

In this study, the fuzzy c-mean (FCM) clustering algorithm [2] is employed to
solve the uncapacitated multi-facility location problem. The FCM allows the
data point belongs to multiple clusters with different degrees of membership.
Combining the degree of membership features and domain expert for the target
region of designing the delivery network, the belonging of demand point to the
cluster will be judged by humans in order to balance the capacity among the
clustering and other area’s characteristics.

##  3 Methodology

Given all the spokes (know as the demand points that will serve the end
customers for a specific area). The methodology aims to find the optimal Hub’s
Location and spoke to the hub assignment with respect to minimize the
transportation cost. The process consisting of three phases: i) Clustering:
detecting the group of nearby spokes to a cluster. ii) Hub’s location: detect
the cluster’s center in which serve as the local DC. iii) Vehicle Routing:
determine the flow distribution scenario among the Hubs’ network, and perform
the vehicle routing for good delivery to the spokes.

###  3.1 Clustering

Formally, a dataset X as a set of N demand points represented as vectors in an
2-dimensional space: $X=\\{x_{1},x_{2},...,x_{N}\\}\subseteq\Re^{\textit{2}}$.
A clustering in X is a set of disjoint clusters that partitions X into K
groups: $C=\\{c_{1},c_{2},...,c_{k}\\}$ where $\cup_{c_{k}\in
C}c_{k}=X,c_{k}\cap c_{l}=\emptyset\forall k\neq l$. The Haversine distance
between objects $x_{i}$ and $x_{j}$ will denote as $d_{h}(x_{i},x_{j})$. The
approximate transportation cost is calculated as

$Cost=\sum_{c_{k}\in C}\left\\{\sum_{c_{l}\in C\setminus
c_{k}}d_{h}(\overline{c_{k}},\overline{c_{l}})+\sum_{x_{i}\in
c_{k}}d_{h}(\overline{c_{k}},x_{i})\right\\}$ (1)

$\overline{c_{k}}=\frac{\sum_{x_{i}\in
X}{w_{k}(x_{i})}^{m}*x_{i}}{\sum_{x_{i}\in X}{w_{k}(x_{i})}^{m}}$ (2)

Where $\overline{c_{k}}$ is the centroid of a cluster $c_{k}$ is defined by
the FCM, m is the hyper-parameter that controls the fuzziness of clustering
and $w_{k}(x)$ is the degree of membership that data point x belong to the
cluster k.

####  3.1.1 Final Cluster Assignment

The degree of membership feature from FCM is utilized for balancing the demand
among the cluster together with humans. Specifically, taking into account the
domain knowledge of the target designed area, those mixture-membership demand
data points will be re-assigned with the domain expert in the loop. Based on
the FCM’s degree of membership feature, the expert will assign each single
demand data point into a target cluster in order to satisfy the capacity
balancing requirements. The number of clusters and algorithm parameters is
searching by incremental methods with the objective is to minimize the
approximate transportation cost in the formula 1. In the Section 4.1, the
experiment results with specific parameters will be conducted.

###  3.2 Hub’s location

After allocating of demand data point to cluster, the center of gravity method
[1] is employed to determine the optimal Hub’s location. The aims of center of
gravity method is to bring the center of cluster closer to demand points with
higher demand. For all $\textit{k}\in\\{1,2,...,K\\}$ and $\nu(x)$ is the
delivery demand for demand point x, the center gravity of each cluster k are
calculated as:

${h_{k}}=\sum_{x\in c_{k}}{\beta(x)*x}$ (3) $\beta(x)=\frac{\nu(x)}{\sum_{x\in
c_{k}}{\nu(x)}}$ (4)

The $\beta(x)$ value serves as the weighted factor within the cluster. The
larger value will pull the hub’s location to its side. The real transportation
cost corresponding to these locations will be shown in Table 2.

###  3.3 Flow Distribution and Vehicle Routing

Figure 1: Flow Distribution Scenarios

Given a designed network, the authors deliberately prototype three flow
distribution scenarios within the network in the Figure 1, and a default
scenario which is using for daily operation in GHN network as follows: $S0$:
Default scenario, single independent DC. $S1$: Each Hub will serve as Urban
Consolidation Center (UCC). $S2$: Each Hub will serve as Urban Distribution
Center (UDC). $S3$: Each Hub will serve as UDC, and drop the delivery packages
at its center for the nearby demand points. Here, the newly designed network
has the ability to offer three different flow distribution scenarios.
Fundamentally, $S1$ only performs the loading/unloading packages to the
designated DC. An upgrade version of $S1$, $S2$ also has the ability to sort
the packages within its coverage and distribute all other packages to other
DC/UDC. As an advantage on multi-depot, the DC locate closer to its demand
points, $S3$ only need to sort packages, and hand-off the packages at its
location, this scenario offers flexibility for the urban area such as reduce
the cost for delivering by trucks, restricted during the rush hours.

##  4 Experiment Result

The numerical results reported in this paper are the real operational data in
Ho Chi Minh City, Vietnam. Including, i) 77 active demand data points. ii) All
transactions data for the entire December 2020 has been extracted and
calculated for two different phases in this paper. First, the average delivery
demand for the whole period would be used for representing the demand point’s
volume that contributes to the Hub’s Location inference in Section 3.2.
Second, transaction data of $22^{nd}December$ has been used to conduct
experimental design for pickup and delivery demand in the designed network.
Essentially, the actual transportation cost in Section 4.2 are computed based
on this day.

###  4.1 Detected Clusters

Regarding to experimental setting, the authors employed the Python’s scikit
fuzzy package222https://pythonhosted.org/scikit-fuzzy/ and the configuration
($\textit{m}=3$, $\textit{error}=0.002$, $\textit{maxiter}=1000$,
$\textit{seed}=12345$, and $\textit{c}=\\{2,3,4,5\\}$). The primary parameter
is c (the desired number of clusters). The best parameter in determined by the
combination of approximate transportation cost in Formula 1, fuzzy
partitioning coefficient, and the variation of delivery demand for each
cluster.

Table 1: Summary of Clustering Results Number of clusters | Approximate transportation cost | Coef
---|---|---
2 | 527.8 | 0.568
3 | 482.8 | 0.424
4 | 490.3 | 0.362
5 | 530.6 | 0.308
Figure 2: The environment and its clusters between the elements in network.

The experimental results of the clustering step shown in Table 1, The
transportation cost with 3 and 4 show the better result than the other values.
Take a deeper look, we can see that the 3 clusters’s cost and Coef slightly
better than the 4. Inspiring by these numbers, the exact assignment of data
points to the cluster will be adjusted by the domain expert due to the fact
that it’s mandatory the satisfy the capacity requirements. Figure 2 shows the
visual outcome of each step for the entire process. The initial setting of 77
demand points display as black dots belong to the Ho Chi Minh region in Sub-
figure 2-a, each color in the right next Sub-figure 2-b represents the 3
detected clusters by virtue of the FCM algorithm and human judgment. A single
triangular in Sub-figure 2-c indicates the Hub’s location for each cluster.
The right-bottom Sub-figure 2-d visually shows the flow distribution from Hub
to Hub and from Hub to Spokes within sub-network.

###  4.2 Actual Transportation Cost

Confidentially, with the goal of comparing the transportation cost among the
scenarios, also not to disclose the commercial confidentiality information.
The values for number of trucks and total cost for the current operation mode
have been changed to v and c respectively. The actual transportation cost
specified in the Formula 1. Due to the fact that the actual situation takes
into account many important and realistic factors like delivering demand at
each demand point, truck capacity, time window requirement, local traffic
condition. Thus, we leverage the capacitated vehicle routing with time window
solutions. In this regard, the optimization software OR-Tool [8] is used to
perform the computational step given the operational constraints and the
conditions of Ho Chi Minh city such as travel distance and duration between
points, time window requirement, delivery demand, and truck’s capacity.

Table 2: Actual Transportation Cost for Scenarios

Scenario | Number of trucks | Total cost | Pickup cost | Delivery cost
---|---|---|---|---
$S0$ |  $v$ * |  $c$ ** | 0.7c | 0.3c
$S1$ | 0.79$v$ | 0.75$c$ | 0.44$c$ | 0.31$c$
$S2$ | 0.8$v$ | 0.72$c$ | 0.47$c$ | 0.25$c$
$S3$ | 0.63$v$ | 0.61$c$ | 0.46$c$ | 0.15$c$

  * *

Number of trucks need for the S0 scenario

  * **

The total amount of money take for the S0 scenario

In the Table 2, we shown the actual transportation cost for all the mentioned
scenarios. By introducing the multi-depot in the delivery network, the
transportation cost reduces 25%, 28% by adding consolidation and distribution
centers respectively. Take a deeper look in the S2 scenario, we can observe
that the delivery cost is cheaper than S1 scenario 19% by reducing the
traveled distance for the delivery trips, S2 also takes more effort in the
pickup trips in order to move the packages to closer DC. On the other side,
about 20% of the number of trucks have been reduced for performing delivery
with the same demand in the network.

##  5 Conclusions

This paper analyzes the advantages of collaborative multi-depot settings in
the urban area, applied the iterative approaches for the logistics network
design problem. By employing the FCM algorithm for clustering and leveraging
the degree of membership feature to reassign the demand point to the cluster
in order to balance the facility’s capacity in each specific sub-network by
virtue of domain expert. Based on the real operational data in Ho Chi Minh
city, the results indicate that 3 clusters achieve the best solution in
approximate transportation cost and reduce the actual transportation cost to
28% and 20% trucks used with 3 distribution centers compared to the original
independent DC. Increasing the truckload utilization for the delivery trips
between Hubs not only reduces the transportation cost and improves operational
efficiency for the firm, but also put an effort into alleviating the traffic
congestion, reducing stress for transportation infrastructure. In the future,
the studies may focus on two directions. First, calculate the system cost for
the designed network including the cost for opening a new facility and
operating that facility with a given demand, and perform the experiment for
multiple urban regions in Vietnam. Second, schedule the dynamic vehicle
routing to deal with the highest demand days and measure its performance.
Prospective studies on these directions would contribute to operating more
efficiently and resilient on both typical and high demand days, also the
promotion of multi-depots setting for the urban regions.

## References

  * [1] Ballou, R.H.: Business logistics management. Upper Saddle River, New Jersey: Prentice-Hall (1999)
  * [2] Bezdek, J.C., Ehrlich, R., Full, W.: Fcm: The fuzzy c-means clustering algorithm. Computers & Geosciences 10(2-3), 191–203 (1984)
  * [3] Cheong, M.L., Bhatnagar, R., Graves, S.C.: Logistics network design with differentiated delivery lead-time: Benefits and insights (2005)
  * [4] Esnaf, Ş., Küçükdeniz, T.: A fuzzy clustering-based hybrid method for a multi-facility location problem. Journal of Intelligent Manufacturing 20, 259–265 (2009)
  * [5] Hu, W., Dong, J., Hwang, B.G., Ren, R., Chen, Z.: Network planning of urban underground logistics system with hub-and-spoke layout: two phase cluster-based approach. Engineering, Construction and Architectural Management (2020)
  * [6] Küçükdeniz, T., Baray, A., Ecerkale, K., Esnaf, Ş.: Integrated use of fuzzy c-means and convex programming for capacitated multi-facility location problem. Expert Systems with Applications 39(4), 4306–4314 (2012)
  * [7] Long, Q.: A multi-methodological collaborative simulation for inter-organizational supply chain networks. Knowledge-Based Systems 96, 84–95 (2016)
  * [8] Perron, L., Furnon, V.: Or-tools, https://developers.google.com/optimization/
  * [9] Roca-Riu, M., Estrada, M., Fernández, E.: An evaluation of urban consolidation centers through continuous analysis with non-equal market share companies. Transportation Research Procedia 12, 370–382 (2016)
  * [10] Wang, Y., Ma, X., Liu, M., Gong, K., Liu, Y., Xu, M., Wang, Y.: Cooperation and profit allocation in two-echelon logistics joint distribution network optimization. Applied Soft Computing 56, 143 – 157 (2017)
  * [11] Wang, Y., Zhang, S., Guan, X., Peng, S., Wang, H., Liu, Y., Xu, M.: Collaborative multi-depot logistics network design with time window assignment. Expert Systems with Applications 140, 112910 (2020)
  * [12] Xu, X.F., Hao, J., Deng, Y.R., Wang, Y.: Design optimization of resource combination for collaborative logistics network under uncertainty. Applied Soft Computing 56, 684–691 (2017)
  * [13] You, M., Xiao, Y., Zhang, S., Yang, P., Zhou, S.: Optimal mathematical programming for the warehouse location problem with euclidean distance linearization. Computers & Industrial Engineering 136, 70–79 (2019)