New Immersed Boundary Method with Irrotational Discrete Delta Vector for Droplet Simulations with Large Density ratio
Chia Rui Ong Hiroaki Miura







New Immersed Boundary Method with Irrotational Discrete Delta Vector for Droplet Simulations with Large Density ratio
    Chia Rui Ong and Hiroaki Miura
     December 30, 2023
=====================================================================================================================






The Immersed Boundary Method (IBM) is one of the popular one-fluid mixed Eulerian-Lagrangian 
methods to simulate motion of droplets. While the treatment of a moving complex 
boundary is an extremely time consuming and formidable task in a 
traditional boundary-fitted fluid solver, the one-fluid methods provide a relatively easier way to 
track moving interfaces on a fixed Cartesian grid since the regeneration of a 
mesh system that conforms to the interface at every time step can be avoided. In the IBM, a series 
of connected Lagrangian markers are used to represent a fluid-fluid interface and the boundary 
condition is enforced by adding a forcing term to the Navier-Stokes equations. To enable this, a 
discrete delta function is necessary for interpolation of velocity field and extrapolation of immersed 
boundary forcing between Eulerian grids and Lagrangian markers.  It is known that the IBM suffers two problems. One is spontaneous generation of unphysical kinetic 
energy, which is known as parasitic currents and also appear in other one-fluid methods as well, and 
the other is spurious reconstruction of interface. These two problems need to be solved for useful 
long-time-scale simulations of droplets with high density ratio and large surface tension. This 
work detects that the discrete delta function is the cause of unphysical parasitic currents. 
Specifically, the irrotational condition is not preserved when the common discrete delta function is 
used to spread the surface tension from Lagrangian markers to Cartesian grid cells. To solve this problem, 
a new scheme that preserves the irrotational condition is proposed to remove the spurious 
currents. Furthermore, for a smooth reconstruction of an interface, a 
B-spline fitting by least squares is adopted to relocate the Lagrangian markers. The 
conventional and new interpolation schemes are implemented in a multigrid finite volume Direct 
Numerical Simulation (DNS) solver and are subjected to standard test cases. It is confirmed that the 
unphysical parasitic currents are substantially reduced and droplet's surface fluctuation is eliminated in the 
new scheme. The new scheme is also applied for simulations of axisymmetric free-fall droplet. In this test, the numerical results 
coincide well with experimental data. 





















			





§ INTRODUCTION


The Direct Numerical Simulation (DNS) of multiphase flow is a challenging problem. The difficulties 
lie in the non-linearity of the interfacial boundary conditions and the correct representation of 
the interface that is freely transported by fluids. 
There are generally two types of methods to tackle these difficulties according to how the location 
of an interface is tracked down: the interface-fitted method and the one-fluid Eulerian-Lagrangian 
models (in short form as one-fluid model). In the interface-fitted model 
<cit.><cit.><cit.><cit.>, a grid is regenerated at 
every time step so that the interfacial grid cells constantly conform to the interfacial boundary. 
The grid regeneration process is highly computational demanding. Furthermore, it has difficulties in 
dealing with the separation and merging of two interfaces. Despite these disadvantages, it has advantages that the interface 
can be accurately represented and discretization of the interfacial boundary condition is straigtforward. 
On the other hand, the one-fluid models have been becoming popular in recent decades because the Navier-Stokes equations are 
solved on the Cartesian grid. The interface is tracked down either by some explicit functions, which 
are transported by the advection equation, or Lagrangian markers. This interface tracking allows a substantial reduction of computational costs. The discretization of equations is much more 
straightforward than the curvilinear mesh used in the interface-fitted model. Although the time 
consuming process of grid regeneration at every time step can be avoided, the adaptive mesh is 
sometimes needed for better accuracy. 

One-fluid models can be further categorized into several methods according to the way to ensure 
that the interfacial boundary conditions are satisfied and interfaces are tracked. There are two schemes 
to treat the interfacial boundary conditions: 1) the sharp interface method 
and the 2) continuous interface method. In the sharp interface method, the interfacial boundary 
conditions are explicitly incorporated into the difference equations when discretizing the 
Navier-Stokes equations. Examples of the sharp interface method are the Immersed Interface 
Method (IIM) <cit.> and the ghost-fluid method <cit.>. In the continuous interface 
method proposed by Brackbill et al <cit.>, a smooth forcing term which has a support of 
several grid spacings is added into the Navier-Stokes equations and the force is spread to the 
Cartesian grid by an extrapolation scheme. A significant advantage of the continuous interface 
method against the sharp interface method is that it is straightforward to implement it in an 
existing fluid solver without changing discretization schemes. Meanwhile, a disadvantage is that a 
force term needs to be extrapolated to the Cartesian cell centers or faces by a regularized discrete 
delta function <cit.><cit.><cit.>. This reduces the accuracy of the model down to the first order due to the smearing of the interface <cit.><cit.><cit.>. 

Depending on the way of tracking interfaces, one-fluid models can be categorized into three main 
groups: (1) The Immersed Boundary Method (IBM) <cit.><cit.><cit.><cit.><cit.> in which 
interfaces are represented by a series of connected or disconnected Lagrangian markers; (2) The 
Volume-Of-Fluid (VOF) <cit.><cit.><cit.> in which a volume function is used to represent the volume fraction of a specific fluid in a grid cell; (3) The Level-Set (LS) <cit.><cit.><cit.> in which a smooth 
level set function is used to represent the location of an interface. There are other one-fluid 
models such as the diffuse interface method <cit.> in which a vastly different strategy is 
adopted. In the diffuse interface method, the density of fluid is determined by the generalized 
diffusion equation. 

In this work, we focus on the implementation of the IBM with the continuous interface method for 
simulations of a droplet. We believe that there are three important conditions that need to be 
satisfied for a successful long-time-scale simulation of a droplet or droplets with moderate-to-large 
density ratio between two fluids. First, the force spreading operator should be irrotational to avoid 
unphysical spurious flow parallel to the interface. Second, the sum of surface tension force over the 
computational domain should be as small as possible (ideally zero). Third, a stable reconstruction of the interface for a smooth 
distribution of curvature. However, we found that the IBM often fails to satisfy one or more 
conditions above. It will be shown in this work that failing to satisfy the first condition leads to 
the generation of the so-called parasitic currents. In other words, to eliminate the parasitic 
currents, the curl of discrete delta vector interpolated onto the Cartesian grid from the Lagrangian 
markers has to be zero strictly. Previous works have already discussed this problem 
in different one-fluid models (see <cit.> and <cit.>), and some 
solutions have been proposed. For example, Francois et al <cit.> have devised a 
balanced-force algorithm, which produces an exact balance between the pressure gradient force and the 
surface tension force, to eliminate it in the context of VOF. Jamet et al <cit.> used the diffuse interface method with 
an energy conserving discretization to remove it. In this work, we proposed a new approach that uses 
an irrotational extrapolation scheme based on the discrete gradient of the Heaviside function in the IBM formulation. 

In the IBM, the curvature of a surface at each marker can be approximated by the absolute value of difference of tangent vectors of adjacent vectors over the distance between them. As the distance between two adjacent markers is variable, the total surface tension force may not be zero. Moreover, the markers often move in a haphazard way and the curvature is very sensitive to the locations of markers due to the property of interpolation spline curve. This leads to unphysical flow fields near the boundary when the curvature fluctuates heavily. As a result, the latter two conditions are often not satisfied and stability is severely impeded. In order to solve these problems, we decided to adopt the global B-spline curve to interpolate and redistribute the markers evenly at every time step.  

The new model proposed above satisfies the first and third conditions while the condition two is only approximately satisfied, because the new spreading method may not conserve the total surface tension force. In the droplet simulations, this is usually a minor issue compared to the instability caused by the parasitic currents and curvature fluctuation. To demonstrate the validity of our new approach, test cases have been done. The test cases are static and free oscillating droplet in both two dimensional and cylindrical axisymmetric domains. A further set of free-fall droplet simulations were performed to demonstrate its capability. 


The overview of this paper is as follows. Section 2 is devoted to the formulation and discretization 
schemes of the finite volume solver and the immersed boundary method. Details of the problem of parasitic currents and implementation of our new scheme are given in section 3. The numerical results are presented in section 4. The conclusion is given in the final section. 





§ CONTINUUM FORMULATION AND NUMERICAL DISCRETIZATION



 §.§ Flow solver

The formulation of the Navier-Stokes equations for incompressible flow with the immersed boundary 
(IB) forcing term is: 

    ρu⃗t + ρ(u⃗u⃗)  = -▽p + 
        (μ∇u⃗)+ f⃗_⃗i⃗b⃗,   and


    u⃗ = 0.

Equations (<ref>) and (<ref>) are the momentum equation and continuity 
equation respectively where u⃗ is velocity, ρ is density, p is pressure, and μ 
is kinematic viscosity. The forcing term f⃗_⃗i⃗b⃗ in equation (<ref>) is added to 
satisfy the interfacial boundary conditions and its specific form will be defined in the next 
section. In the IBM, the simple staggering Cartesian mesh system is used in the finite-volume fluid 
solver, while the fluid-fluid interface is represented by a series of connected Lagrangian markers 
(Figure <ref>). We employ the fractional step method <cit.> to solve equations above. In the fractional step method, the Navier-Stokes equations are solved in two steps by splitting the momentum equation into two parts: 

    ρ^n+1u⃗^* - u⃗^⃗n⃗/Δ t  = RHS + f⃗_⃗i⃗b⃗,   and


    ρ^n+1u⃗^n+1 - u⃗^*/Δ t  = -▽p^n+1.

where u^* is the intermediate velocity, which does not necessarily satisfy the divergence-free 
condition, Δ t is the time step, RHS is the advective and diffusive terms, and the 
superscript n means scalar or vector values at the n^th time step. The density is determined from 
the indicator function which will be given in the next section. The temporal discretization methods 
for the advection term and the viscous stress term are the Adams-Bashforth scheme and the 
Crank-Nicholson scheme, respectively. The spatial discretization of these two terms about a velocity 
component u for the control volume centered at face (i-1/2,j) as shown in Figure <ref> 
is given by (superscript is ignored):

    ∮ρ(u⃗ u)  dV    = ∑_all facesρ u u⃗· ds⃗
       = ρ_i,j ( u_i+1/2,j + u_i-1/2,j/2 )^2 Δ y_i-1/2,j - 
         ρ_i-1,j ( u_i-1/2,j + u_i-3/2,j/2 )^2 Δ y_i-1/2,j
          + ρ̅_i-1/2,j+1/2 ( v_i,j+1/2 + v_i-1,j+1/2/2 ) 
            ( u_i-1/2,j + u_i-1/2,j+1/2 ) Δ x_i-1/2,j
          - ρ̅_i-1/2,j-1/2 ( v_i,j-1/2 + v_i-1,j-1/2/2 ) 
            ( u_i-1/2,j + u_i-1/2,j-1/2 ) Δ x_i-1/2,j


    ∮ (μ ▽ u)  dV    = 
         ∑_all facesμ▽ u · ds⃗
       = μ_i,ju_i+1/2,j - u_i-1/2,j/Δ x_i-1/2,jΔ y_i-1/2,j - 
         μ_i-1,ju_i-1/2,j - u_i-3/2,j/Δ x_i-1/2,jΔ y_i-1/2,j
          + μ_i-1/2,j+1/2u_i-1/2,j+1 - u_i-1/2,j/Δ x_i-1/2,jΔ y_i-1/2,j
          - μ_i-1/2,j-1/2u_i-1/2,j - u_i-1/2,j-1/Δ x_i-1/2,jΔ y_i-1/2,j

where Δ x_i-1/2 and Δ y_i-1/2 are the local horizontal and vertical grid spacings. 
ds⃗ is the normalized normal vector of control volume's surface vector and 
ρ̅_i-1/2,j+1/2 denotes the averaged value of density on grid node (i-1/2,j+1/2) 
(similarly for viscosity coefficient). By taking the divergence of equation (<ref>) and enforcing the 
divergence of velocity field at all time steps is zero (equation (<ref>)), we get:

    u⃗^⃗*⃗/Δ t = ( 1/ρ^n+1▽ p^n+1)

This is a Poisson's equation. The Laplacian term is discretized by the second-order centered difference 
scheme. It is then solved by the multigrid method <cit.><cit.>. Finally, 
velocity field at the next time step can be obtained by substituting the updated pressure 
back into equation (<ref>). The same discretization schemes are adopted for the discretization of the Navier-Stokes equations in 
the axisymmetric cylindrical coordinate system. 




 §.§ Immersed boundary method

The basic formulation of the IBM in this work largely follows Unverdi and Tryggvason 
<cit.>. In the IBM, the fluid-fluid interface is represented by a series of 
Lagrangian markers that are allowed to freely move on the fixed Cartesian mesh. To exchange 
information between the interface and the fixed mesh, interpolation and extrapolation are thus 
required. This is done by the discrete version of delta function. The discrete delta function is in 
the same form for both interpolation and extrapolation and it is defined by:

    δ(x⃗ - X⃗) = 1/Δ x Δ yϕ( dx/Δ x) 
                                  ϕ( dy/Δ y)

where X⃗ is the location of a marker, x⃗ is an arbitrary location at which the value of 
delta function is sought, Δ x and Δ y are the local x and y grid spacings 
respectively, dx and dy are x and y components of vector x⃗ - X⃗, and ϕ is a continuous function with compact support. The grid spacing (in both x and y directions) is assumed to be constant in the vicinity of the fluid-fluid interface in this work. The function ϕ has to satisfy several important properties. For example, the sum over all grid points must be equal to one. Further details are given in <cit.>. There are many candidates for ϕ <cit.><cit.>. In this work, we use the first-order Chebyshev polynomial and its expression is:

    ϕ(r) = 
      1/4Δ s( 1 + cosπ r/2Δ s)    if  0 ≤ |r| < 
        2Δ s
    
        0    if  2Δ s ≤ |r|


Using equation (<ref>), the surface tension force on Lagrangian markers can be spread onto the Eulerian grid and velocity field can be interpolated back into the Lagrangian markers. 

As the sharp interface between two fluids is replaced by a smooth differentiable delta function that spreads from two to four grid spacings, the physical properties such as the density and viscosity inside the smooth transition region become statistical. In other words, physical properties no longer follow thermodynamic or physical law in the transition region, in contrast to the diffuse interface method in which the profiles of physical properties are derived based on the second law of thermodynamics <cit.><cit.>. It is assumed that the distributions of the physical properties in the smearing zone follow the shape of smooth Heaviside step function, i.e. 

    ψ = ψ_1 + (ψ_2 - ψ_1)I,

where I is the indicator function and ψ represents a physical variable and the subscript 
denotes fluid 1 or 2 (Figure <ref>). The value of ψ ranges from zero to one, with being zero in fluid 1 and one in fluid 2. The indicator function satisfies the following Poisson's equation:

    ▽I = δ⃗


The immersed boundary forcing term f_ib has to be defined from the jump conditions across the 
interface. Assuming that the surface tension coefficient is constant, the interfacial jump conditions 
with no mass transfer are:

    n̂× (u⃗_2 - u⃗_1) = 0


    ρ_1(u⃗_1 - u⃗_i) ·n̂ = ρ_2(u⃗_2 - u⃗_i) ·n̂ = 0


    p_1 - τ_1 ·n̂ = p_2 - τ_2 ·n̂ - ∇_s σ + σκn̂

where the subscripts 1, 2, and i denote fluid 1, fluid 2, and the interface respectively, 
ρ is the density, p is the pressure, u⃗ is the velocity, τ is the viscous stress, 
∇_S denotes the surface gradient operator, σ is the surface tension constant, κ 
is the curvature, and n̂ is the normalized vector normal to the surface. The first and second equations above simply mean that the velocity is continuous across the interface. 
The third equation can be rewritten as

    [Π]_1,2 = σκn̂

where we have lumped the pressure and viscous stress jumps to the left hand side of the 
equation and denote it by [Π]_1,2. The square bracket means a discontinuity or jump across the 
interface. The right hand side is surface tension force. In the continuous forcing IBM, the 
interface between two fluids is smeared out and thus the physical properties such as density and 
viscosity are assumed to follow Heaviside distribution. 
Therefore, the bracket in equation (<ref>) is to be replaced by the directional 
derivative of Heaviside function and a corresponding delta function is added to the RHS as follows

    ▽Π = σκδn̂


From this form, it is found that the surface tension force can be interpreted as a flux equivalent to 
the combination of pressure flux and viscous stress. Thus it is physically reasonable to add an extra 
surface tension forcing term to the Navier-Stokes equations to enforce the interfacial jump conditions 
between two fluids. However, it is found that the surface tension force becomes rotational even when the curvature and surface tension coefficient is constant if the discrete delta function (<ref>) is substituted into equation (<ref>), leading to a  generation of parasitic currents. 

The other interfacial boundary condition of continuous velocity can be satisfied by using the discrete delta function (eq (<ref>)) to interpolate velocity from the Eulerian grid onto Lagrangian markers because the velocity is defined on markers. 


    U⃗(X⃗) = Σ_i,ju⃗(x⃗) δ(x⃗ - X⃗) Δ x Δ y

where U⃗(X⃗) is the velocity on markers, u⃗(x⃗) is velocity on Eulerian 
grid, and the summation is over the entire computational domain. However in one-fluid finite volume methods, 
the numerical computations are executed in the discrete Cartesian space and each variable is averaged 
over the size of a grid cell. After the independent variables are updated to the next time step, the 
exact location of the interface is lost because the mapping between the Lagrangian markers and 
volume-averaged variables in the Cartesian grid is not one-to-one. Therefore, the location of the 
interface can only be updated approximately. This may lead to numerical error, or noise, in the 
locations of Lagrangian markers. 
Moreover, discretization errors inevitably lead to some small noise in the interpolated velocity on the Lagrangian markers and subsequently the small noise results in small-scale surface fluctuation when the locations of markers are updated. Furthermore, it is well known that a high-order direct polynomial fitting is very sensitive to data points. These factors give rise to the large fluctuation of curvature if smoothing process is not performed. Eventually the fluctuation leads to unphysical pressure gradient and movements of droplets because the surface tension force is a linear function of the curvature. We consider a perturbed circular interface with coordinates of markers (x_i,y_i) given by the following equations to illustrate the sensitivity of the polynomial fitting to the locations of markers. 

    x_i    = 0.5 + 0.5(1 ± g)sin(2π i/n), 
    
    y_i    = 0.5 + 0.5(1 ± g)cos(2π i/n),

where g is the perturbation constant, n is the total number of markers and i is an integer from 0 to n-1. The number of markers is n=126 and they are equally spaced along the interface. The plus or minus sign in eq. (<ref>) is taken according as i is odd or even. Two plots of curvature against the angle are shown in figure <ref> for g=0.001 and g=0.0001 respectively. The amplitude of the fluctuation of curvature is larger than 4 when g=0.001, which corresponds to just 0.1% perturbation of the diameter. Besides numerical errors, addition or deletion of markers can cause fluctuation too. They are a necessary step because two adjacent markers may move too close or too far to each other as time evolves in actual simulations. In the work of <cit.>, when the distance between two adjacent markers is below half of the grid size, one of the two markers is deleted. If the distance is over 1.5 grid size, an additional marker is added at the midpoint of the two markers. To illustrate the fluctuation, we again consider the circular interface in which a marker is added in the midway between two neighboring markers in every other interval. The plot of curvature is shown in figure <ref>. Similar to the previous case, a strong fluctuation can be observed. To eliminate this instability problem, a smooth enough reconstruction of the interface is needed to suppress high frequency fluctuation of curvature. In the next two sections, we will explain how to define a new spreading operator to remove parasitic currents and employ the global B-spline fitting by least squares to enable smooth surface reconstruction in detail. 



 §.§ Non-inertial frame

When simulating a free fall droplet, it may take a long falling distance until it 
reaches its terminal velocity. In order to avoid preparing a large and dense 
computational domain, the non-inertial framework is adopted in the free-fall droplet tests. 

    u⃗t + a⃗  = J⃗

In equation (<ref>), an extra fictitious force a⃗ is added 
to the momentum equation to account for the non-inertial reference frame centered on the droplet. 
Following the work of Komrakova et al <cit.>, the expression for a⃗ is

    a_y = c_1a_y1 + c_2a_y2,


    a_y1 = (COM_o - COM_n)/Δ t,   and


    a_y2 = (COM_n - COM_n-1)/Δ t.

where COM_o is the initial center of mass of the droplet, COM_n and COM_n-1 are center of mass 
at n^th  and n-1^th time steps. The non-inertial force acts only in the vertical component of 
the momentum equation in the direction of the gravitational acceleration. The two constants c_1 and 
c_2 are both set to 0.1. Once the fictitious acceleration is calculated from the equations above, 
it is added as a source term into the Navier-Stokes equations when solving for intermediate velocity 
at every time step. 



 §.§ Conservation of mass

The advection of markers by the discrete delta function does not follow the conservation law on the 
Eulerian grid. Thus, it is not guaranteed that the mass of the dropelt is always conserved. Following 
the work of Udaykumar et al <cit.>, the area enclosed by the markers is calculated every 20 
time steps. The bisection method is then applied if the deviation from the original volume exceeds a 
discrepancy error of order 1 × 10^-6. The discrepancy error is defined as |A_n - A_o|/A_o, where 
A_n is the area at n^th time step and A_o is the initial area. In other words, markers are 
moved inwards or outwards in the normal direction to the interface when the volume growth is 
overestimated or underestimated. 





§ PROBLEM OF PARASITIC CURRENTS



 §.§ Rotational force spreading

The problem of parasitic currents is commonly found in one-fluid multiphase flow methods 
<cit.><cit.>. In the work of <cit.>, they showed by numerical experiments with the VOF 
that the maximum absolute value of fluid velocity caused by parasitic currents is proportional 
to σ/μ. In order to see the cause of parasitic currents, let us define the spreading 
of a constant force, f⃗, by the discrete delta function as a vector function and call it discrete delta vector:

    δ⃗(x⃗ - X⃗) = f⃗δ(x⃗ - X⃗)

where δ is given in (<ref>). A continuous delta vector should be 
irrotational mathematically. This can be easily seen by writing the delta vector in terms of Heaviside function 
and note that curl of gradient is zero. 

    ∇×δ⃗ = ∇×∇ H = 0

In equation (<ref>), the discrete delta vector is splitted into orthogonal directions following the mesh. 
However, this splitting causes the discrete delta vector fail to satisfy the irrotational property, 
resulting in the generation of parasitic currents. In order to demonstrate this problem, let us 
consider a simple linear differential problem in which a straight interface inclined at an angle with 
a constant forcing spread by the delta vector function defined by equation (<ref>). The 
magnitude of force is one and the force is pointing downwards into the lower region. In mathematical 
form: 

    ∇^2 P = ∇·δ⃗

where P is a scalar function that resembles pressure in physcial sense when the force is surface tension force. 
The analytic solution is simply the Heaviside step function obtained by integrating equation 
(<ref>). After substituting equation (<ref>) into (<ref>) and 
solving it numerically by a standard Poisson solver, the results are shown in figure 
<ref>. The grid size is 0.005. The small discrepancy between the exact discrete delta 
forcing and the numerical solution, i.e. δ⃗ - ∇ P, causes the circulation along the 
interface (figure <ref>). 
Figure <ref> further shows a plot of curl for a circular interface. Since the curl of 
the discrete delta vector is not zero in this case again, the solution obtained from solving the 
Poisson's equation will not be a perfect Heaviside function.
The magnitude of curl is constant at a fixed shortest distance to the interface when the grid spacing 
decreases. And it can be shown that the maximum magnitude of curl is at the order of Δ x by 
brute force calculation. Thus, the numerical solution does converge to the analytic solution at first 
order accuracy. Note that this calculation is based on the particular case in which the discrete delta function ϕ is in the cosine form (eq (<ref>)) and there is a possibility that another discrete delta function might not cause the same kind of problem. Therefore, searching for such delta function can be a way to solve this problem. In this work, however, we are taking another approach introducing a new discrete delta vector that satisfies the irrotational condition. Having understood the cause of parasitic currents, we now propose a new discrete delta vector that is irrotational. 



 §.§ Irrotational discrete delta vector


An easy way to ensure that the discrete delta vector is irrotational is that we construct it from a discrete Heaviside function by defining the delta vector as

    δ⃗ = ∇ H

By discretizing the directional derivative, the discrete delta vector can be written as:

    δ⃗_i-1/2,j   = f_i-1/2,jδ_i-1/2,j = f_i-1/2,j(H_i,j - H_i-1,j) 
    δ⃗_i,j-1/2   = f_i,j-1/2δ_i,j-1/2 = f_i,j-1/2(H_i,j - H_i,j-1)

where H_i,j is the discrete Heaviside function defined on the cell's center. f_i-1/2,j is the surface tension and its calculation will be explained in the next section. In this way, the discrete delta vector is guarenteed to be irrotational, that is(f_i,j-1/2δ_i,j-1/2 - f_i-1,j-1/2δ_i-1,j-1/2)Δ x + (f_i-1/2,j-1δ_i-1/2,j-1 - f_i-1/2,jδ_i-1/2,j)Δ y = 0, if the surface tension force is constant. Since the interface is represented by a series of lines and markers, we can easily construct a discrete Heaviside function based on the linear configuration of the interface. As shown in 
figure <ref>, a polygon area bounded by thick line about each segment along the interface is 
constructed. Inside each polygon, the area is divided into five regions I, II, III, IV, and V. The 
magnitude of discrete Heaviside function of each region is computed based on the integration of 
the discrete delta function (equation (<ref>)) along the direction normal to the interface: 

    H    = ∫_-2Δ s^r_i1/4 Δ s( 1 + cosπ r/2Δ s)   dr 
       = 1/4 Δ s( r_i + 2Δ s/π sin ( π r_i/2Δ s + 1/2Δ s) )

Δ s is a constant because the grid size is constant in the vicinity of the interface, r_i is the distance to the interface, and the interface is assumed to be located at r=0. Hence, by substituting values of r_i according to regions I, II, III, IV, and V into equation (<ref>), we obtain

    H = 
      
        0    Region I 
    
        0.409    Region II 
    
        0.5    Region III 
    
        0.909    Region IV 
    
        1    Region V

Suppose we want to find the discrete Heaviside function H_i,j on a particular cell (i,j) as shown in figure <ref>. There are three piece-wise Heaviside's polygons that overlap with cell (i,j). Each overlapped region contributes to H_i,j by the amount of A_I H_I, where A_I is the overlapped area and H_I is the corresponding piece-wise Heaviside function. So H_i,j inside cell (i,j) is then obtained by adding them. 

In the traditional continuous forcing approach, the surface tension force is added to equation 
(<ref>) when the Navier-Stokes equations are solved by the fractional step method, but 
this causes unnecessarily larger smearing of forcing distribution. To demonstrate this, we have 
applied the IBM to incompressible 1D Navier-Stokes equations in cylindrical coordinates. 
We have run two simulations of timesteps 0.005 and 0.0001. The pressure distributions after one time step are shown in figure <ref>. 
The effective range of surface tension force becomes larger as the time step gets larger 
which apparently extends outside of the original range of the discrete delta function. The 
solution will still converge to the true solution in the limit of Δ t → 0. 
The excessive spreading is due to that the immerbsed boundary force is added into equation (<ref>), and thus the range of spreading becomes proportional to the time step. Therefore, we apply the surface tension force in the projection step, i.e. equation (<ref>) instead of (<ref>), in order to avoid this unnecessary larger smearing of the force. 
    ρ^n+1u⃗^* - u⃗^⃗n⃗/Δ t  = RHS,   and


    ρ^n+1u⃗^n+1 - u⃗^*/Δ t  = -▽p^n+1 + f⃗_⃗i⃗b⃗.




§ B-SPLINE RECONSTRUCTION

Surface tension is a linear function of local curvature. Accurate computation of distribution of 
curvature along the interface is thus important. As it has been shown that a direct polynomial interpolation of 
Lagrangian markers is very sensitive to their locations, we propose to use the method of global 
B-spline curve fitting by least squares <cit.>. 
A B-spline curve is controlled by control points and its knot vector. A knot vector is a non-decreasing sequence, denoted by (t_-d,t_-d+1,...,t_g+d,t_g+d+1). Each knot, t_i, represents a data point in the parametric space. The symbol d denotes the degree of interpolating polynomial and g denotes the number of segments. Assume that a fluid-fluid interface is comprised of n connected Lagrangian markers. We define the i^th knot, t_i, as the distance travelled along the interface (which is a set of connected piece-wise line segments) from the first Lagrangian marker to the i^th marker. The first marker is arbitrarily set. Based on the argument above, the B-spline curve is given by:

    B_i,d(x)    = x - t_i/t_i+d - t_iB_i,d-1 + t_i+d+1 - x/t_i+d+1 - t_i+1B_i+1,d-1
    
      B_i,0(x)    = 1

B_i,d is a polynomial of x of degree d and non-zero only in range between knot t_i and knot t_i+d+1. A global interpolating B-spline polynomial B(x,t,d) is 

    B(x,t,d) = Σ_i=-d^gc_i B_i,d(x)

The coefficients c_i are determined by the least squares method. In the least squares method, the following sum of squares is minimized to determine all the coefficients c_i:

    Δ = Σ_r=1^n[ y_r - Σ_i=-d^gc_i B_i,d(x_r) ]^2

The quality of the B-spline curve is often dependent on the arrangement of knots and its degree. 
However, the determination of knots' locations to best-fit the data points is often difficult and 
depends on the nature of the data itself. Thus, throughout this work, we employ a simple strategy to 
determine the knot sequence. The knot sequence is simply set to an arithmetric progression with the 
common distance equal to 1.5Δ s where Δ s is the grid spacing. The order of B-spline 
is 3 to ensure the continuity of curvature. This fitting process by least squares is performed at 
every time step to redistribute the Lagrangian markers so that the distance between two adjacent 
markers is always a constant. The curvature is defined as the rate of change of tangent vector 
along the interface. In discrete form, it can be written as:

    κ_i = |T⃗_i-2Δ s - T⃗_i+2Δ s/Δ s|,

where T⃗_i-2Δ s is the tangent vector whose location is at i^th 
marker and i+1^th marker. A tangent vector at any point on the curve can be obtained by directly differentiating the B-spline curve as follows:

    T⃗ = (B(x,t,d)x,B(y,t,d)y)/|L|

where |L| is the length of the vector (B(x,t,d)x,B(y,t,d)y). In the axisymmetric case, the curvature is given by the following equation, which can be derived by applying the variational principle on the Gibb's free energy, 

    κ = [ ys/x√(1+(ys)^2) + 
    	         ys/√(1+(ys)^2)]


In this work, this step is no longer needed because all markers are redistributed by the smooth global B-spline curve. However, it may be very difficult to directly extend the current method to a closed 2D surface embedded in 3D domain. Instead of global fitting, one may try local fitting by least-squares to reconstruct a smooth interface iteratively. We leave this examination as a future task. 



§ RESULTS

The numerical results are presented to elucidate the differences between new and 
conventional methods in this section. Static droplet and oscillating droplet simulations, which are 
the standard test cases for multiphase flow simulations <cit.><cit.>, have been done to 
demonstrate the validity of the new method. Numerical simulations of free-fall droplet will be given 
in the last part of this section to show its robustness in long-time-scale simulations. 



 §.§ Static droplet simulation


A static droplet simulation is a very simple test case in which a circular or spherical droplet is 
placed in the computational domain without initial motion to verify that the surface tension force is 
properly distributed onto the Cartesian grid, creating a pressure jump across the interface. This test 
has been doen in both 2D and axisymmetric domains. The pressure jump can be obtained analytically 
from the Young-Laplace's equation:

    Δ P = σ( 1/R_1 + 1/R_2)

where R_1 and R_2 are the principal radii of curvature. 

In our 2D static droplet simulations, a droplet of diameter 0.4 is placed in the center of the 
computational domain. This setup is the same for the axisymmetric case, except that the droplet is 
placed at the center of left boundary of the domain (Figure <ref>). All simulations in the 
static droplet tests were run only for one time step of 1 × 10^-3 to observe the pressure 
difference between the ambient fluid and the droplet. The magnitude of surface tension coefficient 
σ is set to one. The exact pressure jump is 5 for a 2D droplet and 10 for an 
axisymmetric droplet. Numerically, we computed the total pressure jump by taking the difference of 
the averages of pressure inside and outside the droplet excluding the interfacial boundary cells:

    Δ P = 1/N_i∑_N_i P_i - 1/N_a∑_N_a P_a

where N_i and N_a are the number of cells of droplet's fluid and ambient fluid, P_i and P_a 
are the pressure inside and outside the interface respectively. Tables <ref> and <ref> give the error of pressure jump 
for different grid sizes. The results are rounded to 5 significant numbers. The density of the ambient fluid is 0.01 and the droplet's density 
is 1. The discrepancy of pressure jump from the correct value is around 0.2% for the conventional scheme which is 
mainly caused by the conventional discrete delta vector (equation (<ref>)) that does not preserve 
the irrotational property. The discrepancy is generally smaller for the new scheme and it is mainly 
due to small fluctuation of the curvature computed from the 
least squares B-spline interpolation. Moreover, the pressure difference converges properly to the analytic value. 
As shown in Figure <ref>, the parasitic currents 
circulating along the interface in the conventional scheme was successfully eliminated when the new 
scheme is used. In conclusion, the results show that the pressure jump of a static droplet is 
accurately reconstructed by the new scheme and its magnitude converge properly to the analytic value 
with increasing grid resolution. This indicates that the interfacial boundary conditions are satisfied 
in the solutions. 







 §.§ Oscillating droplet simulation

The purpose of simulations of an oscillating droplet in this subsection is to test the performance of 
the new scheme when dynamic movement of markers is included. A circular droplet is slightly perturbed 
to an elliptic shape at the beginning, and then it undergoes a free oscillation. The locus of the 
interface is given by:

    x    = x' + r cos(θ)β
    
      y    = y' + r sin(θ)/β

where n is 2, R is 0.2, and r is the droplet's radius. x' and y' are the coordinates 
of the center of droplet. They are (1,1) and (0,1) respectively for the 2D droplet and the 
axisymmetric droplet simulations. The size of the computational domain is 2 by 2. The droplet's 
density and viscosity are 0.01 and 0.01 while the ambient fluid's density and viscosity are 0.01 and 
0.00001. The surface tension coefficient is one. The initial magnitude of perturbation β is 
set to 1.03. Due to unbalanced distribution of pressure jump along the interface, the droplet will 
start to restore back to a circle by pushing inwards on the side of larger curvature and outwards on 
that of smaller curvature. For a small preturbation, the oscillation period can be predicted 
analytically. The period of a 2D oscillating droplet is given by the following 
equation (see <cit.>): 

    T = 2π√((ρ_droplet + ρ_ambient) r^3/(n^3 - n)σ)

where ρ_droplet is droplet's density and ρ_ambient is ambient fluid's density. 
In the axisymmetric case, the period is given by (for derivation, see <cit.>):

    T = 2π√(((n+1)ρ_droplet + n ρ_ambient)r^3/n(n-1)(n+1)(n+2)σ)

Substituting n by 2, the periods are 0.23057 and 0.19935 respectively. 

Three simulations were performed for different grid resolutions to verify the convergence of the new 
scheme. The grid sizes are 0.02, 0.01, and 0.005 respectively. In other words, there are 20, 
40, and 80 cells across the equivalent diameter of the droplet. The first oscillation cycle is 
observed and compared to the analytic solutions. The results for different grid sizes are summarized 
in Table <ref>. The results show that the new scheme does converge properly to the analytic 
solution as resolution increases. The convergence rate is approximatetly at the first order. Note that the 
error is slightly smaller for the conventional scheme. This may be due to the large parasitic 
currents enhancing the frequency by increasing the restoration energy. There is possibly a second reason due to the smoother surface reconstructed by B-spline least squares fitting. The time evolution of the semi minor diameter for simulations of grid size 0.01 using the new scheme is shown in figure <ref>. 

Figure <ref> shows the time series of total kinetic energy (TKE) for simulations with the grid 
size of 0.01. We ran the simulations until time t=8  s and the total kinetic energy of fluid on a Cartesian grid is defined by: 

    TKE = ∑_i,j (u_i,j^2 + v_i,j^2) Δ V_i,j

where Δ V_i,j is the cell's area (or volume in the axisymmetric case). The TKE is diminishing towards zero with the new scheme while it does not decrease properly for the conventional scheme due to the parasitic currents that persistently exist. Figure <ref> further shows the velocity field of a 2D oscillating droplet at time t=2  s. The parasitic currents are clearly visible when using the conventional scheme. Its magnitude is much larger than that of the physical flow field of the oscillating motion. As the parasitic currents are running in tangential direction to the interface, influence on the oscillating motion of the droplet is limited. However, fictituous kinetic energy due to the parasitic currents is persistently generated, and the total energy is not conserved. Finally, we have run a simulation where the droplet is initially deformed with a large amplitude. In this simulation, the β in equation (<ref>) is set to 1.3. The axis ratio and total absolute value of sum of surface tension force spread on the Eulerian grid over sum of absolute value on Lagrangian markers are shown in figure <ref>. The maximum of total force error is at the order of 10^-4 and it is too small to cause significant spurious movements of the droplet even under such a large deformation. 






 §.§ Free-fall droplet

In this section, we apply our new scheme to the simulations of axisymmetric free-fall droplet with low density ratio and 2D free-fall water droplet through still air with large density ratio of 831.67. In the first set of simulations, the droplet's shape and its Reynolds number at terminal velocity is compared to the experimental data for different Morton numbers and Eötvos numbers. In the second set, a circular water droplet of different diameters falls freely until it reaches its terminal velocity. The motivation of these simulations is to show that the new method is capable of dealing with simulations of long time-scale and large density ratio. According to Burkingharm π theorem, there are four dimensionless numbers governing the shape and terminal velocity of a free-fall droplet. These numbers are defined as: 


    d    = ρ_droplet/ρ_ambient, 
    
      v    = μ_droplet/μ_ambient, 
    Eötvos number    = g |ρ_droplet - ρ_ambient|D^2/σ, 
    Morton number    = g μ_ambient^4|ρ_droplet - ρ_ambient|D^2/ρ_ambient^2σ^3

where d and v are density and viscosity ratios respectively. 
Based on experimental data <cit.>, there exists a relationship between the shape regime, Reynolds number, Eötvos number, and Morton number of a free-fall droplet or bubble with small density and viscosity ratios. Readers can refer to figure 2.5 in <cit.> for the relationship. The shapes drawn in their figure are for bubbles. They should be upside-down for our simulations of a free-fall droplet. According to their figure, there are 6 different shape regimes depending on the Morton number and Eötvos number: Shperical, ellipsoidal, dimpled ellipsoidal-cap, skirted, wobbling, and spherical cap. We performed four free-fall droplet simulations with different sets of dimensionless numbers and compare the droplet's shape and Reynolds number with the data on their figure. These dimensionless numbers are tabulated in table <ref>. For the first set of dimensionless numbers, the corresponding shape regime of a droplet at terminal velocity is spherical with no deformation. For the second, third, and fourth sets, the shapes are ellipsoidal, ellipsoidal with larger deformation, and dimpled ellipsoidal-cap respectively. These numbers are so chosen to demonstrate that stable simulations of a free-fall droplet with different magnitudes of deformation, from spherical shape to ellipsoidal-cap shape, at terminal velocity can be achieved by the new method. 

The diameter is non-dimensionalized to one in the simulations. The grid size is 0.0125 for the first two sets of 
parameters and 0.00625 for the other sets. The droplet initially has a perfect spherical shape. The shape of the 
droplet for each set of parameters when it reaches its terminal velocity is shown in figure <ref>. Figure 
<ref> further shows the plot of time evolution of 
Reynolds number against nondimensionial distance that the droplet has travelled. The droplet's shape matches the corresponding shape regime provided in <cit.> and the Reynolds number at terminal velocity is also in good agreement with their experimental data qualitatively. 

In the simulations of 2D free-fall water droplet, we have chosen four different diameters 0.1  mm, 
0.5  mm, 0.8  mm, 1.0  mm to test the performance of the new scheme for high density ratio and 
surface tension force. The grid size is 0.0125 for the simulation of droplet of diameter 1  mm and 0.25 for the others. A plot of falling velocity against distance travelled by the droplet is shown in figure <ref>. The variation in velocity is less than 0.01% at the end of each simulation. This indicates that the numerical solutions converge properly to steady state solutions. From the figure, we also see that droplet with larger diameter travels longer before it reaches its terminal velocity. However, as no experimental 
data exist for free-fall 2D water droplets, no comparisons can be done to verify our numerical 
solutions. In figure <ref>, we plot the velocity field for 0.1  mm water droplet at terminal velocity using the conventional IBM and our new IBM. The grid resolution is 0.0025  mm. The internal circulation is not visible with the conventional method due to the strong parasitic currents. The maximum magnitude of the parasitic currents is as large as 0.2  ms^-1 which is about 40% of the droplet's terminal velocity. Parasitic currents are weaker for a larger droplet because it has a smaller surface tension force. Further application of our new scheme to axisymmetric and 3D free-fall water droplet will be performed for comparisons with the experimental data and published in the near future. 



§ CONCLUSION


In this work, we have shown that preserving the irrotational condition in the discretization of 
the delta function for the force spreading is the key to the elimination of parasitic currents. A new immersed boundary 
method with the discrete delta vector that satisfies the irrotational condition is developed. 
Moreover, we have proposed B-spline curve fitting of the Lagrangian markers by least squares. This enables 
smooth reconstruction of an interface at every time step to filter out unphysical fluctuations of 
curvature due to the spline interpolation that is very sensitive to locations of markers. 
The implementation of this new IBM has been verified through the static droplet and oscillating 
droplet test-cases. The results show that the parasitic currents are eliminated properly and the 
new method maintain a first order accuracy. The application of the irrotational delta vector is not limited to 
the IBM, but can be extended to other one-fluid models such as the VOF and the LS methods. 
Free-fall axisymmetric droplet and 2D water droplet simulations have subsequently been 
performed as an application of the new IBM. The results show good agreement with experimental data, 
demonstrating that the new method is capable of long-time-scale and high density ratio simulations. We hope that the new 
spreading operator that we proposed in this work can be a good alternative to the traditional delta function spreading operator 
when dealing with numerical problems where the internal or surface circulations and mixing are 
important. We are planning to extend it to a 2D interface for general 3D simulations of free-fall 
water droplet in the near future. 



unsrt



§ FIGURES










































