*[1]  *[1]  *[1]

# A level set-based solver for two-phase incompressible flows: Extension to
magnetic fluids

Paria Makaremi-Esfarjani , Andrew J. Higgins , Alireza Najafi-Yazdi
Department of Mechanical Engineering, McGill University,
Montréal, Québec, Canada Email address for correspondence:
<EMAIL_ADDRESS>address for correspondence:
<EMAIL_ADDRESS>
###### Abstract

Development of a two-phase incompressible solver for magnetic flows in the
magnetostatic case is presented. The proposed numerical toolkit couples the
Navier–Stokes equations of hydrodynamics with Maxwell’s equations of
electromagnetism to model the behaviour of magnetic flows in the presence of a
magnetic field. To this end, a rigorous implementation of a second-order two-
phase solver for incompressible nonmagnetic flows is introduced first. This
solver is implemented in the finite-difference framework, where a fifth-order
conservative level set method is employed to capture the evolution of the
interface, along with an incompressible solver based on the projection scheme
to model the fluids. The solver demonstrates excellent performance even with
high density ratios across the interface (Atwood number $\approx 1$), while
effectively preserving the mass conservation property. Subsequently, the
numerical discretization of Maxwell’s equations under the magnetostatic
assumption is described in detail, utilizing the vector potential formulation.
The primary second-order solver for two-phase flows is extended to the case of
magnetic flows, by incorporating the Lorentz force into the momentum equation,
accounting for high magnetic permeability ratios across the interface. The
implemented solver is then utilized for examining the deformation of
ferrofluid droplets in both quiescent and shear flow regimes across various
susceptibility values of the droplets. The results suggest that increasing the
susceptibility value of the ferrofluid droplet can affect its deformation and
rotation in low capillary regimes. In higher capillary flows, increasing the
magnetic permeability jump across the interface can further lead to droplet
breakup as well. The effect of this property is also investigated for the
Rayleigh–Taylor instability growth in magnetic fluids.

Keywords:
Two-phase flows
Incompressible flow
Magnetostatic
Magnetic flows
Level set method
Ferrofluid droplets
Rayleigh–Taylor instability

##  1  Introduction

Modelling multi-phase (interfacial) flows involves simulating systems with the
presence of two or more immiscible fluids with different physical properties
and distinguishable interfaces. Multi-phase flows, particularly two-phase
flows, have garnered substantial interest in various applications, including
spray atomization [1], bubbly flows [2], and nuclear reactors [3]. The
complexity of analytically and experimentally studying the physics of two-
phase flows [4] has necessitated the development of accurate, cost-effective,
and consistent numerical methods. Numerous studies have been conducted to
establish a comprehensive numerical toolkit in the field. Despite these
efforts, developing a numerical solver to study two-phase flows remains a
challenging task. The difficulty originates from modelling fluid property
discontinuities across thin interfaces, particularly when large density ratios
are present. The inability to accurately capture discontinuities in fluid
properties can result in numerical instabilities. Errors arising from the
inadequate discretization of fluid property discontinuities can become
particularly pronounced for large density ratios, thereby limiting simulations
to low density ratios. However, most realistic problems of interest involve
large density ratios, such as the formation and dynamics of bubbles, molten
metal flows in atmospheric air, gas entrainment in liquid phases, and the
aerodynamic effects of gas on the liquid phase [5]. Enforcing mass and
momentum conservation is also crucial to obtain physical results and to avoid
numerical instabilities in simulations. Therefore, the numerical solver must
be able to provide consistent mass and momentum exchange across interfaces
throughout the simulation. It is essential for the two-phase solver to be able
to address topology changes of the interface and accommodate a wide range of
time and length scales as well. In numerical solvers for two-phase flows, it
is crucial to properly couple the governing equations of fluids with an
appropriate interface-tracking method.

Existing techniques for simulating two-phase flows can be grouped into four
categories: the lattice Boltzmann method (LBM), smoothed particle
hydrodynamics (SPH), two-fluid, and one-fluid models. The first two methods
examine the behaviour of the fluid by representing it as a collection of
particles [6, 7, 8, 9], while the last two approaches assume the fluid as a
continuum medium that can be described by solving the Navier–Stokes equations.
The two-fluid method treats the two phases as separate fluids that interact
with each other. This model has demonstrated success in simple problems but
has been shown to be inadequate for more complex scenarios [4]. In this study,
we adopt the one-fluid formulation and implement it on a fixed Eulerian grid
to develop a two-phase incompressible numerical solver. The Navier–Stokes
equations are solved for the entire computational domain, accounting for the
density and viscosity jump at the interface while implicitly imposing
appropriate boundary conditions across the interface separating the two fluid
regions. Among different numerical schemes for treating high density ratios
across the interface and modelling surface tension forces, the ghost fluid
method (GFM) [10] and the continuum surface force (CSF) [11] method stand out
as robust solutions. The GFM is based on a generalized Taylor series expansion
and explicitly accounts for the density jump at the interface; as a result, it
is not sensitive to the amplitude of the density jump. The surface tension
force is also directly incorporated in the pressure jump condition, leading to
a sharp numerical treatment of this singular term [1]. However, in the CSF
approach, instead of including the surface tension force directly in the
pressure jump condition, this force is represented as a volumetric force
spread over a few grid points surrounding the interface. While this approach
may result in a slightly less accurate interface representation, particularly
in cases with small front structures, it is generally considered to be less
numerically challenging. Furthermore, discretizing the viscous terms using the
GFM can be difficult and complex to implement numerically, making it less
desirable for some applications [1]. Thus, many researchers employ the CSF
approach to discretize the viscous term. In this study, we used the CSF method
to model both the surface tension and the viscosity terms. Albeit slightly
less accurate as compared to GFM for surface tension modelling, the CSF
approach is more straightforward to implement and can provide robust and
accurate results.

The available methods to numerically transport an interface can be divided
into two categories: interface-tracking and interface-capturing [12]. One
well-known approach in the interface-tracking category is the front-tracking
method introduced by Unverdi and Tryggvason [13]. This method involves
breaking down the fluid interface into discrete material points, referred to
as front-tracking points, which are then transported using a moving mesh that
follows a Lagrangian approach. While this approach benefits from purely
Lagrangian transport, it faces difficulties in preserving liquid volume due to
the requirement for frequent mesh rearrangements [1]. Additionally,
parallelization of the front-tracking method presents a significant challenge.
Furthermore, any break-up or merging of the interface should be addressed
manually due to the inability of this technique to inherently handle topology
changes. As a result, front-tracking methods are not well-suited for
simulations with frequent topological changes, such as primary atomization
[1].

Interface-capturing methods such as the volume-of-fluid (VOF) method [14] and
level set method [15] implicitly capture the interface and can robustly
address complex topological changes in the simulation. The VOF method employs
a liquid volume fraction transport equation to depict the interface, ensuring
mass conservation. However, since the VOF scalar is discontinuous across the
interface, specific numerical treatments are required for the discretization
of the transport equation. The discontinuous nature of the VOF scalar presents
difficulties in computing interface properties such as normal and curvature
values as well.

The level set method, introduced by Sethian [15] in the field of image
processing and computer graphics, represents an interface implicitly using the
iso-level of a smooth function, i.e., the signed distance function. The
smoothness of the level set function is maintained with the re-initialization
process, and the Eulerian scalar transport equation can be solved using high-
order numerical schemes. In addition, parallelization of the solver can be
accomplished efficiently, and interface characteristics such as normal and
curvature values are easily calculable due to the smoothness of the level set
function. Despite all the mentioned advantages of the level set method, this
method does not inherently conserve mass during the simulations, leading to
potentially significant errors. Various hybrid methods have been introduced to
overcome the stated drawback of the level set method, such as the coupled-
level-set and volume-of-fluid method (CLSVOF) by Sussmann and Puckett [16].
The CLSVOF method incorporates the mass conservation property of the VOF
method with the smoothness of the level set function [17]. Another hybrid
method is the hybrid particle level set method (HPLS) proposed by Enright et
al. [18]. This method updates the interface location computed using the
Eulerian transport equation through the use of Lagrangian markers, resulting
in improved mass conservation. Although all of these hybrid methods improved
the mass conservation property of the original level set method, they lack the
main benefits of the original level set method, i.e., the cost-effective,
straightforward implementation of the Euler transport equations using
different existing high-order schemes [12].

Furthermore, several studies have explored the use of mesh refinement
techniques to mitigate errors in mass conservation. For instance, Herrmann
[19] proposed the refined level set grid (RLSG) method, in which the level set
equation is solved on an auxiliary high-resolution grid. Another approach is
the standard arbitrary mesh refinement (AMR) method, in which the mesh is made
finer near the interface [20, 21]. Although mesh refinement techniques offer
improved mass conservation, they can be computationally expensive, difficult
to implement in parallel systems, and constrained by small time steps due to
the finer mesh resolution.

Studies by Olsson and Kreiss [22] and later Olsson et al. [23] addressed the
conservation issue of the classical level set method by proposing a
modification while maintaining its simplicity. They replaced the traditional
signed distance level set function with the diffuse interface profile defined
by the hyperbolic tangent function and solved the transport and re-
initialization equations in a conservative form. This approach showed an
improvement in mass conservation by an order of magnitude compared to results
using the signed distance function [22]. In this study, we will use the
conservative level set (CLS) approach to capture the interface between two
flows.

The implemented interface-capturing scheme should then be coupled with a
proper incompressible flow solver to simulate the physics of two-phase
incompressible flows. Here, we will utilize the projection method introduced
by Chorin [24] to model the behaviour of incompressible flows. In this
approach, the momentum equation is split into two parts. The first one solves
the momentum equation while ignoring the pressure term to calculate the
intermediate velocity field, which does not necessarily satisfy the
divergence-free constraint. Subsequently, the second equation uses pressure to
project the intermediate velocity field into a divergence-free velocity field.

Magnetic fields exert a considerable influence on the behaviour of conducting
fluids. The motion of these fluids in the presence of magnetic fields is
described through the coupling of the Navier–Stokes equations with Maxwell’s
equations of electromagnetics. This coupling gives rise to a set of equations
known as the magnetohydrodynamics (MHD) equations. The interaction between
electromagnetic fields and incompressible conducting fluids finds applications
in fusion reactors, the metallurgical industry, MHD generators, and aluminum
reduction cells [25]. To investigate the physics of these problems, a two-
phase MHD solver is required. While numerous numerical studies have been
conducted on simulating one-phase incompressible or compressible MHD problems
[26, 27, 28], the existing two-phase MHD solvers are highly limited due to
their complexity and multi-physics nature. For instance, Huang et al. [29]
developed a three-dimensional free-surface MHD solver to simulate the
evolution of the liquid lithium film free-surface due to the existing magnetic
forces in a fusion reactor known as NSTX (National Spherical Torus
Experiment). Their implemented MHD model is based on the magnetic field
induction equation, and the free boundary is tracked using the concept of a
fractional volume of fluid method. On the other hand, Gao et al. [30]
simulated the motion of a liquid lithium droplet under a strong non-uniform
magnetic field in a vacuum environment without the influence of gravity. In
their approach, the VOF model is incorporated to capture the interface, and
the CFS method is used to account for the surface tension. Later, Tagawa [31]
developed a numerical solver to investigate the movement of a falling droplet
of liquid metal into a pool of liquid metal under a uniform magnetic field in
the cylindrical geometry. The primary level set approach of Sussman et al.
[32] has been used in the study by Tagawa [31] to capture the interface
without the re-initialization step. In that study, the mass conservation,
convergence, and consistency of the solver have not been verified for other
benchmarks. Additionally, various studies have been also conducted to simulate
the two-phase MHD flows in the finite-element framework. For example, Yang et
al. [33] proposed a diffuse interface model to numerically simulate the two-
phase MHD flows and studied the performance of their solver for two-phase
Hartmann flows, which is the MHD version of the classical Poiseuille flows.

Despite various numerical efforts in MHD flows in compressible liquid, there
remains a need for a general and systematic implementation of a numerical
framework with a higher order of accuracy. This is crucial for effectively
capturing the formation of instabilities at the interface, which is of high
importance in different applications, particularly when dealing with abrupt
changes in magnetic properties across the interface. In this study, the
magnetostatic case of Maxwell’s equations is studied and integrated into the
implemented two-phase incompressible solver. Two-phase magnetostatic solvers
are widely employed to simulate the deformation of ferrofluid droplets in
various flow fields, with applications in different fields, including
biomedicine and rheology [34, 35, 36]. In the presence of a magnetic field in
two-phase magnetic flows, magnetic permeability experiences a discontinuity
across the interface, leading to the induction of the Lorentz force. This
force significantly influences the evolution of the interface. Therefore, it
is essential to incorporate the role of the Lorentz force into the governing
equations. The proposed two-phase magnetostatic solver adequately addresses
the magnetic permeability jumps across the interface, imposes proper boundary
conditions for the magnetic field at the interface, and satisfies the
divergence-free constraint for the magnetic field.

The objective of this study is two-fold. First, it primarily aims to present a
detailed second-order numerical toolkit for simulating the physics of two-
phase incompressible flows and its extension to magnetic flows. The
mathematical formulation and numerical grid implementation will be given in
Sec. 2. While investigating the surface instabilities necessitates the use of
higher-order numerical solvers, most numerical studies in this area are
limited to first-order accuracy. Thus, a fifth-order mass-conservation level
set approach is presented in Sec. 3, which includes a conservative re-
initialization step to minimize the mass loss. The accuracy and robustness of
the implemented level set solver are also analyzed using several benchmarks
existing in the literature. A high-order conservative incompressible solver
based on the projection method is discussed in Sec. 4 to solve the momentum
equation under the divergence-free constraint of the velocity field. The level
set and incompressible solvers are then coupled in Sec. 5, and the detailed
implementation of a second-order two-phase incompressible solver is
demonstrated. This numerical approach solves the governing equations in a
conservative form and effectively handles high density ratios and viscosity
jumps across the interface without introducing numerical instabilities.
Furthermore, a consistent method for calculating the interface curvature is
utilized, which is cost-effective and straightforward to implement. The
robustness and accuracy of this two-phase solver are examined through three
benchmarks. This effort is followed by introducing a two-phase
magnetohydrodynamics solver under the magnetostatic assumption, achieved by
extending the implemented second-order two-phase solver to the magnetic case.
To the knowledge of the authors, studies on developing a high-order two-phase
solver for magnetic flows are limited, and there is no specific study focusing
on the development of a two-phase magnetostatic solver in the finite-
difference framework using a high-order level set method to capture the
interface between two fluids. Thus, in Sec. 6, a procedure for adding the
magnetic terms to the solver is established. The introduced solver
successfully accounts for significant magnetic permeability variations across
the interface while satisfying the divergence-free condition of the magnetic
field. The performance of the solver is evaluated by proposing three test
cases: The deformation of a ferrofluid droplet in quiescent and shear flow
regimes and the magneto-Rayleigh–Taylor instability in magnetic fluids.

The second and principal contribution of this paper is the investigation of
sheared ferrofluid droplet deformation for various values of the droplet’s
susceptibility in both low and high capillary flow regimes. While previous
studies have examined the effects of various factors on droplet deformation,
such as the viscosity ratio between the droplet and the surrounding medium and
the strength and direction of the imposed magnetic field, exploring
deformations for different susceptibility values remains a crucial avenue to
explore. Therefore, in Sec. 6, this paper also investigates the impact of
droplet susceptibility values on its deformation, rotation, and potential
breakup.

##  2  Grid arrangement and mathematical formulation

The use of a staggered grid in an incompressible solver ensures the accurate
coupling between the velocity and pressure fields. The staggered arrangement,
as illustrated in Fig. 1, eliminates oscillations in the pressure field and
avoids the checker-board problem, a common issue for incompressible numerical
simulations [37]. In this computational grid system, scalar values such as
pressure are defined at the cell centers and velocity components are defined
at the cell faces. As a result, the continuity equation is solved at cell
center points, while the momentum equation corresponding to each velocity
component is defined at cell faces.

Figure 1:  Staggered grid system in Cartesian coordinates. In the staggered
grid arrangement, the values of scalar fields such as pressure ($p$), level
set function ($\psi$), density ($\rho$), dynamic viscosity ($\mu$), magnetic
permeability ($\mu_{\mathrm{m}}$), and $z-$component of the vector potential
($A_{z}$), are defined at cell centers. Velocity components, $u$ and $v$, as
well as magnetic field components, $B_{x}$ and $B_{y}$, are defined at cell
faces.

In this study, we have employed similar notation of the conservative centred
high-order finite-difference scheme of Morinishi et al. [37] and Desjardins et
al. [38], briefly introduced in this section for the sake of completeness.
According to their notation, the second-order finite-difference operator with
the stencil size $n$ for a variable $\phi$ in the $x_{1}-$direction in the
computational domain $\mathbf{x}=\left(x_{1},x_{2},x_{3}\right)$ is defined as

$\frac{\delta_{n}\,\phi}{\delta_{n}\,x_{1}}=\frac{\phi\left(x_{1}+n\,\Delta
x/2,\,x_{2},\,x_{3}\right)-\phi\left(x_{1}-n\,\Delta
x/2,\,x_{2},\,x_{3}\right)}{n\,\Delta x}.$ (1)

The second-order differential operator with respect to the $x_{2}-$ and
$x_{3}-$directions, $\frac{\delta_{n}\,\phi}{\delta_{n}\,x_{2}}$ and
$\frac{\delta_{n}\,\phi}{\delta_{n}\,x_{3}}$, can be defined in the same
manner. The second-order interpolation of a quantity $\phi$ defined on the
computational domain $\mathbf{x}=\left(x_{1},x_{2},x_{3}\right)$ with the
stencil size $n$ in the $x_{1}-$direction is given as

${\overline{\phi}}^{n\,x_{1}}=\frac{\phi\left(x_{1}+n\,\Delta
x/2,\,x_{2},\,x_{3}\right)+\phi\left(x_{1}-n\,\Delta
x/2,\,x_{2},\,x_{3}\right)}{2},$ (2)

and is defined similarly in the $x_{2}-$ and $x_{3}-$directions.

The $n^{\mathrm{th}}-$order central finite-difference operator in
$x_{i}-$direction is defined as

$\frac{\delta_{n{\mathrm{th}}}\,\phi}{\delta_{n{\mathrm{th}}}\,x_{i}}=\sum_{l=1}^{n/2}\alpha_{l}\frac{\delta_{\left(2l-1\right)}\,\phi}{\delta_{\left(2l-1\right)}\,x_{i}},$
(3)

where weight values, $\alpha_{l}$, are calculated as

$\sum_{l=1}^{n/2}\left(2l-1\right)^{2\left(k-1\right)}\alpha_{l}=\delta_{kl},\;\;k\in[1,\,n/2].$
(4)

Additionally, the $n^{\mathrm{th}}-$order interpolation in the
$x_{i}-$direction is given as

${\overline{\phi}}^{n{\mathrm{th}}\,\,x_{i}}=\sum_{l=1}^{n/2}\alpha_{l}\,{\overline{\phi}}^{\left(2l-1\right)\,x_{i}}.$
(5)

The introduced $n^{\mathrm{th}}-$order finite-difference and interpolation
schemes will be later used to discretize the governing equations.

##  3  Implementation of level set

The traditional level set function introduced by Sethian [15] is a smooth
signed distance function given as

$|\phi\left(\mathbf{x},t\right)|=|\mathbf{x}-\mathbf{x}_{\Gamma}|,$ (6)

where variable $\mathbf{x}_{\Gamma}$ denotes the closest point on the
interface from point $\mathbf{x}$. The positive and negative values of
$\phi(\mathbf{x},t)$ indicate the location of a point relative to the
interface, with convention determining which side is positive or negative.
Therefore, the zero iso-contour of the defined signed distance function,
$\phi(\mathbf{x},t)$=0, corresponds to the interface itself. The level set
motion under the velocity field $\mathbf{u}$ can be described using the
following transport equation

$\frac{\partial\phi}{\partial t}+\mathbf{u}\cdot\nabla\phi=0.$ (7)

Advecting the interface employing Eq. (7) for a few time steps can cause the
$\phi$ function to lose its signed distance property, becoming distorted and
losing its smoothness, thereby causing issues in the simulation [1]. To
prevent this issue, different re-initialization methods have been introduced
to reconstruct the $\phi$ function to be a smooth signed distance function
during the simulation. One of the well-known re-initialization techniques is
solving a Hamilton–Jacobi equation, introduced by Sussman et al. [32], which
can be solved using high-order numerical discretization schemes and recover
the distance function accurately. This method is proven to have various
limitations, thoroughly discussed in the literature [1, 15], which are outside
the scope of this paper. However, the main disadvantage of using this approach
for simulating two-phase flows is that both the transport equation and re-
initialization process fail to conserve the volume of the region enclosed by
the zero iso-contour. This can result in mass gain or loss in numerical
simulations, leading to unphysical results. In this study, the conservative
level set method (CLS) is employed to address this issue.

In the CLS approach [22], the interface between two immiscible flows is
defined using a diffuse profile in the form of a hyperbolic tangent function,
$\psi$, as below

$\psi(\mathbf{x},t_{0})=\frac{1}{2}\left(\tanh\left(\frac{\phi(\mathbf{x},t_{0})}{2\epsilon}\right)+1\right),$
(8)

where $\epsilon$ is a parameter to indicate the interface thickness and is
commonly defined as a function of the mesh resolution. In the hyperbolic
tangent definition, the iso-contour 0.5, $\psi\left(\mathbf{x},t\right)$=0.5,
specifies the location of the interface. For an incompressible flow,
$\nabla\cdot\mathbf{u}$ = 0, the transport equation can then be re-written as

$\frac{\partial\psi}{\partial t}+\nabla\cdot(\mathbf{u}\psi)=0.$ (9)

In order to recover the hyperbolic tangent form of the level set profile,
maintain the interface thickness, and prevent diffusion and smearing of the
interface during the simulation, a re-initialization step should be
introduced. The derived conservative re-initialization step by Olsson [22] is
given as

$\frac{\partial\psi}{\partial\tau}+\nabla\cdot\left(\psi\left(1-\psi\right)\mathbf{n}\right)=\nabla\cdot\left(\epsilon\nabla\psi\right),$
(10)

where the variable $\tau$ is the pseudo-time, $\mathbf{n}$ is the normal
vector at $\tau=0$, calculated as $\mathbf{n}=\nabla\psi/|\nabla\psi|$, and
the equation is solved until convergence is reached. In the proposed re-
initialization step, the compression flux,
$\psi\left(1-\psi\right)\mathbf{n}$, is included to maintain the resolution of
the interface and sharpen the profile which may smear due to the numerical
diffusion occurring during the simulation of the transport equation. Also, in
order to make sure that the level set profile remains of thickness $\epsilon$
and avoids the formation of discontinuities at the interface, the diffusion
flux, $\epsilon\nabla\psi$, with a small amount of viscosity is added to the
re-initialization equation. In Eq. (10), the diffusion in the normal direction
to the interface would be balanced by the compression term. However, diffusion
might also occur in the direction tangential to the interface, causing the
interface to move [23]. For that reason, in the later study by Olsson et al.
[23], the diffusion term has been modified as
$\nabla\cdot\left(\epsilon\left(\nabla\psi\cdot\mathbf{n}\right)\mathbf{n}\right)$
to avoid any tangential movement of the interface due to diffusion, improving
the re-initialization process.

###  3.1  High-order level set transport

Employing the notation introduced in Section 2, the $n^{\mathrm{th}}-$order
level set transport can be discretized as

$\nabla\cdot\left(\mathbf{u}\psi\right)=\sum_{i=1}^{3}\left(\frac{\delta_{2{\mathrm{nd}}}}{\delta_{2\mathrm{nd}}\,x_{i}}\left[u_{i}\overline{\psi}^{n\mathrm{th}\,\,x_{i}}\right]\right)$
(11)

where $\overline{\psi}^{n\mathrm{th}\,\,x_{i}}$ is a $n^{\mathrm{th}}-$order
interpolation of the variable $\psi$ to the cell face in $i-$direction and
$u_{i}$ is the $i^{\mathrm{th}}-$component of the velocity vector. Instead of
using central interpolation schemes, it is better to use an upwind-biased
scheme to prevent numerical oscillations appearing around $\psi=0.5$. Thus,
different upwind-biased approaches can be employed in order to calculate
interpolated $\psi$ values at the cell faces, such as WENO-type schemes [39,
40, 41] or High Order Upstream Central (HOUC) schemes [42]. The use of higher
order schemes improves the conservation of the transported level set, reducing
the need for re-initialization and yielding more accurate results. In this
study, we utilize the fifth-order WENO interpolation method, described in
detail in Appendix A. This approach is non-oscillatory and aims to mitigate
unwanted oscillations, especially near sharp gradients, providing a smooth
solution at each time step. The third-order total variation diminishing (TVD)
Runge–Kutta scheme is used for the temporal integration, presented in Appendix
B. Employing TVD spatial and temporal numerical schemes suppresses the
formation of unwanted oscillations in the numerical simulation of the $\psi$
profile, which is one of the main considerations in solving Eq. (9).

###  3.2  Level set conservative re-initialization

As mentioned earlier, it is essential to maintain a consistent thickness of
the hyperbolic tangent profile during the transport of the level set function.
However, numerical schemes may diffuse the interface, resulting in a violation
of mass conservation in the simulation. To address this issue, the re-
initialization step is incorporated into the level set transport equation. The
implemented re-initialization equation in this study, given as

$\frac{\partial\psi}{\partial\tau}+\nabla\cdot\left(\psi\left(1-\psi\right)\mathbf{n}\right)=\nabla\cdot\left(\epsilon\left(\nabla\psi\cdot\mathbf{n}\right)\mathbf{n}\right),$
(12)

is taken from the work by Olsson et al. [23], although their discretization
method is based on the finite-element approach. Therefore, we will establish a
robust, consistent finite-difference approach inspired by the study of
Desjardins et al. [1] to numerically discretize Eq. (12). To proceed, we
denote the diffusive and compressive fluxes as
$\mathbf{F}_{\mathrm{D}}=\epsilon\left(\nabla\psi\cdot\mathbf{n}\right)\mathbf{n}$
and $\mathbf{F}_{\mathrm{C}}=\psi\left(1-\psi\right)\mathbf{n}$, respectively.
According to the analysis shown by Desjardins et al. [1], employing a more
compact computational stencil to discretize Eq. (12) will lead to a more
accurate and robust reconstruction of the interface while eliminating the
appearance of spurious oscillations; hence, the second-order discretization is
employed. For the sake of clarity, we rewrite Eq. (12) for a two-dimensional
case as below

$\displaystyle\frac{\partial\psi}{\partial\tau}+\overbrace{\frac{\partial}{\partial
x}\left(\psi\left(1-\psi\right)n_{x}\right)+\frac{\partial}{\partial
y}\left(\psi\left(1-\psi\right)n_{y}\right)}^{\mathrm{compression}}=\underbrace{\epsilon\frac{\partial}{\partial
x}\left(\frac{\partial\psi}{\partial
x}{n_{x}}^{2}+\frac{\partial\psi}{\partial
y}n_{x}n_{y}\right)+\epsilon\frac{\partial}{\partial
y}\left(\frac{\partial\psi}{\partial x}n_{x}n_{y}+\frac{\partial\psi}{\partial
y}{n_{y}}^{2}\right)}_{\mathrm{diffusion}}.$ (13)

In order to update the cell center $\psi$ values, compression and diffusion
fluxes should be calculated at cell faces, and components of the normal vector
should be found at cell faces. To this end, normal values at cell faces, i.e.,
face normals [1], are calculated, determining $x-$ and $y-$components of the
normal vector at cell faces in both $x-$ and $y-$directions. The
$x_{i}-$component of $\psi$ gradient across the $x_{i}-$face is given as

$\left(\nabla\psi\right)_{i}^{x_{i}-\mathrm{face}}=\frac{\partial\psi}{\partial
x_{i}}\Bigg{|}_{x_{i}-\mathrm{face}}=\frac{\delta_{\mathrm{2nd}}\psi}{\delta_{\mathrm{2nd}}\,x_{i}},$
(14)

while for the $x_{j}-$component, a second-order interpolation of $\psi$ in the
$x_{i}-$direction is needed, and the face gradient value is calculated as

$\left(\nabla\psi\right)_{j}^{x_{i}-\mathrm{face}}=\frac{\partial\psi}{\partial
x_{j}}\Bigg{|}_{x_{i}-\mathrm{face}}=\frac{\delta_{\mathrm{2nd}}{\overline{\psi}}^{2\mathrm{nd}\,x_{i}}}{\delta_{\mathrm{2nd}}\,\,x_{j}}.$
(15)

Normalized face gradient values will correspond to normal vector values at
cell faces, defining the normal vector at $x_{i}-$face as
$\mathbf{n}^{x_{i}}=\left(\nabla\psi\right)^{x_{i}-\mathrm{face}}/|\left(\nabla\psi\right)^{x_{i}-\mathrm{face}}|$.
Therefore, the discrete version of diffusion and compression terms can be
written as

$\nabla\cdot\left(\epsilon\left(\nabla\psi\cdot\mathbf{n}\right)\mathbf{n}\right)=\epsilon\sum_{i=1}^{3}\left(\frac{\delta_{\mathrm{2nd}}}{\delta_{\mathrm{2nd}}\,x_{i}}\left[\mathbf{n}_{i}^{x_{i}}\left(\mathbf{n}^{x_{i}}\cdot\left(\nabla\psi\right)^{x_{i}-\mathrm{face}}\right)\right]\right),$
(16)

and

$\nabla\cdot\left(\psi\left(1-\psi\right)\mathbf{n}\right)=\sum_{i=1}^{3}\left(\frac{\delta_{\mathrm{2nd}}}{\delta_{\mathrm{2nd}}\,x_{i}}\left[\mathbf{n}_{i}^{x_{i}}{\overline{\psi\left(1-\psi\right)}}^{2\mathrm{nd}\,x_{i}}\right]\right),$
(17)

respectively. According to [23], the re-initialization equation converges
quickly. In their study, it has been shown that for the case of $\Delta
x\sim\Delta t\sim\epsilon$, the solution of the conservative re-initialization
will converge within one or two steps.

Desjardins et al. [1] demonstrated that by employing the length scale
analysis, the proper $\Delta\tau$ and the number of steps needed to obtain the
steady state solution of the re-initialization equation can be determined.
They showed that by taking the CFL number of the convection equation, Eq. (9),
to be $n$ times greater than the CFL number corresponding to the compression
term in the re-initialization equation, Eq. (12),
$n\,\mathrm{CFL}_{\mathrm{comp}}=\mathrm{CFL}_{\mathrm{conv}}$, the solution
of the re-initialization equation converges after $n$ steps. Thus, there is no
requirement to evaluate the convergence criteria during the simulation. Our
simulations utilize this approach and produce decent, robust results for the
re-initialization step, further explained in the following section.

###  3.3  Level set test cases

In Appendix C, the numerical order of accuracy of the implemented level set
solver is demonstrated through the rotating circle test case. In this section,
the accuracy, consistency, and robustness of the level set solver along with
the re-initialization procedure are verified and discussed using two
additional test cases, i.e., the circle in a deformation field and Zalesak’s
disk problems.

####  3.3.1 Circle in a deformation field

The main purpose of this test case is to evaluate the ability of the solver to
properly resolve thin filament structures, mainly appearing in stretching and
tearing flows [18]. Here, the initial center of the circle interface is
located at $(0.5,0.75)$, with the radius of $r=0.15$. The interface thickness
equals to $\epsilon=(\Delta x)^{0.7}/2$, and the simulation is performed on
the computational domain $[x,y]\in[0,1]\times[0,1]$, with the grid resolution
of $256\times 256$. The velocity field is defined as

$u=\sin\left(\pi x\right)\sin\left(\pi x\right)\sin\left(2\pi y\right),$ (18)

and

$v=-\sin\left(\pi y\right)\sin\left(\pi y\right)\sin\left(2\pi x\right),$ (19)

causing the circle interface to stretch out into a long, thin fluid element
that continuously wraps around itself. The CFL values are set as $0.5$ and
$0.25$ for transport and re-initialization equations, respectively, and the
re-initialization process is applied every five steps. Figure 2(a) displays
the evolution of the interface, $\psi=0.5$, at eight different time steps:
$t=0.5,1,1.5,2,2.5,3,3.5$, and 4 in black. Figure 2(a) visually confirms that
the implemented solver is capable of sustaining thin, elongated filament
structures of the interface for a long simulation time. The calculated maximum
area loss during the simulation, from $t=0$ to $4$, is equal to $2.8\%$,
demonstrating the area-conserving character of the level set solver, even
after stretching the vortex for a long time and causing the trailing ligament
thickness to become of the order of the mesh resolution. Additionally, in
order to evaluate how the thickness of the transition layer changes,
$\psi=0.05$ and $\psi=0.95$ contours are also depicted in blue and purple,
respectively, for the first six time steps, as for the final stages of the
vortex stretching, the thickness becomes too thin, and visualization of the
transition layer is a bit difficult. According to the presented contours, the
thickness of the transition layer remains almost constant even after drastic
changes in the interface. However, at the tail of the stretched vortex, it can
be seen that the transition layer is not constant. This is mainly due to the
fact that the thickness of the stretched circle becomes similar to the
thickness of the interface. This behavior was also observed in other studies
like the one by Olsson and Kreiss [22] and is called a pinch-off effect. The
pinch-off is a numerical effect and can be prevented if the interface’s
thickness is smaller than the distance between two interfaces.

Figure 2(b) presents the interface location of the vortex at $t=4.5$ in black,
while the blue iso-contour depicts the solution at the same time without
applying the re-initialization step during the simulation, from $t=0$ to
$4.5$. Figure 2(b) makes evident that incorporating the re-initialization
process in the simulation leads to a smoother solution with better area
conservation. The zoomed-in insets in Fig. 2(b) better illustrate that
employing the re-initialization step during the simulation results in a better
reconstruction of the thin filaments and prevent the tearing of the interface
in the solution, leading to a more conservative solution.

(a)

(b)

Figure 2:  (a) The interface location, $\psi=0.5$ (in black), of the single
vortex in the velocity field defined by Eqs. (18-19), from $t=0$ to $t=4.0$,
with the grid resolution of $256\times 256$ and CFL$=0.5$. For $t=0$ to $2.5$,
$\psi=0.05$ and $\psi=0.95$ contours are also plotted in blue and purple,
respectively. (b) Interface location of the vortex in the deformation field,
black: with re-initialization process applied every five time steps, and blue:
without re-initialization step.

To achieve a better quantitative analysis and error calculation, the vortex
field can be reversed in time by multiplying the velocity components by
$\cos\left(\pi t/T\right)$, where $T$ denotes the time that the circle returns
to its initial state, and, therefore, the interface location is known. In our
simulation, the periodicity is set to $T=2$, and the simulation is run from
$t=0$ to $t=4$, hence, the interface goes through two complete rotations.
Figure 3(a) illustrates the interface from $t=0$ to $4$ with increments of
$0.5$, showing that the circle interface goes through two complete rotations
and returns to its original place at $t=4$. The calculated root mean square
(rms) of the error for the level set field at $t=2$ and $4$, are $4.8\times
10^{-3}$ and $7\times 10^{-3}$, respectively.

The area loss during the simulation and the convergence of the numerical
result for increasing mesh resolutions are studied by repeating this test case
for four other resolutions, $16\times 16$, $32\times 32$, $64\times 64$, and
$128\times 128$. In order to investigate the mass conservation property of the
implemented level set, ${A}/{A}_{0}$ is calculated, where ${A}_{0}$ represents
the initial area and ${A}$ is the computed area at each time step, for all
five grid resolutions. In Fig. 3(b), the conservative behaviour of the
implemented level set is demonstrated by a decrease in area loss as the mesh
is refined and the value of ${A}/{A}_{0}$ becomes closer to one as the grid
resolution is increased. Figure 3(b) also compares the numerical and
analytical results for different mesh resolutions. Furthermore, increasing the
grid resolution leads to the convergence of the numerical result to its
analytical counterpart. The conservative property noticeably improved when the
mesh resolution is increased from $16\times 16$ to $32\times 32$. For mesh
resolutions $64\times 64$, $128\times 128$, and $256\times 256$, the
calculated $\psi=0.5$ contour is close to the analytical interface, and for
$128\times 128$ and $256\times 256$, the computed interface is
indistinguishable from the exact one. In order to evaluate the accuracy of the
solver, the $L_{\infty}$ error of the area loss during the simulation from
$t=0$ to $4$ is calculated for all five mesh resolutions. In Fig. 3(c), the
calculated error values are plotted against the computational grid size using
a logarithmic scale. According to Fig. 3(c), the error decreases as the mesh
is refined, and a second-order convergence is achieved. This convergence
analysis indicates that the solver is able to maintain second-order accuracy
for the mass conservation property even for more pronounced interface
deformations over prolonged periods of simulation. It is worth noting that the
complexity of the interface evolution may necessitate more frequent re-
initialization steps to ensure optimal mass conservation during the
simulation. The re-initialization step is second-order accurate and can affect
the global order of accuracy of the solver. That is why the convergence rate
here is expected to be less than that of the test case presented in Appendix
C.

(a)

(b)  (c)

Figure 3:  (a) The interface location of the circle in the vortex field with
the periodicity of $T=2$ at nine different time steps,
$t=0,0.5,1,1.5,2,2.5,3,3.5$, and $4$. The simulation is performed in the two-
dimensional computational domain, $[x,y]\in[0,1]\times[0,1]$, with the grid
resolution of $256\times 256$, and CFL$=0.5$. The re-initialization step is
applied every five time steps. (b) The temporal evolution of the normalized
area for the circle in the vortex field is compared for five different mesh
resolutions of $16\times 16$, $32\times 32$, $64\times 64$, $128\times 128$,
and $256\times 256$. Additionally, the interface location of the circle,
$\psi=0.5$, is shown for all five different grid resolutions. (c) Order of
convergence analysis of the area loss property. The $L_{\infty}$ error values
of the area loss are computed for five different grid resolutions and plotted
against the grid size.

####  3.3.2 Zalesak’s disk

The Zalesak’s disk problem, widely used in the literature to indicate the
robustness of a solver towards diffusion errors [43], is studied in this
section. In this test case, the ability of the solver to transport sharp
corners and thin structures without introducing noticeable diffusion errors
can be examined. This test case includes a notched circle in solid body
rotation under the velocity field of $u=2\pi\left(y-0.5\right)$ and
$v=-2\pi\left(x-0.5\right)$. Initially, a notched circle of radius $r=0.15$ is
centred at $(0.5,0.5)$ with a notch of width $0.05$ and height $0.25$ in the
computational domain of $(x,y)\in[0,1]\times[0,1]$. At $t=1$, the disk goes
through one complete rotation and returns to its primary location. The level
set thickness and CFL values are similar to the previous case, and a Dirichlet
boundary condition is imposed at all boundaries. Figure 4 displays the
solution at $t=1$ for three different grid resolutions $64\times 64$,
$128\times 128$, and $256\times 256$. As can be seen in Fig. 4, increasing the
grid resolution reduces the diffusion of the interface, and the final solution
is close to the initial description of the interface, recovering the edges and
corners of the disk more accurately (see zoomed-in insets). The calculated
area loss for the three grid resolutions from coarse to fine are $0.43\%$,
$0.25\%$, and $0.07\%$, respectively, showing the convergence of the solution
by increasing the grid resolution.

Figure 4:  The interface of the Zalesak’s disk, $\psi=0.5$, for three
different grid resolutions, $64\times 64,128\times 128$, and $256\times 256$,
from left to right. The black contour shows the initial interface while the
blue one represents the calculated interface at $t=1$. The numerical solution
converges to the exact solution by increasing the mesh resolution, accurately
capturing the corners and structure of the interface.

##  4 Implementation of incompressible Navier–Stokes solver

In this section, the implementation of a high-order incompressible solver is
presented. The density of each fluid particle does not change as it moves in
the incompressible flow regime. Therefore, the mass conservation equation
simplifies to a divergence-free condition for the velocity field, which must
be satisfied while solving the momentum equation. To this end, the projection
method is adopted in our algorithm, and two test cases are examined to verify
the accuracy and robustness of the implemented solver.

###  4.1  Governing equations

The incompressible form of the Navier–Stokes equations are given as

$\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)=-\frac{1}{\rho}\nabla
p+\frac{1}{\rho}\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)+\mathbf{g},$
(20a) $\nabla\cdot\mathbf{u}=0,$ (20b)

where $\mathbf{u}$ is the velocity vector, $\rho$ is the density field, and
$p$ is the pressure field. Variables $\mu$ and $\mathbf{g}$ denote the dynamic
viscosity and gravitational acceleration, respectively. It is noteworthy to
mention that the momentum equation is solved in a conservative form, ensuring
momentum conservation and avoiding unphysical numerical solutions. The
presented incompressible numerical solver employs a spatial staggered-variable
formulation, described in Section 2. Compared to a collocated formulation,
staggering offers the benefit of localized derivative stencils in space that
enhance the precision of the stencil. In order to solve the momentum equation,
the pressure gradient is required at cell faces, where the velocity values are
defined. In the staggered arrangement, the second-order pressure gradient can
be calculated using pressure values that are one cell apart. However, in the
collocated formulation, the pressure values that span three cells are needed
to compute the second-order pressure derivative. The localized stencils
accessible in staggered approaches are considerably more precise than the
broader stencils that would be used in a collocated approach on the same mesh
[44]. Furthermore, the staggered grid system offers a strong coupling between
the pressure and velocity field when a Poisson equation is solved for the
pressure compared to the collocated grids [37]. As a result, the unphysical
solution obtained in the collocated grid system due to weak coupling is
prevented.

###  4.2  Projection method

The main difficulty of numerically solving the incompressible Navier–Stokes
equations is the lack of an explicit time-derivative term in the continuity
equation. Consequently, the velocity divergence-free constraint must be
satisfied by implicitly coupling the mass equation with the pressure term in
the momentum equation [45]. To this end, this study incorporates the
projection method initially introduced by Chorin [24], which is based on the
Helmholtz–Hodge decomposition that states any vector field can be decomposed
into two components: solenoidal (divergence-free), and irrotational (curl-
free). In this operator splitting approach, the momentum equation is divided
into two separate equations, given as

$\frac{{\tilde{\mathbf{u}}}^{n+1}-\mathbf{u}^{n}}{\Delta
t}=-\nabla\cdot\left(\mathbf{u}{{}^{n}}\mathbf{u}{{}^{n}}\right)+\frac{1}{\rho}\nabla\cdot\mu\left(\nabla\mathbf{u}{{}^{n}}+(\nabla\mathbf{u}{{}^{n}})^{\mathrm{T}}\right)+\mathbf{g},$
(21a) and $\frac{\mathbf{u}^{n+1}-\tilde{\mathbf{u}}^{n+1}}{\Delta
t}+\frac{\nabla p^{n+1}}{\rho}=0,$ (21b)

where $\tilde{\mathbf{u}}$ denotes the intermediate velocity which does not
necessarily satisfy the divergence-free constraint. Adding Eq. (21a) and Eq.
(21b) recovers the original momentum equation, i.e., Eq. (20a). Equation
(21a), also known as the intermediate step, is straightforward to solve since
only one unknown is present in the equation, that is, the intermediate
velocity, $\tilde{\mathbf{u}}$. However, Eq. (21b) has two unknown variables
and cannot be solved in its present form. To address this issue, by taking the
divergence of Eq. (21b) and knowing that the velocity field solution at time
step $n+1$ should be divergence-free, $\nabla\cdot\mathbf{u}^{n+1}=0$, the
following equation will be obtained

$\nabla\left(\frac{1}{\rho}\nabla p^{n+1}\right)=\frac{1}{\Delta
t}\nabla\cdot\tilde{\mathbf{u}}^{n+1},$ (22)

which is a Poisson equation for pressure. Since the density is constant, Eq.
(22) can be simplified as

$\Delta p^{n+1}=\frac{\rho}{\Delta t}\nabla\cdot\tilde{\mathbf{u}}^{n+1}.$
(23)

In this way, the intermediate velocity field, which was initially computed
without forcing incompressibility, is projected into the divergence-free
field, satisfying Eq. (20b). Solving the Poisson equation for pressure results
in the solution for the pressure field. Lastly, the velocity field at time
step $n+1$ can be updated by rearranging Eq. (21b) to be read as

$\mathbf{u}^{n+1}=\tilde{\mathbf{u}}^{n+1}-\left(\frac{\Delta
t}{\rho}\right)\nabla p^{n+1}.$ (24)

Here, the projection method is explained using the first-order Euler scheme to
discretize time derivative terms. The projection method can be easily extended
to other numerical temporal-integration schemes, such as the third-order
Runge–Kutta method employed in this study.

Using the notation introduced in Section 2, the $n^{\mathrm{th}}-$order
spatial discretization of the diffusion term for the $x_{1}-$component can be
written as

$\displaystyle\left[\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)\right]_{x_{1}}^{n\mathrm{th-
order}}=\frac{\delta_{n\mathrm{th}}}{\delta_{n\mathrm{th}}\,x_{1}}\left[2\mu\left(\frac{\delta_{n\mathrm{th}}\,u_{1}}{\delta_{n\mathrm{th}}\,x_{1}}\right)\right]+\frac{\delta_{n\mathrm{th}}}{\delta_{n\mathrm{th}}\,x_{2}}\left[\overline{\overline{\mu}^{2\mathrm{nd}\,x_{1}}}^{2\mathrm{nd}\,x_{2}}\left(\frac{\delta_{n\mathrm{th}}\,u_{1}}{\delta_{n\mathrm{th}}\,x_{2}}+\frac{\delta_{n\mathrm{th}}\,u_{2}}{\delta_{n\mathrm{th}}\,x_{1}}\right)\right]$
$\displaystyle+\frac{\delta_{n\mathrm{th}}}{\delta_{n\mathrm{th}}\,x_{3}}\left[\overline{\overline{\mu}^{2\mathrm{nd}\,x_{1}}}^{2\mathrm{nd}\,x_{2}}\left(\frac{\delta_{n\mathrm{th}}\,u_{1}}{\delta_{n\mathrm{th}}\,x_{3}}+\frac{\delta_{n\mathrm{th}}\,u_{3}}{\delta_{n\mathrm{th}}\,x_{1}}\right)\right],$
(25)

where $\mathbf{u}=\left(u_{1},u_{2},u_{3}\right)$. In this study, the second-
order finite-difference scheme is implemented to discretize the diffusion
term. Hence, the $x_{1}-$ and $x_{2}-$components of the diffusion term,
$\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)$,
can be computed as

$\displaystyle\left[\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)\right]_{x_{1}}^{2\mathrm{nd-
order}}=\frac{\delta_{2}\,T_{x_{1}x_{1}}}{\delta_{2}\,x_{1}}+\frac{\delta_{2}\,T_{x_{1}x_{2}}}{\delta_{2}\,x_{2}},$
(26a) and
$\displaystyle\left[\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)\right]_{x_{2}}^{2\mathrm{nd-
order}}=\frac{\delta_{2}\,T_{x_{1}x_{2}}}{\delta_{2}\,x_{1}}+\frac{\delta_{2}\,T_{x_{2}x_{2}}}{\delta_{2}\,x_{2}},$
(26b) where
$T_{x_{1}x_{1}}=2\mu\left(\frac{\delta_{2}\,u_{1}}{\delta_{2}\,x_{1}}\right),$
(26c)
$T_{x_{1}x_{2}}=\mu\left(\frac{\delta_{2}\,u_{1}}{\delta_{2}\,x_{2}}+\frac{\delta_{2}\,u_{2}}{\delta_{2}\,x_{1}}\right),$
(26d) and
$T_{x_{2}x_{2}}=2\mu\left(\frac{\delta_{2}\,u_{2}}{\delta_{2}\,x_{2}}\right).$
(26e)

The discretization can be simply extended to the three-dimensional case,
$\mathbf{x}=(x_{1},x_{2},x_{3})$.

Most previous numerical studies have employed a second-order discretization
for the convection term, $\nabla\cdot(\mathbf{uu})$, in the context of two-
phase flows. This choice is made because discontinuities are expected to arise
in the velocity fields across the interface, and transitioning to a higher-
order finite difference scheme introduces numerical challenges. In this study,
instead of utilizing second-order interpolation, we adopted fifth-order WENO
interpolation for discretizing the convection term. This approach has been
demonstrated to yield approximately third-order accuracy for the single-phase
incompressible solver (please refer to Appendix D) and provides second-order
accuracy for the two-phase solver. This improvement is notable when compared
to other second-order based numerical discretizations. The employed
discretization of $\nabla\cdot(\mathbf{uu})$ term for
$\mathbf{x}=(x_{1},x_{2})$ in this study is given as

$\displaystyle\left[\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)\right]_{x_{1}}=\frac{\delta_{2}}{\delta_{2}\,x_{1}}\left(\overline{u_{1}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{1}}\,\,\overline{u_{1}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{1}}\right)+\frac{\delta_{2}}{\delta_{2}\,x_{2}}\left(\overline{u_{1}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{2}}\,\,\overline{u_{2}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{1}}\right),$
(27a) and
$\displaystyle\left[\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)\right]_{x_{2}}=\frac{\delta_{2}}{\delta_{2}\,x_{1}}\left(\overline{u_{1}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{2}}\,\,\overline{u_{2}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{1}}\right)+\frac{\delta_{2}}{\delta_{2}\,x_{2}}\left(\overline{u_{2}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{2}}\,\,\overline{u_{2}}^{5^{\mathrm{th}}\mathrm{WENO}\,\,x_{2}}\right).$
(27b)

The Poisson equation for the pressure is discretized using the second-order
central finite-difference scheme. The solution of the Poisson equation is
computed at each iteration of the simulation by employing Krylov methods and
multi-grid preconditioning implemented in the PETSc package library [46].

In the staggered grid arrangement shown in Fig. 1, boundary conditions for the
velocity components normal to the boundary can be readily implemented.
However, imposing the proper boundary condition for the tangential components
of the velocity can be more challenging. Thus, ghost points are used to apply
the desired boundary condition for velocity components tangential to the
boundary [47]. For example, consider a case where the no-slip boundary
condition is needed to be implemented for the velocity component in the
$y-$direction, $v$, at the left boundary (see Fig. 5). In order to update the
value of $v_{i,j+1/2}$, which has a half-cell offset from the left boundary,
and to impose a proper boundary, the value of $v_{i-1,j+1/2}$ is needed that
is located outside of the domain. The tangential velocity at this ghost point
can be obtained by using linear interpolation as
$v_{\mathrm{boundary}}=\left(v_{i,j+1/2}+v_{i-1,j+1/2}\right)/2$, where
$v_{\mathrm{boundary}}$ is the velocity value at the boundary, equal to the
wall tangential velocity for the no-slip case. Therefore, the value at the
ghost point can be calculated as
$v_{i-1,j+1/2}=2v_{\mathrm{boundary}}-v_{i,j+1/2}$, and the no-slip boundary
condition is imposed. The same approach can be taken for other boundary
conditions as well. For instance, for the slip boundary condition where the
derivative of the tangential velocity is zero, the value of the ghost point is
given as $v_{i,j+1/2}=v_{i-1,j+1/2}$. In Appendix D, the order of accuracy and
performance of the introduced incompressible solver are examined. Interested
readers can refer to this section for more details.

Figure 5:  The proper boundary condition for the tangential velocity is
implemented by employing ghost points located outside of the computational
domain, shown in blue. The tangential velocity of these ghost points is
specified to enable the linear interpolation to produce the intended tangent
velocity at the wall.

##  5  Implementation of two-phase solver for nonmagnetic flows

This section presents the methodology implemented for simulating
incompressible two-phase flows. The introduced solver is based on coupling the
conservative level set (CLS) approach, detailed in Section 3 with the
incompressible solver presented in Section 4. In two-phase flows, material
properties such as density and viscosity experience a jump across the
interface, requiring special numerical treatments to avoid the appearance of
numerical instabilities. In addition, proper boundary conditions must be
imposed across the interface to obtain accurate physical results. This
section, thus, discusses these issues and introduces a complete solution
procedure for modelling two-phase flows.

###  5.1  One-fluid formulation

The governing equation for the one-fluid formulation approach to describe the
two-phase incompressible flows is given as

$\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)=-\frac{1}{\rho}\nabla
p+\frac{1}{\rho}\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)+\mathbf{g}+\frac{1}{\rho}\mathbf{F}_{\mathrm{v}},$
(28)

where $\mathbf{F}_{\mathrm{v}}$ represents any volume force that may be
present. The solution of Eq. (28) should also satisfy the incompressibility
constraint. In each phase, the material properties are constant, that is,
$\rho=\rho_{\mathrm{l}}$ and $\mu=\mu_{\mathrm{l}}$ for the liquid phase,
while $\rho=\rho_{\mathrm{g}}$ and $\mu=\mu_{\mathrm{g}}$ for the gas phase.
However, across the thin interface, denoted by $\Gamma$, fluid properties
experience a jump that can be written as
$[\rho]_{\Gamma}=\rho_{\mathrm{l}}-\rho_{\mathrm{g}}$ and
$[\mu]_{\Gamma}=\mu_{\mathrm{l}}-\mu_{\mathrm{g}}$ for the density and dynamic
viscosity, respectively. Since there is no mass exchange between the two
phases, the normal component of the velocity should be continuous across the
interface, i.e.,

$[\mathbf{u}\cdot\mathbf{n}]_{\Gamma}=\mathbf{u}_{\mathrm{l}}\cdot\mathbf{n}-\mathbf{u}_{\mathrm{g}}\cdot\mathbf{n}=0,$
(29)

where $\mathbf{n}$ is the normal vector to the interface. For viscous flows,
the tangential component of the velocity should also be equal for the two
phases at the interface. Thus, the velocity field should be continuous across
the interface, and the boundary condition for the velocity can be written as

$[\mathbf{u}]_{\Gamma}=0.$ (30)

Additionally, applying the momentum conservation principle to a control volume
located at the interface leads to the following boundary condition for the
pressure jump

$[\left(-p\mathbf{I}+\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)\right)\cdot\mathbf{n}]_{\Gamma}=0,$
(31)

where $\mathbf{I}$ is the identity tensor. If the surface tension force is
considered, the boundary condition for the pressure jump across the interface
is modified as

$[\left(-p\mathbf{I}+\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)\right)\cdot\mathbf{n}]_{\Gamma}=\sigma\kappa\mathbf{n},$
(32)

where variable $\sigma$ denotes the surface tension coefficient. The curvature
of the interface, $\kappa$, is calculated as
$\kappa=-\nabla\cdot\mathbf{n}=-\nabla\cdot\left(\nabla\psi/|\nabla\psi|\right)$.
In the present study, the corresponding jump condition in the pressure
gradient is modelled by employing the continuum surface force (CSF) method of
Brackbill et al. [11]. The CSF approach defines surface tension as a volume
force spreading over the finite interface width, expressed as
$F_{\sigma}=\sigma\kappa\nabla\psi$. Therefore, proper calculation of the
interface curvature is essential for accurate and robust surface tension
modelling, which will be discussed in more detail in the next section. The
viscous term is also discretized using the CSF method.

###  5.2  Projection method and discretization

As previously introduced, the projection method is employed to solve the
incompressible Navier–Stokes equations. This method involves two steps, known
as prediction and correction. During the prediction step, Eq. (28) is solved,
ignoring the pressure term, to advance the velocity field $\mathbf{u}^{n}$ to
an intermediate velocity $\tilde{\mathbf{u}}^{n+1}$. In the correction step,
the intermediate velocity is projected to its divergence-free solution,
$\mathbf{u}^{n+1}$, using the solution of the pressure Poisson equation. The
projection method has been discussed in Sec. 4.2 for the single incompressible
flow case. However, for two-phase flows, special considerations and numerical
treatments must be taken into account. The discretization of the convective
term, $\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)$, can be performed
similarly to the one demonstrated earlier in Sec. 4.2. The main challenges
associated with treating the interface discontinuities in two-phase
incompressible flows include properly modelling the viscosity discontinuity
across the interface and accurately calculating the curvature. The viscosity
term in the momentum equation must be appropriately discretized, particularly
in the presence of a dynamic viscosity jump, to ensure the accurate
calculation of kinetic energy dissipation near the interface. Additionally, a
robust and accurate method to assess the evolution of the interface curvature
should be employed to model the surface tension force. Another issue in two-
phase flow simulations pertains to the numerical discretization of the
pressure gradient term. The pressure gradient component in Eq. (28) includes
the density term as well, and since the density has a jump across the
interface, specific consideration is required. We will conclude this section
by discussing these issues in more detail and presenting the discretization
used to model viscosity, surface tension, and pressure jump terms.

#### Viscosity and surface tension modelling

The viscosity term,
$\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)$,
can be discretized using Eq. (4.2), which displays the general form of the
$n^{\mathrm{th}}-$order spatial discretization of this term, considering a
non-constant dynamic viscosity field. The implemented two-phase solver employs
the second-order finite-difference scheme and second-order linear
interpolation to discretize this term. Furthermore, the density and viscosity
are assumed to smoothly vary over the interface [22]. Therefore, the density
and viscosity fields can be represented using the level set function as

$\rho\left(\mathbf{x},t\right)=\rho_{\mathrm{g}}+(\rho_{\mathrm{l}}-\rho_{\mathrm{g}})\psi\left(\mathbf{x},t\right)$
(33)

and

$\mu\left(\mathbf{x},t\right)=\mu_{\mathrm{g}}+(\mu_{\mathrm{l}}-\mu_{\mathrm{g}})\psi\left(\mathbf{x},t\right),$
(34)

respectively. As a result, the discontinuity in the density and viscosity
fields is smoothed out across the interface within a few layers of cells,
which is the function of the interface thickness. Thus, numerical
instabilities that may appear due to the sharp jump in fluid properties are
avoided in the solution.

According to the CSF model introduced by Brackbill et al. [11], the surface
tension force per unit interfacial area between two fluids with a constant
surface tension coefficient, $\sigma$, can be written as

$\mathbf{F}_{\mathrm{sa}}\left(\mathbf{x}_{\Gamma}\right)=\sigma\kappa\left(\mathbf{x}_{\Gamma}\right)\mathbf{n}\left(\mathbf{x}_{\Gamma}\right),$
(35)

where $\mathbf{x}_{\Gamma}$ is an arbitrary point located on the interface.
The introduced interfacial surface tension force can be recast as

$\mathbf{F}_{\mathrm{sv}}(\mathbf{x})=\sigma\overbrace{\left(-\nabla\cdot\frac{\nabla\psi}{|\nabla\psi|}\right)}^{\kappa}\nabla\psi,$
(36)

which represents the surface tension volume force at any point $\mathbf{x}$
[22]. The introduced force $\mathbf{F}_{\mathrm{sv}}(\mathbf{x})$ results in
the same total force as that of
$\mathbf{F}_{\mathrm{sa}}(\mathbf{x}_{\Gamma})$, but it is distributed over
the width of the interface. This approximation is only valid for small
interface thicknesses. Therefore, it is essential to avoid excessively wide
interfaces and to maintain a constant thickness for the interface during the
simulation, which has been addressed in the implementation of the level set
(see Sec. 3.2). In order to calculate the surface tension force, the solver
should be able to robustly and accurately compute the curvature value.
Employing high-order numerical schemes to calculate curvature results in
oscillations appearing in the curvature field, leading to an unphysical
solution for the velocity field, known as spurious currents. To tackle this
issue, first- or second-order schemes are usually used to calculate curvature
values. Various methods have been developed, such as utilizing height
functions or implementing the least-squares minimization approach [48], to
formulate a robust and consistent framework for curvature calculation.
However, most of these approaches are computationally expensive and difficult
to implement. In this study, curvature computation is based on using the
calculated face normals introduced in Sec. 3.2. Therefore, the curvature field
for the two-dimensional case is given as

$\kappa=\frac{\delta_{2\mathrm{nd}}\,(\mathbf{n}_{x})^{x-\mathrm{face}}}{\delta_{2\mathrm{nd}}\,x}+\frac{\delta_{2\mathrm{nd}}\,(\mathbf{n}_{y})^{y-\mathrm{face}}}{\delta_{2\mathrm{nd}}\,y}.$
(37)

Equation (37) is second-order accurate and can be easily implemented,
requiring no additional computational effort. Various test cases have been
conducted to evaluate the results obtained from the Eq. (37) curvature
calculation, demonstrating its robustness. Notably, the results exhibit no
discernible oscillations, which will be expounded upon in the subsequent
section. Finally, for the two-dimensional case of $\mathbf{x}=(x,y)$, the $x-$
and $y-$components of the surface tension force spatial discretization can be
written as

$[\sigma\kappa\nabla\psi]_{x-\mathrm{comp}}=\sigma\overline{\left(\frac{\delta_{2\mathrm{nd}}\,(\mathbf{n}_{x})^{x-\mathrm{face}}}{\delta_{2\mathrm{nd}}\,x}+\frac{\delta_{2\mathrm{nd}}\,(\mathbf{n}_{y})^{y-\mathrm{face}}}{\delta_{2\mathrm{nd}}\,y}\right)}^{2\mathrm{nd}\,x}\frac{\delta_{2\mathrm{nd}}\,\psi}{\delta_{2\mathrm{nd}}\,x}$
(38a) and
$[\sigma\kappa\,\nabla\psi]_{y-\mathrm{comp}}=\sigma\overline{\left(\frac{\delta_{2\mathrm{nd}}\,(\mathbf{n}_{x})^{x-\mathrm{face}}}{\delta_{2\mathrm{nd}}\,x}+\frac{\delta_{2\mathrm{nd}}\,(\mathbf{n}_{y})^{y-\mathrm{face}}}{\delta_{2\mathrm{nd}}\,y}\right)}^{2\mathrm{nd}\,y}\frac{\delta_{2\mathrm{nd}}\,\psi}{\delta_{2\mathrm{nd}}\,y},$
(38b)

respectively.

#### Poisson equation

In the two-phase incompressible solver, since the density field is not
constant in the computational domain, the Poisson equation cannot be
simplified as Eq. (23), and the variable coefficients should be considered
while discretizing the Poisson equation. Therefore, the second-order
discretization of the Poisson equation can be written as

$\displaystyle\nabla\left(\frac{1}{\rho}\nabla
p\right)\Bigg{|}_{i,j}=\frac{\partial}{\partial
x}\left(\frac{1}{\rho}\frac{\partial p}{\partial
x}\right)\Bigg{|}_{i,j}+\frac{\partial}{\partial
y}\left(\frac{1}{\rho}\frac{\partial p}{\partial
y}\right)\Bigg{|}_{i,j}=\frac{\left(\frac{1}{\rho}\frac{\partial p}{\partial
x}\right)_{i+1/2,j}-\left(\frac{1}{\rho}\frac{\partial p}{\partial
x}\right)_{i-1/2,j}}{\Delta x}+\frac{\left(\frac{1}{\rho}\frac{\partial
p}{\partial y}\right)_{i,j+1/2}-\left(\frac{1}{\rho}\frac{\partial p}{\partial
y}\right)_{i,j-1/2}}{\Delta y}$
$\displaystyle=\frac{\left(\frac{2}{\rho_{i+1,j}+\,\rho_{i,j}}\frac{p_{i+1,j}-p_{i,j}}{\Delta
x}\right)-\left(\frac{2}{\rho_{i,j}+\,\rho_{i-1,j}}\frac{p_{i,j}-p_{i-1,j}}{\Delta
x}\right)}{\Delta
x}+\frac{\left(\frac{2}{\rho_{i,j+1}+\,\rho_{i,j}}\frac{p_{i,j+1}-p_{i,j}}{\Delta
y}\right)-\left(\frac{2}{\rho_{i,j}+\,\rho_{i,j-1}}\frac{p_{i,j}-p_{i,j-1}}{\Delta
y}\right)}{\Delta y}.$ (39)

It is noteworthy to mention that second-order interpolation is used to
calculate density values at cell faces, e.g.,
$(1/\rho)_{i+1/2,j}=2/\left(\rho_{i+1,j}+\rho_{i,j}\right)$. The same
interpolation should also be applied while calculating the $1/\rho$
coefficients for pressure gradient, viscosity, and surface tension terms in
Eq. (28). Furthermore, it is imperative to use the same finite-difference
scheme when computing the gradient of pressure and level set gradient for
determining the surface tension force. This ensures that the gradient operator
provides the proper force balance between pressure gradient and surface
tension across the interface, thereby reducing the formation of unphysical
spurious velocities in the solution [48].

###  5.3  Solution procedure

The complete solution procedure for the two-phase incompressible solver can be
summarized as follows:

  1. 1.

The conservative level set (CLS) approach is used to implicitly advance the
location of the interface from $t^{n}$ to $t^{n+1}$, employing the velocity
field at $t^{n}$.

  2. 2.

The updated location of the interface, $\psi^{n+1}$, is then utilized to
obtain the density and viscosity fields at $t^{n+1}$. Thus, $\rho^{n+1}$ and
$\mu^{n+1}$ fields are computed using Eqs. (33) and (34).

  3. 3.

The intermediate velocity field at time step $t^{n+1}$ is calculated by
solving Eq. (28), while ignoring the pressure term.

  4. 4.

The obtained intermediate velocity is then projected into a divergence-free
field by solving the Poisson equation for the pressure based on the
discretization introduced in Eq. (5.2).

  5. 5.

The correct velocity field is calculated at time $n+1$ by solving Eq. (24).
The obtained velocity at $t^{n+1}$ is then employed to advance the level set
profile for the next time step.

Based on the CFL condition, the stability constraint for the time step due to
convection, viscosity, and surface tension terms is given as [48, 49]

$\Delta t\leq\mathrm{min}\left(\frac{\Delta
x}{\mathrm{max}(||\mathbf{u}||)},\sqrt{\frac{{\Delta
x}^{3}\left(\rho_{\mathrm{l}}+\rho_{\mathrm{g}}\right)}{4\pi\sigma}},\frac{1}{4}\frac{{\Delta
x}^{2}}{\mathrm{max}(\nu_{\mathrm{l}},\nu_{\mathrm{g}})}\right).$ (40)

Usually, the surface tension limits the time step; however, for large density
ratios, the viscosity term may be more restrictive.

###  5.4  Two-phase solver test cases

To investigate the accuracy, robustness, and performance of the implemented
two-phase solver, three test cases are studied: the static droplet, damping
surface wave, and Rayleigh–Taylor instability. The static droplet test case,
which aims to evaluate the capability of the two-phase solver to accurately
calculate curvature and model surface tension forces, is presented in Appendix
E, and interested readers are encouraged to refer to this section for more
detail. The latter two test cases are examined in this section. The damping
surface wave benchmark is employed to investigate the order of accuracy of the
implemented solver and the interaction between viscosity and surface tension
terms in a simulation. Lastly, the Rayleigh–Taylor instability is investigated
to assess the performance of the solver in handling the evolution of a complex
interface and also in treating cases where a high-density ratio exists across
the interface.

####  5.4.1 Damping surface wave

The viscous damping of a surface wave is a well-known test case in the
literature, widely used to assess the capability of the implemented solver to
accurately model the interaction between viscosity and surface tension forces.
In this test case, two superimposed fluids with density $\rho_{1}$ and
$\rho_{2}$ are separated by a sinusoidal interface with wavelength $\lambda$
and amplitude ${A}_{0}$, and thus the perturbed interface profile is given as

$y=y_{0}-{A}_{0}\cos(2\pi x/\lambda).$ (41)

For the case where both fluids have the same kinematic viscosity, $\nu$, and
constant surface tension, $\sigma$, the analytical solution for the evolution
of the wave amplitude with time is derived by Prosperetti [50] by employing
the initial value theorem. The initial value solution is obtained as [50]

$\displaystyle{A}_{\mathrm{exact}}(t)=\frac{4(1-4\beta)k^{4}\nu^{2}}{8(1-4\beta)k^{4}\nu^{2}+{\omega_{0}}^{2}}\,{A}_{0}\,\mathrm{erfc}\left(\sqrt{\nu
k^{2}t}\right)+\sum\limits_{i=1}^{4}\frac{z_{i}}{Z_{i}}\left(\frac{{\omega_{0}}^{2}\,A_{0}}{{z_{i}}^{2}-\nu
k^{2}}\right)\exp[\left({z_{i}}^{2}-\nu
k^{2}\right)t]\mathrm{erfc}\left(z_{i}\sqrt{t}\right),$ (42)

where $z_{i}$ are the roots of

$\displaystyle
z^{4}-4\beta\sqrt{k^{2}\nu}\,z^{3}+2\left(1-6\beta\right)k^{2}\nu\,z^{2}+4(1-3\beta)(k^{2}\nu)^{\frac{3}{2}}\,z+(1-4\beta)k^{4}\nu^{2}+{\omega_{0}}^{2}=0.$
(43)

The dimensionless parameter $\beta$ is defined as
$\beta=\rho_{1}\rho_{2}/\left(\rho_{1}+\rho_{2}\right)^{2}$, while the
inviscid frequency of the wave oscillation is given by
$\omega_{0}=\sqrt{\frac{\sigma k^{3}}{\rho_{1}+\rho_{2}}}$, and
$Z_{i}=\Pi_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{4}\left(z_{j}-z_{i}\right)$.

Here, the simulation is performed in a computational domain of
$[0,1]\times[0,1]$, with periodic boundary conditions along the $x-$direction
and slip wall boundary along the $y-$direction, and $y_{0}$ is set to
$y_{0}=0.5$. The wavelength of the perturbation is set to $\lambda=1$, and the
initial amplitude of the wave is $A_{0}=0.01\lambda$. Three scenarios are
investigated, namely, density ratios of 1, $10$, and $1000$, for four
different grid resolutions of $8\times 8$, $16\times 16$, $32\times 32$, and
$64\times 64$, with the interface thickness being set to $\epsilon=\Delta
x/2$. For the first case of unity density ratio, the surface tension
coefficient is $\sigma=2$, with the constant kinematic viscosity for both
fluids being set to $\nu=0.0647$, and $\rho_{1}=\rho_{2}=1$. The simulation
was run for all different mesh resolutions with the constant time step $\Delta
t=0.0005$. Figure 6[top-left] displays the time evolution of the normalized
wave amplitude, i.e., ${A}/\lambda$, for all four meshes, along with the
analytical solution. This figure visually confirms that, as mesh resolution is
increased, the numerical solution converges to the analytical one. For a
better quantitative comparison, the rms value of the error in the wave
amplitude is plotted against the mesh resolution in Fig. 6[top-right].

As can be concluded from the computed error values, although for the coarsest
mesh resolution, the solver is able to capture the proper behaviour of the
surface wave, a small error in the numerical oscillation period of the wave
has caused a noticeable error value for this mesh. As the mesh resolution is
increased, more accurate results are obtained, and according to Fig. 6[top-
right], close to a second-order convergence rate is observed. It is noteworthy
to mention that for the case of $\rho_{2}/\rho_{1}=1$, there is no dynamic
viscosity or momentum jump across the interface. Therefore, the numerical
errors solely originate from the solution of the level set transport equation
and the curvature computation. However, by increasing the density ratio, due
to the existence of the density and dynamic viscosity jump across the
interface, both convection and viscosity terms of the momentum equation affect
the solution. Therefore, the numerical solver is also evaluated for density
ratios of 10 and 100 to investigate whether the numerical model is capable of
addressing these jump conditions and momentum transfer across the interface in
the presence of a surface tension force. For both density ratios of $10$ and
$100$, the kinematic viscosity of the fluids is set to $\nu=0.00647$, and time
steps are set to $\Delta t=0.001$ and $0.003$, respectively. As expected, as
the mesh resolution increases, the numerical results become more accurate for
both cases, as shown in Fig. 6[middle-left, bottom-left]. According to the
convergence study, the approximate convergence rate of 1.6 is achieved for
both cases, see Fig. 6[middle-right, bottom-right], demonstrating that in the
presence of a density and viscosity jump, the solver tends to preserve its
order of accuracy.

Figure 6:  The solutions of the damped surface wave problem for three density
ratios of [top] unity, [middle] 10, and [bottom] 100, are represented for four
different mesh resolutions of $8\times 8$, $16\times 16$, $32\times 32$, and
$64\times 64$. The accuracy analysis of the implemented two-phase solver is
also illustrated for all three density ratios.

####  5.4.2 Rayleigh–Taylor instability

The Rayleigh–Taylor instability occurs when a layer of liquid is superimposed
to another less dense liquid layer in such a way that by interchanging the
fluids, the energy of the system can be reduced. The Rayleigh–Taylor
instability has been widely studied in classical hydrodynamics and is studied
as a benchmark to evaluate the performance of the implemented two-phase solver
for simulating problems containing the highly nonlinear and multi-scale nature
of fluid dynamics. To simulate the Rayleigh–Taylor instability, we consider a
two-dimensional rectangular domain $[x,y]\in[0,1]\times[0,4]$, with a fluid
interface parallel to the horizontal axis at $y_{0}=2.0$. The fluid interface
is initialized with a small sinusoidal perturbation, whose wavelength and
amplitude are $2\pi$ and 0.1, respectively. Following the study by Huang et
al. [51], the material properties are set to $\rho_{1}=1$, $\rho_{2}=3$,
$\mu_{1}=\mu_{2}=0.01$, $\sigma=10^{-12}$, and the gravity is acting downwards
with a magnitude of unity. The simulation is performed for four mesh
resolutions of $16\times 64$, $32\times 128$, $64\times 256$, and $128\times
512$. For all cases, the time step is set to $\Delta t=5\times
10^{-4}/\sqrt{\mathrm{At}}$, where the Atwood number is defined as
$\mathrm{At}=\left(\rho_{2}-\rho_{1}\right)/\left(\rho_{1}+\rho_{2}\right)$.

Figure 7(a) depicts the results for the four different mesh resolutions at
time $t\sqrt{\mathrm{At}}=0$ to 2.5 with an increment of $0.25$. It is
visually evident that even for the coarsest mesh resolution case, the main
features of the Rayleigh–Taylor instability are captured. However, refining
the mesh allows for a better representation of the evolution and growth of
bubbles and spikes in the solution. For a more comprehensive quantitative
analysis, we validated our results with four different previous studies by
Huang et al. [51], Ding et al. [52], Guermond and Quartapelle [53], and
Tryggvason [47]. To this end, the transient location of the spike tip and the
interface location at the left (right) edge during the simulation are compared
to results from previous studies, shown in Fig. 7(b). As can be observed from
Fig. 7(b), by increasing the mesh resolution, our results converge to ones
from earlier studies, especially for the mesh resolution of $128\times 512$,
the results closely match those reported in the literature. Note that the
study performed by Tryggvason [47] only considered the inviscid case, which
explains the slight deviation from the other results. Assuming the solution of
the finest grid, $128\times 512$, as an analytical solution, the convergence
rate of the solution is calculated. To this end, the $L_{2}$ norm of the spike
and bubble locations during the simulation is computed. Figure 7(c) displays
the order of accuracy, which is around 1.5. The obtained convergence rate is
close to the expected second-order, confirming that the solver is robust and
accurate for more complex problems as well.

(a)

(b)  (c)

Figure 7:  (a) The interface of the Rayleigh–Taylor instability with a density
ratio of 3 ($\mathrm{At}=0.5$) at
$t\sqrt{\mathrm{At}}=0,1,1.25,1.5,1.75,2,2.25$, and 2.5, from left to right,
for four different mesh resolutions. (b) The transient location of the spike
tip and bubble is depicted during the simulation for four different mesh
resolutions of $16\times 64$, $32\times 128$, $64\times 256$, and $128\times
512$. The results are compared to four different studies by Huang et al. [51],
Ding et al. [52], Guermond and Quartapelle [53], and Tryggvason [47]. (c)
$L_{2}$ error in spike tip and bubble locations during the Rayleigh–Taylor
simulation with the density ratio of three for the mesh resolutions of
$16\times 64$, $32\times 128$, and $64\times 256$.

Finally, the performance of the solver is also evaluated by simulating the
Rayleigh–Taylor instability for higher density ratios of 30, 1000, and 3000.
The mesh resolution is set as $64\times 256$ and $\Delta t=5\times
10^{-4}/\sqrt{\mathrm{At}}$. Figure 8 displays the results for density ratios
of 30 ($\mathrm{At}=0.935$), 1000 ($\mathrm{At}=0.998$), and 3000
($\mathrm{At}=0.999$), at time $t\sqrt{\mathrm{At}}=0$ to 2 with an increment
of 0.25. As expected, the interface evolves faster for higher density ratios,
and the interface structures are simpler. Therefore, for the minimum density
ratio of 30, the rate of penetration of the heavier fluid to the lighter one
is less as a result of the reduced growth rate of the Rayleigh–Taylor
instability, and the mushroom structure of the Rayleigh–Taylor instability can
be observed. Since the Atwood number is comparable for two cases of density
ratios of $1000$ and $3000$, the results of the surface evolution are similar
(see Fig. 8). However, in the case of the density ratio of $3000$, the tip of
the spike has penetrated the lower density region a bit further, and its
structure is sharper at the tip (see the insets of Fig. 8). The greater
density ratio of 3000 indicates a more significant disparity in fluid
densities, leading to stronger buoyancy forces. This stronger buoyancy causes
the spikes at the interface to penetrate farther into the lower-density region
due to the accelerated downward motion of the denser fluid. Additionally, a
greater density ratio corresponds to a higher Atwood number, which accelerates
the exponential growth of Rayleigh–Taylor instability, further promoting the
deeper penetration of spike tips into the low-density region. However, because
the Atwood numbers for these cases are relatively close, the difference in
results remains relatively small.

Figure 8: The interface of the Rayleigh–Taylor instability for different
density ratios of 30 ($\mathrm{At}=0.935$), 1000 ($\mathrm{At}=0.998$), and
3000 ($\mathrm{At}=0.999$) at
$t\sqrt{\mathrm{At}}=0,0.25,0.5,0.75,1,1.25,1.5,1.75$, and $2$, from left to
right, for the mesh resolution of $64\times 256$.

##  6  Implementation of two-phase solver for magnetic flows

In this section, the previously introduced two-phase solver is extended to
account for magnetic flows. To achieve this, the effect of the Lorentz force
is incorporated into the momentum equation, taking into consideration the
magnetic permeability jump across the interface. The governing equations of
two-phase flows along with the magnetostatic equation are solved, and the
detailed numerical discretization is presented.

###  6.1  Discretization of governing equations for incompressible flows under
magnetic fields

The Lorentz force quantifies the force experienced by fluids due to
electromagnetic fields, given as $\mathbf{J}\times\mathbf{B}$, where
$\mathbf{J}$ and $\mathbf{B}$ are electric current and magnetic flux
densities, respectively. This force can be written in the form of the Maxwell
stress tensor, $\tau^{\mathrm{M}}$, given as [25]

$\nabla\cdot\tau_{ij}^{\mathrm{M}}=\nabla\cdot\left(\frac{B_{i}B_{j}}{\mu_{\mathrm{m}}}-\frac{\mathbf{B}^{2}}{2\mu_{\mathrm{m}}}\delta_{ij}\right),$
(44)

where variable $\mu_{\mathrm{m}}$ denotes the magnetic permeability. This
force can be treated as a body force acting on a fluid, and, hence, the
updated momentum equation is expressed as

$\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)=-\frac{1}{\rho}\nabla
p+\frac{1}{\rho}\nabla\cdot\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}}\right)+\mathbf{g}+\frac{1}{\rho}\mathbf{F}_{\mathrm{sv}}+\frac{1}{\rho}\nabla\cdot\tau_{ij}^{\mathrm{M}}.$
(45)

When the field quantities do not change with time, Maxwell’s equations are
reduced to the electrostatic and magnetostatic case, which are given as

$\nabla\cdot\mathbf{B}=0,\ \ \ \ \ \ \ \ \ \ (\text{Gauss's law})$ (46a)
$\nabla\times\mathbf{E}=0,\ \ \ \ \ \ \ \ (\text{Faraday's law})$ (46b)
$\nabla\times\mathbf{H}=\mathbf{J},\ \ \ \ \ \ \ \ (\text{Maxwell--Amp\\`{e}re
law})$ (46c) and $\nabla\cdot\mathbf{J}=0,\ \ \ \ \ \ \ \ \ \ (\text{equation
of continuity})$ (46d)

where variables $\mathbf{H}$ and $\mathbf{E}$ denote the magnetic field and
electric field intensities, respectively. The magnetic flux density,
$\mathbf{B}$, is related to the magnetic field intensity, $\mathbf{H}$, using
the magnetic permeability, $\mu_{\mathrm{m}}=\mathbf{B}/\mathbf{H}$. In the
magnetostatic case, the behaviour of the magnetic field can be studied in the
absence of electric currents, since the electric charges are either at rest or
moving very slowly, so that the magnetic field induced by them can be
neglected. Consequently, there is no interaction between electric and magnetic
fields, and an electrostatic case or a magnetostatic case can be studied
separately.

Under the magnetostatic assumption, Eqs. (46a) and (46c) explain the evolution
of the magnetic field. One approach to solving the Maxwell–Ampère equation
while satisfying the magnetic field divergence-free constraint is the vector
potential formulation. In this method, the magnetic field, $\mathbf{B}$, is
defined as the curl of an auxiliary vector, $\mathbf{A}$, with the gauge
condition of $\nabla\cdot\mathbf{A}=0$, as
$\mathbf{B}=\nabla\times\mathbf{A}$. As a result, Gauss’s law of magnetism is
automatically satisfied, and Eq. (46c) is recast as

$\nabla\times\mathbf{H}=\nabla\times\left(\frac{1}{\mu_{\mathrm{m}}}\mathbf{B}\right)=\nabla\times\left(\frac{1}{\mu_{\mathrm{m}}}\nabla\times\mathbf{A}\right).$
(47)

For the two-dimensional case, $\mathbf{B}=\left(B_{x},B_{y},0\right)$, the
vector potential is reduced to $\mathbf{A}=\left(0,0,A_{z}\right)$, and Eq.
(47) will be simplified as

$\nabla\times\left(\frac{1}{\mu_{\mathrm{m}}}\mathbf{B}\right)=\nabla\times\left(\frac{1}{\mu_{\mathrm{m}}}\nabla\times\mathbf{A}\right)=\nabla\times\left(\frac{1}{\mu_{\mathrm{m}}}\frac{\partial
A_{z}}{\partial y}\mathbf{i}-\frac{1}{\mu_{\mathrm{m}}}\frac{\partial
A_{z}}{\partial
x}\mathbf{j}\right)=\nabla\left(-\frac{1}{\mu_{\mathrm{m}}}\nabla
A_{z}\right)=J_{z}.$ (48)

Equation (48) is a Poisson equation with variable coefficients and can be
discretized similarly to the pressure Poisson equation discussed earlier, Eq.
(5.2), provided that $\mu_{\mathrm{m}}$ and $A_{z}$ are defined at cell
centers (see Fig. 1).

Employing the obtained $A_{z}$ from Eq. (48), the components of the magnetic
field are given as

$B_{x}=\frac{\partial A_{z}}{\partial
y}=\frac{\delta_{\mathrm{2nd}}\,\,A_{z}}{\delta_{\mathrm{2nd}}\,\,y},$ (49a)
and $B_{y}=-\frac{\partial A_{z}}{\partial
x}=-\frac{\delta_{\mathrm{2nd}}\,\,A_{z}}{\delta_{\mathrm{2nd}}\,\,x},$ (49b)

at cell faces in $y-$ and $x-$directions, respectively (see Fig. 1). Finally,
for the two-dimensional case of $\mathbf{x}=\left(x,y\right)$, the components
of the Lorentz force are discretized as

$\left[\nabla\cdot\tau_{ij}^{\mathrm{M}}\right]_{x-\mathrm{comp}}=\frac{\delta_{2\mathrm{nd}}}{\delta_{2\mathrm{nd}}\,x}\left(\frac{\left(\overline{B_{x}}^{2\mathrm{nd}\,\,y}\right)^{2}-\left(\overline{B_{y}}^{2\mathrm{nd}\,\,x}\right)^{2}}{2\mu_{\mathrm{m}}}\right)+\frac{\delta_{2\mathrm{nd}}}{\delta_{2\mathrm{nd}}\,y}\left(\frac{\overline{B_{x}}^{2\mathrm{nd}\,\,x}\,\,\overline{B_{y}}^{2\mathrm{nd}\,\,y}}{\overline{\overline{\mu_{\mathrm{m}}}^{2\mathrm{nd}\,\,x}}^{2\mathrm{nd}\,\,y}}\right),$
(50a)
$\left[\nabla\cdot\tau_{ij}^{\mathrm{M}}\right]_{y-\mathrm{comp}}=\frac{\delta_{2\mathrm{nd}}}{\delta_{2\mathrm{nd}}\,x}\left(\frac{\overline{B_{x}}^{2\mathrm{nd}\,\,x}\,\,\overline{B_{y}}^{2\mathrm{nd}\,\,y}}{\overline{\overline{\mu_{\mathrm{m}}}^{2\mathrm{nd}\,\,x}}^{2\mathrm{nd}\,\,y}}\right)+\frac{\delta_{2\mathrm{nd}}}{\delta_{2\mathrm{nd}}\,y}\left(\frac{\left(\overline{B_{y}}^{2\mathrm{nd}\,\,x}\right)^{2}-\left(\overline{B_{x}}^{2\mathrm{nd}\,\,y}\right)^{2}}{2\mu_{\mathrm{m}}}\right).$
(50b)

The solution procedure introduced in the previous section for two-phase
nonmagnetic flows can be applied to the magnetic case as well with some
modifications. In step (2), the magnetic permeability field at $t^{n+1}$,
$\mu_{\mathrm{m}}^{n+1}$, should be computed according to the updated location
of the interface as well. The obtained magnetic permeability field will then
be employed to calculate the magnetic field at time step $n+1$, solving Eq.
(48), which in turn will be utilized to determine the Lorentz force in the
momentum equation. Figure 9 summarizes the complete procedure for implementing
the two-phase incompressible solver for magnetic flows.

Figure 9:  A flowchart representing the complete procedure for implementing
the two-phase incompressible numerical toolkit for magnetic flows.

###  6.2  Magnetic two-phase test cases

Three test cases are conducted in this section, namely, the deformation of
both a static and a sheared ferrofluid droplet as well as Rayleigh–Taylor
instability in magnetic fluids, to evaluate the performance and accuracy of
the implemented solver. The static droplet test case is designed to assess the
capability of the solver to accurately simulate the behaviour of the Lorentz
force at the interface for various magnetic field strengths. The numerical
results are validated by comparing them with experimental and analytical data.
In the second test, the deformation of a droplet in a shear flow is
investigated, considering both low and high capillary flow regimes under
varying magnetic field conditions. This test also involves comparing results
with theoretical solutions, particularly in the context of low magnetic field
values. Furthermore, within this test case, the impact of the magnetic
permeability ratio between the ferrofluid droplet and the surrounding flow on
its deformation and rotation is examined. The third benchmark is employed to
evaluate the solver’s performance in modelling the evolution of a complex
interface in the presence of different magnetic field densities and high
magnetic permeability jumps across the interface. Additionally, the impact of
the magnetic field on the growth rate of the Rayleigh–Taylor instability is
investigated and compared with the results obtained from linear analysis. It
is noteworthy to mention that since ferrofluids do not conduct electric
current, and in our test cases, no external current is imposed, the right-hand
side value of Eq. (48), $J_{z}$, is set to zero in the following numerical
simulations.

####  6.2.1 Deformation of a stationary magnetic droplet

In this test case, a liquid droplet of diameter $D=1$ is considered at the
center of a two-dimensional domain of $[0,4]\times[0,4]$ filled with gas, in a
stationary velocity field. In the case of zero gravity, similar to the static
droplet test case presented in Appendix E, the droplet remains at rest since
the pressure and capillary forces are balanced. However, if the gas and liquid
phases have different values of magnetic permeability, in the presence of a
magnetic field, the induced Lorentz force at the interface affects the
deformation of the droplet; the competition between the Lorentz force and the
surface tension force dictates the evolution of the droplet’s interface.
Suppose the density and viscosity are constant for both phases,
$\rho_{1}=\rho_{2}=1$ and $\mu_{1}=\mu_{2}=0.001$. A uniform magnetic field,
$\mathbf{B}_{0}=2~{}\mathrm{mT}$, is imposed from bottom to top and
$\mu_{\mathrm{m,l}}=9\,\mu_{\mathrm{m,g}}$. The capillary force attempts to
maintain the interface of the droplet in its initial shape. Nonetheless, since
the magnetic field lines will be distorted around the droplet’s interface, the
created Lorentz force causes the droplet to deform and stretch. In order to
better interpret the deformation of the droplet for different scenarios and
magnetic strengths, scale analysis is employed. To this end, the following
non-dimensional variables are introduced

$\mathbf{u^{*}}=\mathbf{u}/u_{0},\ \ \mathbf{B^{*}}=\mathbf{B}/B_{0},\
\,\mathbf{l}^{*}=\mathbf{l}/l_{0},\ \ t^{*}=t\,u_{0}/l_{0},\ \
p^{*}=p/\rho_{0}\,u_{0}^{2},$ (51)

where the zero subscripts refer to the initial value, and $l_{0}$ is the
length scale of the problem. Rewriting the momentum equation in the non-
dimensional form will then result in

$\displaystyle\frac{\partial\mathbf{u}^{*}}{\partial
t^{*}}+\frac{\partial}{\partial\mathbf{l}^{*}}\left(\mathbf{u}^{*}\mathbf{u}^{*}\right)=-\frac{\partial
p^{*}}{\partial\mathbf{l}^{*}}+\overbrace{\frac{\mu}{u_{0}\,l_{0}\,\rho_{0}}}^{\frac{1}{\mathrm{Re}}}\frac{\partial}{\partial\mathbf{l}^{*}}\left(\frac{\partial\mathbf{u}^{*}}{\partial\mathbf{l}^{*}}+\frac{\partial{\mathbf{u}^{*}}^{\mathrm{T}}}{\partial\mathbf{l}^{*}}\right)-\overbrace{\frac{\sigma}{\rho_{0}\,u_{0}^{2}\,l_{0}}\left(\nabla\cdot\mathbf{n}\right)\nabla\psi}^{\text{surface
tension}}$
$\displaystyle+\underbrace{\frac{{B_{0}^{2}}}{\rho_{0}\,\mu_{\mathrm{m}}\,{u_{0}}^{2}}\frac{\partial}{\partial\mathbf{l}^{*}}\left(\mathbf{B}^{*}\mathbf{B}^{*}-\frac{|\mathbf{B^{*}}|^{2}}{2}\delta_{i,j}\right)}_{\text{Lorentz
force}}.$ (52)

The ratio between the Lorentz force and the surface tension force can be
qualified by a non-dimensional number defined as

$\frac{\text{Lorentz force}}{\text{surface tension
force}}=\frac{\frac{B_{0}^{2}}{\rho_{0}\,\mu_{\mathrm{m}}\,{u_{0}}^{2}}}{\frac{\sigma}{\rho_{0}\,u_{0}^{2}\,l_{0}}}=\frac{B_{0}^{2}\,l_{0}}{\mu_{\mathrm{m}}\sigma}.$
(53)

Hence, for the case of $B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\ll 1$, the
surface tension force overcomes the Lorentz force and the droplet is expected
to retain its shape. Conversely, in the case of
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\gg 1$, since the Lorentz force is
greater than the capillary force, the droplet deforms and stretches. If
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma$ is of the order of one, the
Lorentz and surface tension forces are of the same magnitude. As a
consequence, oscillations in the deformation of the droplet could be observed,
as the Lorentz force stretches the interface and the surface tension force
tries to prevent the deformation.

This finding is then used to validate the behaviour of the implemented solver.
In the numerical setup, the thickness of the interface is set to
$\epsilon=\Delta x/2$, with the grid resolution of $101\times 101$, where a
no-slip boundary is imposed at all boundaries. Three cases of
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\ll 1$,
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\gg 1$, and
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\approx 1$ are investigated, where
the surface tension coefficient is set to 0.01, 1, and 100, respectively, and
the results are represented in Fig. 10. It can be appreciated in Fig. 10 that
for the case of $B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\gg 1$, the Lorentz
force is vertically stretching the droplet while for the
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\ll 1$ case, the droplet remains
unchanged in time. Additionally, when the Lorentz and capillary forces are of
the same order of magnitude, there is an oscillation in the deformation of the
droplet’s interface, evident in Fig. 10.

Figure 10: The evolution of the droplet under the constant magnetic field
imposed from bottom to top of the computational domain, with
$\mu_{\mathrm{m,l}}/\mu_{\mathrm{m,g}}=9$. Three cases of
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\gg 1$,
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\approx 1$, and
$B_{0}^{2}\,l_{0}/\mu_{\mathrm{m}}\,\sigma\ll 1$ are represented from left to
right, respectively.

The non-dimensional number introduced in Eq. (53) is the magnetic Bond number,
denoted as $\mathrm{Bo_{m}}$, mostly introduced as
$\mathrm{Bo}_{\mathrm{m}}=l_{0}\,\mu_{\mathrm{m,0}}\,H_{0}^{2}/2\sigma$ in the
literature. This parameter plays a critical role in various applications,
including the study of the dynamics and deformation of ferrofluid droplets in
the presence of a magnetic field. Ferrofluids are colloidal suspensions of
nanoscale magnetic particles, typically around 10 nm in size, dispersed in a
base fluid [36]. They were initially introduced by NASA in 1963, and since
then, ferrohydrodynamics has become a subject of significant interest in the
field of fluid mechanics [54]. Ferrofluids have found applications across
various fields, including microfluidics [55], biomedical applications such as
the treatment of retinal detachment and targeted drug delivery[56, 57],
droplet generation from nozzles [58], and heat transfer augmentation [59].
Understanding the behavior of ferrofluid droplets in the presence of magnetic
fields is essential for their practical applications. When subjected to a
uniform magnetic field, a ferrofluid droplet suspended in a viscous medium
elongates along the direction of the field, ultimately reaching a stable
equilibrium configuration. Researchers have conducted numerous studies to
investigate the deformation and dynamics of single ferrofluid droplets under
the influence of magnetic fields, employing analytical solutions, numerical
simulations, and experimental observations [60, 61, 34, 35].

To further validate our implemented two-phase solver for magnetic flows
quantitatively, we have leveraged the theoretical work presented by Afkhami et
al. [35], who explored the deformation of a ferrofluid droplet in a quiescent
fluid subjected to a uniform magnetic field. In their study, they established
a relation between the deformation of a ferrofluid and the magnetic Bond
number. This theoretical solution is derived under the assumption that the
droplet maintains its ellipsoidal shape during elongation due to the presence
of the magnetic field [35]. Consequently, the extent of deformation of the
droplet is quantified by introducing the aspect ratio denoted as $b/a$, where
$2a$ and $2b$ represent the major and minor axes of the droplet, respectively,
after it has undergone deformation and reached a steady state. Afkhami et al.
[35] modeled the magnetization of the ferrofluid droplet, denoted as
$\mathbf{M}_{\mathrm{ferrofluid}}$, as a linear function of the applied
magnetic field, given as $\mathbf{M}_{\mathrm{ferrofluid}}=\chi\mathbf{H}$,
where $\mathbf{H}$ represents the external uniform magnetic field strength,
and $\chi$ is the magnetic susceptibility of the ferrofluid droplet [36, 35].
Magnetic susceptibility, $\chi$, is a material property that quantifies the
magnetization response of the ferrofluid droplet to an applied magnetic field.
It is expressed as
$\chi=\left(\mu_{\mathrm{m,droplet}}/\mu_{\mathrm{m},0}-1\right)$ and is
assumed to remain constant in each phase in the analytical solution of Afkhami
et al. [35]. Consequently, the magnetic induction field is calculated as
$\mathbf{B}=\mu_{\mathrm{m},0}\left(\mathbf{H}+\mathbf{M}_{\mathrm{ferrofluid}}\right)=\mu_{\mathrm{m},0}\left(1+\chi\right)\mathbf{H}$
[35]. The theoretical finding regarding the droplet deformation in an external
magnetic field reported by Afkhami et al. [35] is presented as

$\mathrm{Bo}_{\mathrm{m}}=\left(\frac{1}{\chi}+k\right)^{2}\left(\frac{b}{a}\right)^{\frac{1}{3}}\left[2\frac{b}{a}-\left(\frac{b}{a}\right)^{-2}-1\right],$
(54)

where $k$ is the demagnetizing factor that is calculated as

$k=\left(\frac{1-E^{2}}{2E^{3}}\right)\left(\ln{\frac{1+E}{1-E}}-2E\right),\ \
\ \ \mathrm{with}\ \ \ \ E=\sqrt{1-a^{2}/b^{2}}.$ (55)

Here, a test case similar to the previous test is conducted to compare our
results with the theoretical solution of Eq. (54) and other numerical and
experimental results existing in the literature. In this test case, a circular
ferrofluid droplet with a radius of $R_{0}=0.5$ is placed at the center of a
computational domain of $[0,2]\times[0,6]$ filled with gas. The mesh
resolution of $100\times 300$ and time step $\Delta t=0.001$ were used for all
simulations. The density and viscosity are set to be the same for both phases,
$\rho_{1}=\rho_{2}=1$ and $\mu_{1}=\mu_{2}=0.01$, and a constant surface
tension value of $\sigma=1$ was employed. The magnetic permeability of the gas
is set to $\mu_{\mathrm{m,g}}=\mu_{\mathrm{m},0}=1$, therefore, to change the
magnetic permeability ratio between two phases, the magnetic permeability of
the ferrofluid droplet is altered. Figure 11(a)[left] illustrates the initial
magnetic field configuration for the case where $\chi=2$ and
$\mathrm{Bo}_{\mathrm{m}}=3$, along with the initial shape of the ferrofluid
droplet at $t=0$, as an example. As can be observed in this figure, the
magnetic field lines near the droplet interface are distorted due to the
varying magnetic permeability of the two phases. These distorted magnetic
field lines induce magnetic forces at the droplet interface, resulting in the
deformation of the ferrofluid droplet. The deformed droplet for this specific
case at steady state along with the corresponding magnetic field are
represented in Fig. 11(a)[middle]. As expected, the droplet has elongated
along the magnetic field lines due to the presence of magnetic forces. Figure
11(a)[right] demonstrates the forces acting on the ferrofluid droplet
interface, namely, the magnetic force (depicted in black) and the capillary
force (shown in red). It is visually evident from this figure that the
magnetic force exhibits a higher amplitude, enabling it to overcome the
capillary force and induce elongation. The surface tension force, primarily
concentrated at the poles of the droplet where high curvature is present,
opposes the magnetic forces, attempting to preserve the initial shape of the
droplet.

Figure 11(b) displays the droplet deformation for eight different magnetic
Bond numbers of $\mathrm{Bo}_{\mathrm{m}}=0.25$, 0.5, 1, 2, 3, 6, 8, and 10,
with the susceptibility value of $\chi=2$ for all cases. It can be seen that
for higher magnetic Bond numbers, the droplet deformation is more pronounced,
owing to the increased magnetic forces, since the surface tension coefficient
remains consistent across all cases. As the magnetic Bond number increases,
indicating stronger magnetic forces, the magnetic force becomes more effective
in overcoming capillary forces, further deforming the droplet.

In Fig. 11(c), the deformation of the ferrofluid droplet is explored under
three different magnetic susceptibility values of $\chi=2$, 5, and 20, across
various magnetic Bond numbers. The results are compared with analytical
solutions, demonstrating a close agreement between numerical and theoretical
predictions. As expected, an increase in the magnetic Bond number leads to
more significant changes in the aspect ratio. The maximum error between the
numerical and analytical results in the studied cases is approximately 4.2%.
This discrepancy can be attributed to the theoretical assumption that the
droplet is axisymmetric and retains an elliptical shape during deformation
[35]. However, this assumption may not hold true, particularly for cases with
higher susceptibility and magnetic Bond numbers [35], where more significant
discrepancies with the analytical solution are evident. Additionally, the
obtained results for $\chi=2$, 5, and 20 cases are compared with numerical
results from Afkhami et al. [35], showing good agreement, with some
discrepancy due to the differences in the nature of simulations (axisymmetric
in [35] vs. two-dimensional in this study).

Furthermore, in the case where $\chi=20$, we compare our results with the
experimental data obtained by Bacri and Salin [62]. While a reasonable
agreement was observed, discrepancies, particularly for higher magnetic Bond
numbers, can be attributed to the three-dimensional nature of the experiments
and the assumption of constant surface tension in our simulations. In reality,
it is observed that the interfacial tension of the ferrofluid droplet depends
on the magnetic field, especially for high magnetic field strengths [35].

(a)

(b)

(c)

Figure 11: (a) Configuration of the ferrofluid droplet along with magnetic
field at [left] $t=0$, [middle] steady state, and [right] depiction of the
magnetic force (in black) and capillary force (in red) acting at the interface
for the case $\chi=2$ and $\mathrm{Bo}_{\mathrm{m}}=3$. (b) droplet
deformation across eight different magnetic Bond numbers with a susceptibility
value of $\chi=2$. (c) Comparison of analytical, numerical, and experimental
results for three different magnetic susceptibility values of $\chi=2$, 5, and
20, under varying magnetic Bond numbers.

####  6.2.2 Deformation of a sheared magnetic droplet

As previously mentioned, ferrofluid droplets have garnered significant
attention from researchers due to their diverse applications in areas such as
biomedicine, microfluidics, and rheology [63]. In some applications, these
droplets are not only subjected to an external magnetic field, as discussed in
the previous test case, but also experience the presence of hydrodynamic flow.
This flow can influence the droplet’s deformation, inclination, and potential
breakup into smaller droplets. Consequently, studying the deformation of
sheared droplets becomes crucial across various industrial applications that
utilize emulsions. The deformation of the sheared droplet depends on various
parameters, such as surface tension, shear rate, magnetic field strength, and
the viscosity ratio between the droplet and the suspending fluid. In addition
to the magnetic Bond number, $\mathrm{Bo}_{\mathrm{m}}$, another non-
dimensional parameter known as the capillary number, $\mathrm{Ca}$, which
quantifies the ratio between shear and surface tension forces, also influences
the droplet deformation. Hassan et al. [64] conducted a numerical study on the
deformation of sheared droplets under a uniform external magnetic field using
the finite element method. Their study investigated the effects of shear rate,
magnetic field strength, and magnetic field direction on ferrofluid
deformation. They found that in the low capillary regime, magnetic forces
dominate over shear forces, exerting primary control over droplet dynamics
[64]. Additionally, increasing the magnetic field strength, i.e., increasing
the magnetic Bond number, enhances droplet deformation [64]. In a subsequent
study by Hassan and Wang [65], the impact of viscosity ratio and magnetic
field direction on droplet deformation and breakup was explored under the
constraint of a low Reynolds number ($\mathrm{Re}\leq 0.03$). Their results
indicated that when the magnetic field is applied at an angle of $45^{\circ}$,
droplet elongation is more pronounced, leading to accelerated droplet breakup.
Conversely, varying the magnetic field angle to $0^{\circ}$ or $135^{\circ}$
suppresses droplet breakup. To our knowledge, no previous study has
investigated sheared droplet deformation by varying the susceptibility value
of the ferrofluid droplet. Consequently, in this test case, we investigate
sheared droplet deformation for low capillary number to validate the
implemented numerical solver. Additionally, we explore the effect of magnetic
susceptibility ratio between the ferrofluid droplet and its surrounding medium
on the dynamics and deformation of the ferrofluid droplet in both low and high
capillary regimes.

Figure 12 illustrates the schematic of this test case, including a ferrofluid
droplet with radius $R_{0}=0.5$ suspended in another nonmagnetic viscous
fluid. The droplet is positioned at the center of a square computational
domain with dimension $W_{\mathrm{domain}}=6$, and a velocity profile of
$u=\dot{\gamma}y$ is imposed, where $\dot{\gamma}$ represents the
corresponding shear rate. A uniform magnetic field of $\mathbf{H}_{0}$ at the
angle of $90^{\circ}$ is applied. It is worth noting that numerical studies
model simple shear flows by moving two confining walls, which can introduce
confinement effects on the deformation of the droplet. Studies by Kennedy [66]
and Guido [67] on the droplet deformation in a simple shear flow suggest that
the confinement effect on droplet deformation is negligible when
$2R_{0}/W_{\mathrm{domain}}<0.4$. In our simulations,
$2R_{0}/W_{\mathrm{domain}}$ is equal to 0.16, effectively eliminating the
confinement effect on droplet deformation. In the presented simulations, the
density and viscosity of both phases are assumed to be equal, $\rho=1$ and
$\mu=0.1$. The top and bottom boundaries of the computational domain are
treated as moving walls with velocities of
$+\dot{\gamma}W_{\mathrm{domain}}/2$ and $-\dot{\gamma}W_{\mathrm{domain}}/2$,
respectively, producing constant shear rate of $\dot{\gamma}$, while the left
and right boundaries are periodic. By setting surface tension coefficient and
shear rate to 1 and 0.4, respectively, the capillary number is adjusted to
$\mathrm{Ca}=\mu_{\mathrm{m},0}\,R_{0}\,\dot{\gamma}/\sigma=0.02$, and the
deformation of the droplet is investigated for different values of
$\mathrm{Bo_{m}}$ and $\chi$. To quantify droplet deformation, Taylor’s
deformation parameter [68, 69] is calculated, given as

$D=\frac{L-B}{L+B},$ (56)

where $L$ and $B$ represent the major and minor axes of the deformed droplet
(see Fig. 12). The pioneering work of Taylor indicating the deformation
parameter of a buoyant droplet suspended in another viscous fluid under a
shear flow in the Stokes flow limit [68, 69] was extended by Jesus et al. [70]
to a sheared ferrofluid droplet under an external magnetic field in the limit
of both small $\mathrm{Ca}$ and $\mathrm{Bo}_{\mathrm{m}}$ numbers. According
to the asymptotic theory derived by Jesus et al. [70], the deformation of a
sheared ferrofluid droplet in the presence of a uniform magnetic field is
given as

$D=\frac{\sqrt{\left[\alpha\left(\nu_{\mathrm{r}}\right)\mathrm{Ca}\right]^{2}+\left[\beta\left(\chi_{\mathrm{droplet}}\right)\mathrm{Bo}_{\mathrm{m}}\right]^{2}}}{2+\frac{1}{3}\beta\left(\chi_{\mathrm{droplet}}\right)\mathrm{Bo}_{\mathrm{m}}},$
(57)

where $\alpha\left(\nu_{\mathrm{r}}\right)$ and
$\beta\left(\chi_{\mathrm{droplet}}\right)$ are the functions of the viscosity
ratio between the droplet and the surrounding fluid, $\nu_{\mathrm{r}}$, and
the susceptibility of the ferrofluid droplet, $\chi_{\mathrm{droplet}}$,
respectively. These functions are calculated as

$\alpha\left(\nu_{\mathrm{r}}\right)=\frac{19\nu_{\mathrm{r}}+16}{8(\nu_{\mathrm{r}}+1)}$
(58a) and
$\beta\left(\chi_{\mathrm{droplet}}\right)=\frac{3\,\chi_{\mathrm{droplet}}\left(2\,\chi_{\mathrm{droplet}}+1\right)}{4\left(\chi_{\mathrm{droplet}}+3\right)^{2}}.$
(58b)

Figure 12: Schematic of a ferrofluid droplet suspended in another nonmagnetic
viscous fluid in a simple shear flow under a uniform external magnetic field,
$\mathbf{H}_{0}$.

Figure 13[left] presents the droplet deformation parameter with
$\mathrm{Ca}=0.02$ for different magnetic bond numbers, ranging from 0 to 6,
for three different droplet susceptibility values of $\chi=1$, 2, and 3. Since
the simulations are performed in the low capillary regime, it is expected that
magnetic forces dominate over shear forces and as the magnetic Bond number
increases, the droplet deformation is anticipated to increase. Figure 13[left]
visually confirms that the numerical results capture this expected behaviour.
As shown in this figure, the results for $\mathrm{Bo_{m}}<1$ are in good
agreement with the theoretical solution. Please refer to the inset in this
figure, which depicts analytical and numerical solutions for the magnetic Bond
number ranging from 0 to 0.8. However, as the magnetic Bond number increases,
the numerical results deviate from the theoretical solution derived by Jesus
et al. [70]. This behavior is expected since the proposed theoretical solution
is only valid for small magnetic Bond numbers ($\mathrm{Bo_{m}}\ll 1$).
Furthermore, it is evident that by increasing the susceptibility value, the
numerical results differ more pronouncedly from the small perturbation
analysis of the theoretical solution. According to the numerical results,
increasing the magnetic susceptibility value of the droplet significantly
affects its deformation. For example, for $\chi=3$ and
$\mathrm{Bo}_{\mathrm{m}}=6$, the deformation parameter is approximately 2
times higher than the same case with $\chi=1$. Figure 13[middle] depicts the
steady-state droplet at $\mathrm{Bo}_{\mathrm{m}}=1$ for the three different
magnetic susceptibility values of $\chi=1$, 2, and 3. As can be seen from this
figure, increasing the magnetic susceptibility value by three times results in
an elongated, thin droplet shape compared to the $\chi=1$ case. This behavior
can be explained by investigating the existing Maxwell stresses at the droplet
interface for different susceptibility values. In Fig. 14, the magnetic field
for three different values of $\chi=1$, 2, and 3 is shown at $t=0$ and the
steady state. According to this figure, a larger $\chi$, or in other words, a
larger magnetic permeability discontinuity between the ferrofluid droplet and
the surrounding fluid, will result in a greater distortion of magnetic field
lines around the droplet interface. This higher bending in the magnetic field
lines induces a larger Lorentz force at the interface, which in turn causes
greater droplet elongation, as the surface tension (restoring force) is
constant for all cases. This conclusion emphasizes that changing the
susceptibility of the ferrofluid droplet also plays an important role in its
dynamics and deformation. In the study by Hu et al. [71], it was also shown
that the magnetic field lines inside the ferrofluid droplet at $t=0$ are
aligned with the externally imposed magnetic field. However, the magnetic
field will be distorted near the interface of the droplet. Hu et al. [71]
found the analytical solution for the magnetic field intensity inside and
outside of the circular cylinder within an externally vertical uniform
magnetic field. According to their analytical solution, the magnetic field
intensity inside the cylinder will be given as $\mathbf{H}=A\mathbf{e}_{y}$,
where
$A=H_{0}\,(2\mu_{\mathrm{m},2})/(\mu_{\mathrm{m},1}+\mu_{\mathrm{m},2})$.
Variables $\mu_{\mathrm{m},1}$ and $\mu_{\mathrm{m},2}$ denote the magnetic
permeability of the surrounding gas and the ferrofluid droplet, respectively.
We measured the numerical magnetic field intensity for the three different
cases of $\chi=1$, $2$, and $3$, shown in Fig. 14, and compared it to the
analytical solution. According to our results, the magnetic field intensity
for the cases $\chi=1$, $2$, and $3$ are $1.32H_{0}$, $1.51H_{0}$, and
$1.62H_{0}$, respectively. These results are in good agreement with the
numerical findings of $1.33H_{0}$, $1.5H_{0}$, and $1.6H_{0}$ for the cases
$\chi=1$, $2$, and $3$, respectively.

Figure 13: [left] Comparison of analytical and numerical results for three
different magnetic susceptibility values of $\chi=1$, 2, and 3 under varying
magnetic Bond numbers and fixed capillary number, $\mathrm{Ca=0.02}$. [middle]
Ferrofluid droplet at steady state for three cases of $\chi=1$, 2, and 3, with
capillary number and magnetic Bond number of 0.02 and 1, respectively. The
dashed line connects the poles of the droplet for each case. [right]
Ferrofluid droplet inclination, $\theta$, for three different magnetic
susceptibility values of $\chi=1$, 2, and 3 under varying magnetic Bond
numbers and fixed capillary number, $\mathrm{Ca=0.02}$.

In Fig. 14, the acting magnetic and surface tension forces at the interface
are shown in black and red, respectively, at the steady-state situation. Since
we are in the low capillary regime, shear forces play a minor role in droplet
deformation and are neglected in the figure. For $\chi=1$, since the droplet
is still in a circular shape, the surface tension is present around the entire
interface. However, for $\chi=2$ and 3, as the droplet has elongated, the
surface tension force is mainly focused at the poles, competing with the
magnetic force. Increasing the magnetic Bond number requires a longer
simulation time for the droplet reaches a steady state, as previously reported
by Jesus et al. [70].

Figure 14: [top]Configuration of the ferrofluid droplet along with magnetic
field at $t=0$, [middle] steady state, and [bottom] depiction of the magnetic
force (in black) and capillary force (in red) acting at the droplet interface
for three cases of $\chi=1$, 2, and 3, with $\mathrm{Bo}_{\mathrm{m}}=2$ and
$\mathrm{Ca}=0.02$.

Another important parameter to consider when investigating ferrofluid droplet
deformation is the rotation of the droplet under various conditions,
quantified by the angle measured counterclockwise from the positive
$x-$direction to the major axis of the droplet, $\theta$. Previous studies
have explored the rotational angle of the droplet under different magnetic
Bond numbers in both low and high capillary regimes. For instance, Hassan et
al. [64] demonstrated that for $\mathrm{Ca}\approx 0.2$ and an imposed
perpendicular magnetic field, increasing the magnetic Bond number leads to an
increase in $\theta$ angle. This effect is attributed to the combined
influence of the magnetic field and shear flow. However, the impact of
magnetic susceptibility on droplet rotation has not yet been discussed in
detail. In Fig. 13[right], the rotation of the droplet for Ca=0.02 under
different magnetic Bond numbers is presented. For a susceptibility value of 1,
similar to the behavior reported in previous studies [64], an increase in the
magnetic Bond number leads to an increase in the angle $\theta$ until it
reaches a constant value. The shear flow attempts to rotate the droplet at an
angle of 45∘, while the external magnetic field aims to elongate the droplet
along its vertical field lines. The superposition of these two effects results
in the rotation of the droplet. In the low capillary regime, as the magnetic
Bond number increases, magnetic forces overcome the shear effect, causing the
$\theta$ angle to increase along with the droplet deformation parameter.
According to the numerical results, for a constant magnetic Bond number less
than 2, the higher the droplet’s susceptibility, the greater the rotational
angle. This can be attributed to the increased Maxwell stress for higher
susceptibility values, which further prevents the droplet from deflecting at
the 45∘ angle. However, for magnetic Bond numbers of 3 and greater, the
rotational angle becomes almost the same for all different susceptibility
values, with no noticeable difference observed between them. These results
suggest that not only does increasing the susceptibility value significantly
affect droplet deformation, but also for lower magnetic Bond numbers, higher
susceptibility leads to a greater rotational angle, causing the droplet to
become more vertical. This is an important observation that can be further
utilized when studying sheared droplet deformation in various applications.

(a)  (b)

(c)

Figure 15: (a) Comparison of theoretical and numerical results for two
different magnetic susceptibility values, $\chi=1$ and 2, under varying
magnetic Bond numbers, while maintaining a fixed capillary number,
$\mathrm{Ca=0.2}$. The inclination of the ferrofluid droplet, $\theta$, for
different cases is also included in this figure. (b) Visualization of a
ferrofluid droplet at steady state for two cases with $\chi=1$ and 2, each at
a capillary number of 0.2, across six different magnetic Bond numbers (0, 0.5,
1, 2, 3, and 4). The red contour represents the initial droplet configuration
at $t=0$. Notably, for the case with $\chi=2$ and magnetic Bond numbers of 3
and 4, the ferrofluid droplet has not reached a steady state, leading to a
loss of its ellipsoidal shape as it continues to stretch and thin in the
middle. (c) The sheared droplet deformation and breakup for the case $\chi=2$
with a magnetic Bond number of 4 at $t=0$, 0.5, 1, 1.5, 2, 2.2, 2.4, and 2.5,
presented from left to right, respectively.

The effect of the value of susceptibility on the deformation of a ferrofluid
droplet is also investigated in the high capillary regime, Ca=0.2, achieved by
increasing the shear rate. In Fig. 15(a), the deformation parameters of
deformed droplets with two different susceptibilities, $\chi=1$ and 2, are
represented for magnetic Bond numbers ranging from 0 to 4, while keeping
Ca=0.2 constant. The theoretical results are also included in this figure. As
expected, since we are in the higher capillary regime, the theoretical
approach does not accurately predict the $D$ values. The theoretical analysis
is valid only for $\mathrm{Ca}\ll 1$, and its predictions become less reliable
as the magnetic Bond number increases. According to this figure, for a
constant magnetic Bond number, increasing the $\chi$ values results in an
increase in the deformation parameter, similar to the observations in the
previous section. However, it is notable that the obtained deformation
parameters in this case are greater than those in the previous case with
Ca=0.02. This is because the greater shear rate allows the external shear flow
to more effectively influence the droplet’s deformation. Moreover, it is
evident that by increasing the magnetic Bond number, the angle $\theta$
increases and approaches an almost constant value. However, in this case, due
to the stronger shear flow, the droplet tends to deflect more toward the 45∘
shear angle (please compare Fig.13[right] and Fig. 15(a)). The calculated
deflection angles in the case of $\chi=2$ for magnetic Bond numbers below 2,
are higher than those for $\chi=1$. This is because for $\chi=2$, the magnetic
field lines are deflected further near the droplet’s interface due to a higher
magnetic permeability discontinuity across the interface, resulting in a
stronger force to elongate the droplet along the field lines and leading to
higher $\theta$ angles. It is worth noting that the results for $\chi=2$ under
magnetic Bond numbers of 3 and 4 are absent in this figure. This is because
for these magnetic Bond number values, the ferrofluid droplet has not reached
a steady-state condition, as shown in Fig. 15(b), and it continues deforming
until it undergoes breakup. Figure 15(c) illustrates the sheared droplet’s
deformation and breakup for the case where $\chi=2$ and $\mathrm{Bo_{m}}=4$ at
eight different time steps, up to the moment when breakup occurs. The magnetic
Bond number of approximately 3 is the critical value below which the droplet
does not rupture. In previous studies, other properties such as the effect of
viscosity ratio on droplet breakup were investigated [65]. In this test case,
we have demonstrated that increasing the magnetic susceptibility value of the
droplet can also cause droplet breakup, introducing the concept of a critical
susceptibility value above which, for a constant Ca and $\mathrm{Bo_{m}}$
values, breakup occurs. Therefore, the magnetic permeability ratio between the
droplet and its surrounding flow not only affects its deformation but can also
lead to a breakup mechanism.

####  6.2.3 Rayleigh–Taylor instability in magnetic fluids

The Rayleigh–Taylor instability, which occurs when the interface between
fluids of different densities is accelerated towards the direction of the
heavy fluid, is widely studied for conducting fluids in the presence of a
magnetic field as well. Magneto-Rayleigh–Taylor instability is of importance
in different applications, e.g., Z-pinch fusion reactors [72], magnetic flux
compression [73], and laser-driven inertial confinement fusion [74]. In this
test case, we investigate the Rayleigh–Taylor instability for magnetic fluids
in the presence of a tangential quasi-static magnetic field to validate the
performance of the implemented solver. By employing the linear analysis on the
Rayleigh–Taylor instability, Awasthi [75] derived the following quadratic
dispersion relation for two viscous, incompressible, and magnetic fluids in
the presence of a tangential magnetic field

$\overbrace{\left[\rho_{1}+\rho_{2}\right]}^{\alpha}\omega^{2}+\overbrace{\left[4k^{2}(\mu_{1}+\mu_{2})\right]}^{\beta}\omega+\overbrace{\left[\left(\rho_{1}-\rho_{2}\right)gk+\sigma
k^{3}+\frac{k^{2}H_{0}^{2}\left(\mu_{\mathrm{m,2}}-\mu_{\mathrm{m,1}}\right)^{2}}{\left(\mu_{\mathrm{m,1}}+\mu_{\mathrm{m,2}}\right)}\right]}^{\gamma}=0.$
(59)

In Eq. (59), $\omega$ and $k$ variables are the growth rate and wave number of
a small perturbation at the interface, respectively, and subscript 1 (2)
refers to the lighter (heavier) fluid. According to the instability criteria
of Routh–Hurwitz, the stability condition is given as $\alpha>0$, $\beta>0$,
and $\gamma>0$ [75]. Coefficients $\alpha$ and $\beta$ are always positive,
and, therefore, the stability condition is reduced to $\gamma>0$. The last
term of the coefficient $\gamma$ represents the effect of the imposed magnetic
field on the growth rate of the Rayleigh–Taylor instability. This term is
always positive since the value of all the components in this term is
positive. Hence, the tangential magnetic field always has a stabilizing effect
on the Rayleigh–Taylor instability growth, similar to the effect of surface
tension. For the sake of a more consistent analysis, by using the following
dimensionless variables as

$\hat{k}=kL,\ \ \hat{\rho}=\frac{\rho_{1}}{\rho_{2}},\ \
\hat{\mu}=\frac{\mu_{1}}{\mu_{2}},\ \
\hat{\mu}_{\mathrm{m}}=\frac{\mu_{\mathrm{m,1}}}{\mu_{\mathrm{m,2}}},\ \
\hat{\sigma}=\frac{\sigma}{\rho_{2}\,g\,L^{2}},\ \
\hat{H}=\frac{H_{0}\sqrt{\hat{\rho}}}{\sqrt{\left(1-\hat{\rho}\right)\,g\,L}}\sqrt{\frac{\mu_{\mathrm{m,2}}}{\rho_{2}}},$
(60)

where $L$ is the characteristic length scale of the domain, the non-
dimensional stability condition is expressed as

$1-\frac{\hat{\sigma}\hat{k}^{2}}{\left(1-\hat{\rho}\right)}-\frac{\hat{k}\hat{H}^{2}\left(1-\hat{\mu}_{\mathrm{m}}\right)^{2}}{\hat{\rho}\left(1+\hat{\mu}_{\mathrm{m}}\right)}\leq
0.$ (61)

For the state of the marginal stability, when Eq. (61) is zero, the critical
value of the transverse magnetic field, $\hat{H}_{\mathrm{c}}$, is given as

${\hat{H}_{\mathrm{c}}}^{2}=\frac{\left[\left(1-\hat{\rho}\right)-\hat{\sigma}\hat{k}^{2}\right]\hat{\rho}\left(1+\hat{\mu}_{\mathrm{m}}\right)}{\hat{k}\left(1-\hat{\rho}\right)\left(1-\hat{\mu}_{\mathrm{m}}\right)^{2}}.$
(62)

Hence, depending on whether the imposed transverse magnetic field is larger or
smaller than the critical magnetic field, the system becomes stable or
unstable, respectively.

To investigate the performance of the implemented solver for two-phase
magnetic flows, the evolution of the Rayleigh–Taylor instability under
different magnetic field intensities is simulated. To this end, the same
initial condition as those presented in Sec. 5.4 for the Rayleigh–Taylor
instability test case, with the grid resolution of $32\times 128$,
$\hat{\rho}=1/3\,(\rho_{1}=1,\rho_{2}=3)$, and
$\hat{\mu}_{\mathrm{m}}=0.01\,(\mu_{\mathrm{m,1}}=10^{-3},\mu_{\mathrm{m,2}}=10^{-1})$
is considered. Three cases of
$\hat{H}_{\mathrm{0}}<\hat{H}_{\mathrm{c}}\,({H}_{\mathrm{0}}=2)$,
$\hat{H}_{\mathrm{0}}\approx\hat{H}_{\mathrm{c}}\,({H}_{\mathrm{0}}=5)$, and
$\hat{H}_{\mathrm{0}}>\hat{H}_{\mathrm{c}}\,({H}_{\mathrm{0}}=10)$ are
examined 111The critical value of the transverse magnetic field,
${H}_{\mathrm{c}}$, is expected to be around 4.. The purpose is to
qualitatively compare the numerical results with the ones obtained from linear
perturbation analysis. Based on Fig. 16(a), it is visually evident that the
presence of a transverse magnetic field, regardless of its magnitude, tends to
suppress the growth of the Rayleigh–Taylor instability. This observation is
consistent with the results obtained from other studies [76, 75], stating that
the Lorentz force at the interface acts as a restoring force. However,