# CLSVOF Method to Study the Formation Process of Taylor Cone in Crater-Like Electrospinning of Nanofibers

Authors: Yong Liu; Jia Li; Yu Tian; Xia Yu; Jian Liu; Bao-Ming Zhou
Journal: Journal of Nanomaterials
 (2014)

---

## Abstract

The application of two-phase computational fluid dynamics (CFD) for simulating crater-like Taylor cone formation dynamics in a viscous liquid is a challenging task. An interface coupled level set/volume-of-fluid (CLSVOF) method and the governing equations based on Navier-Stokes equations were employed to simulate the crater-like Taylor cone formation process. The computational results of the dynamics of crater-like Taylor cone slowly formed on a free liquid surface produced by a submerged nozzle in a viscous liquid were presented in this paper. Some experiments with different air pressures were carried out to evaluate the simulation results. The results from both CFD and experimental observations were compared and analyzed. The numerical results were consistent with the experimental results. Our study showed that the CLSVOF method gave convincing results, and the computational method is robust to extreme variations in interfacial topology.


---

## Body

## 1. Introduction

Electrospinning is a simple, versatile, and effective approach to fabricate nanofibers from various synthetic and natural polymers [1–5]. To date, various nanofibers produced by electrospinning have been applied successfully in many fields, for example, tissue engineering, biotechnology, environmental engineering, filters, and sensors [6–10].A traditional and typical electrospinning setup, which has a thin needle as a spinneret, still has some limitations due to the relatively low yield of the nanofibers. Recently, many novel free-liquid electrospinning techniques, such as needleless electrospinning [11], roller electrospinning [12], wire coil electrospinning [13], porous-tube electrospinning [14], and bubble electrospinning [15, 16], have been invented and have the potential to solve the current problem. Crater-like electrospinning, which creates crater-like solution blowup as a Taylor cone, is also a novel electrospinning process developed on the basis of bubble-electrospinning by our group [17]. This electrospinning process has shown feasibility and potential application for mass production of nanofibers [18].The crater-like electrospinning system consists of a vertical solution reservoir with a gas tube feeding from the reservoir bottom, like a bubble-electrospinning system described in many literatures [15]. One of the differences between the two electrospinning processes is that the crater-like Taylor cone is continuous, but the bubble Taylor cone is not. As shown in Figure 1, multiple jets were ejected from a bubble or crater-like Taylor cone. The former had a process of breaking and producing a new bubble again, while the latter could maintain stable spinning, which is considered as the main reason of a higher yield of nanofibers in this process.Photos of bubble electrospinning (a) and crater-like electrospinning (b).
(a)
(b)In previous studies, the effects of different process parameters on the fabrication of nanofibers in crater-like electrospinning have been investigated [19, 20]. Recently, a multiphysics coupled FEM method was tentatively suggested to simulate the formation of crater-like Taylor Cone [21]. However, there is little literature in which the crater-like solution blowup structure has been studied systematically. The reason might be that the distributed quantities in the fluid field and extreme variations in interfacial topology are difficult to be measured. As a result, the application of two-phase computational fluid dynamics (CFD) for simulating crater-like Taylor cone formation dynamics in a viscous liquid is a challenging task.In the past decades, various numerical methods have been developed to simulate the two-phase flow problem, such as front tracking method [22], marker particle method [23], VOF method [24], and the level set (LS) method [25]. The VOF method forms a building block of computations involving two fluids separated by a sharp interface. The VOF method satisfies compliance with mass conservation extremely well, but sometimes it is difficult to capture the geometric properties of the complicated interface [26]. The LS method was another interface-capturing method, which was first introduced by Osher and Sethian [25]. This method captures the interface very accurately, but in some cases it may violate mass conservation [26]. To achieve mass conservation as well as capture the interface accurately, the CLSVOF method, in which the LS method was coupled with the VOF method, was identified as a better method [27]. In the CLSVOF method, the LS function is used only to compute the geometric properties at the interface, while the void fraction is calculated using the VOF method. Compared with other methods, CLSVOF method can effectively overcome the calculation errors.In this paper, we presented an interface coupled level set/volume-of-fluid (CLSVOF) method for investigating the forming process of crater-like Taylor cone in crater-like electrospinning process. According to the numerical results, the formation mechanism and the formation dynamics of the crater-like Taylor cone in a viscous liquid could be understood for deep explanation of the process.


## 2. Computational Methods

### 2.1. Governing Equations

Generally, for two incompressible and immiscible fluids separated by a moving surface, the Navier-Stokes equations were formulated to describe the motion of both fluids. The whole fluid fills a domain, which may be decomposed into any number of subdomains filled with the individual phases. For the interface between domains, some quantities, such as velocities, are required to be continuous, whereas others, such as pressure, are required to have specific jumps [28]. The discontinuities in physical properties such as density, viscosity, and surface tension at the interface can be treated in different ways; accordingly, the governing equations are reformulated into two useful types: whole-domain formulations and jump-condition formulations [28, 29].


#### 2.1.1. Navier-Stokes Equation: Whole-Domain Formulation

The whole fluid domain can be divided into several subdomains which are occupied by individual phases. At the interface, the discontinuous physical properties, such as density, viscosity, and surface tension, can be smoothed over a transition region of finite thickness. Thus, the whole flow domain can be described by a single set of momentum and continuity equations within the one-field formulation approach, where different fluid properties are considered in each individual phase. Appropriate stress conditions at the interface between different phases are enforced implicitly. With the assumption that the fluid properties are constant in both phases, the mass and momentum conservation equations for the incompressible Newtonian fluids for the liquid and air phases can be written as(1)
∇
·
V
⃑
=
0
,
(2)
ρ
(
∂
V
⃑
∂
t
+
∇
·
(
V
⃑
V
⃑
)
)
=
-
∇
P
+
2
∇
·
μ
s
+
ρ
g
⃑
+
F
b
⃑
,
where V
→ is the velocity, ρ the solution density, p the pressure, g
→ the acceleration of gravity, F
b
→ the body force, μ the dynamic viscosity, and s the strain rate tensor. There are two physical properties of the bulk phases that may have discontinuities across phase boundaries. The strain rate tensor s can be given by
(3)
s
=
1
2
[
(
∇
V
⃑
)
+
(
∇
V
⃑
)
T
]
.The effective densityρ and viscosity μ at each grid point can be expressed as
(4)
ρ
=
ρ
g
(
1
-
H
)
+
ρ
l
H
,
μ
=
μ
g
(
1
-
H
)
+
μ
l
H
,
where subscripts g and l represent gas and liquid, respectively; H is the Heaviside function to prevent numerical instability arising from steep density gradients [29, 30], defined as
(5)
H
=
{
1
ϕ
>
ε
,
0
ϕ
≤
ε
,
1
2
+
ϕ
2
ε
+
1
2
π
[
sin
(
π
ϕ
ε
)
]
|
φ
|
≤
ε
,
where ϕ is the level set function and ε is the interface numerical thickness, which we have taken in our simulations as 1.5
Δ
x (Δ
x refers to the size of a mesh cell). By using the smoothed Heaviside function, one effectively assigns the interface a fixed finite thickness of a small parameter of the order ε, over which the phase properties are interpolated [26].Equation (2) can be discretized as [29]
(6)
ρ
n
V
⃑
n
+
1
-
V
⃑
n
δ
t
=
-
∇
·
(
V
⃑
V
⃑
)
n
ρ
n
-
∇
P
n
+
1
+
∇
·
(
2
μ
s
)
n
+
ρ
n
g
⃑
n
+
F
b
n
⃑
,
where the superscripts n and n
+
1 denote the value of the variable at consecutive time steps and δ is a delta function concentrated at the interface. In the above equation, the pressure P is an implicit term. Other physical quantities, such as gravity, surface tension, and viscosity, are well approximated with t
n values.


#### 2.1.2. Navier-Stokes Equation: Jump-Condition Formulation

The discontinuous quantities, including density, viscosity, and surface tension, at the domain interface are treated in a continuum, where another useful formulation of the momentum-balance equation, jump-condition formulation, can be used.The jump conditions at the interface include the pressure boundary condition written in tensor form as [29]
(7)
p
v
-
p
+
σ
κ
=
-
2
μ
n
k
∂
u
k
∂
n
for the normal direction, where σ is the surface tension coefficient and κ the curvature of the interface. And the velocity boundary condition is
(8)
μ
(
t
i
∂
u
i
∂
n
+
n
k
∂
u
k
∂
s
)
=
∂
σ
∂
s
for the tangential direction, where ∂
/
∂
s and ∂
/
∂
n are the surface and normal derivatives, respectively. For a constant surface tension coefficient, the term on the right-hand side of (8) will be zero.In the fluid domain, the Navier-Stokes equation can be written as(9)
∂
V
→
∂
t
ρ
+
∇
·
(
V
→
V
→
)
ρ
=
-
∇
p
+
2
∇
·
μ
s
+
ρ
g
→
and the continuity equation remains the same as (1). Equation (9) is then discretized as
(10)
V
→
n
+
1
-
V
→
n
δ
t
ρ
n
=
-
∇
·
(
V
→
V
→
)
n
ρ
n
-
∇
p
n
+
1
+
∇
·
(
2
μ
s
)
n
+
ρ
n
g
→
n
.
When proper jump conditions were applied at the interface, (13) could be solved by the same two-step projection method as discussed in [29].


### 2.2. Interface Tracking

The problem in the work is the formation of crater-like solution blowup with large flow distortions and topological changes, for which the Eulerian-based methods are better suited [29]. Two Eulerian-based methods, the volume-of-fluid (VOF) method and the level set (LS) method, have been widely used in this field. The two methods both use phase functions to track the interface implicitly: volume fraction for the VOF method and distance function for the LS method. One of the advantages of the two methods is that they can easily deal with flow problems with large topological changes and interface deformations such as liquid ligament breakup and bubble merging and bursting, which is particularly well suited to this study. However, each method has its own weaknesses. For example, the characteristics of level set determine its flexibility to describe changes in topology and boundary capture, but in the process of solving a serious loss of quality, resulting in decreased accuracy, while VOF lacks accuracy on the normal and curvature calculations. In order to eliminate such weaknesses of the two methods, the coupled level set and volume-of-fluid (CLSVOF) method was employed in this study.


### 2.3. CLSVOF Method

The function of the CLSVOF method is to combine the advantages of the LS method with the VOF method. A flow chart for the CLSVOF algorithm is shown in Figure2. Generally, the interface is reconstructed from the VOF function and the interface normal vector computed from the LS function. On the basis of the above reconstructed interface, the level set functions are redistanced to achieve mass conservation. By combining the advantages of VOF and LS, the CLSVOF method is capable of computing the normal and curvature more accurately while maintaining mass conservation. The coupling of the LS and VOF methods occurs at the interface reconstruction and the redistancing of the level set function.Figure 2
A flow-chart of CLSVOF method.


#### 2.3.1. Volume-of-Fluid Method

In the volume-of-fluid method, the interface is tracked by the VOF function, which is defined as the liquid volume fraction in a cell with its value between zero and one in a surface cell and at zero and one in air and liquid, respectively. One has(11)
F
(
x
⃑
,
t
)
=
{
1
in the fluid
0
<
F
<
1
at
the interface
0
in the void.
An example for the VOF functions representing a delta-shaped region is shown in Figure 3.Figure 3
A sample VOF data on the mesh representing a triangular interface.In Figure3, the number in each cell represents the volume fraction occupied by the liquid. The void fraction F is introduced as the volumetric fraction of the liquid inside a control volume (cell), with the void fraction taking the values 0 for pure gas cell, 1 for pure liquid cell, and between 0 and 1 for a two-phase cell. The VOF functions can be written as
(12)
∂
F
∂
t
+
∇
·
(
V
→
F
)
=
0
.
Equation (12) is rewritten in the conservative form:
(13)
∂
F
∂
t
+
∇
·
(
V
→
F
)
=
F
(
∇
·
V
→
)
.
Equation (13) is discretized temporally and decomposed into two fractional steps [29]:
(14)
F
~
-
F
n
δ
t
+
∂
∂
x
(
u
F
n
)
=
F
~
∂
u
∂
x
,
F
n
+
1
-
F
~
δ
t
+
∂
∂
y
(
v
F
~
)
=
F
n
+
1
∂
v
∂
y
,
where F
~ is the intermediate VOF function. On the staggered grid, the VOF function, F, is located at the cell center and velocities, u and v, are stored at the cell edges, as shown in Figure 4. In Figure 4, the x- and y-velocities are located at the vertical and horizontal cell faces, respectively, and the pressure, the VOF function, and the level set function are stored at the cell centers.Figure 4
Diagram of the discrete variables,u, p, F, and ϕ, in relation to the computational cell.Discretizing the above equations spatially and integrating over a computational cell (i
,
j) yield [29]
(15)
F
~
i
,
j
=
F
i
,
j
n
δ
x
i
δ
y
j
-
δ
t
δ
y
j
(
flu
x
i
+
(
1
/
2
)
,
j
-
flu
x
i
-
(
1
/
2
)
,
j
)
δ
x
i
δ
y
j
-
δ
t
δ
y
j
(
u
i
+
(
1
/
2
)
,
j
-
u
i
-
(
1
/
2
)
,
j
)
,
F
i
,
j
n
+
1
=
F
~
i
,
j
δ
x
i
δ
y
j
-
δ
t
δ
x
i
(
flu
x
i
,
j
+
(
1
/
2
)
-
flu
x
i
,
j
-
(
1
/
2
)
)
δ
x
i
δ
y
j
-
δ
t
δ
x
i
(
v
i
,
j
+
(
1
/
2
)
-
v
i
,
j
-
(
1
/
2
)
)
,
where flu
x
i
±
(
1
/
2
)
,
j
=
(
u
F
n
)
i
±
(
1
/
2
)
,
j and flu
x
i
,
j
±
(
1
/
2
)
=
(
v
F
~
)
i
,
j
±
(
1
/
2
). They denote VOF fluxes across the edges of the computational cell.


#### 2.3.2. Level Set Method

In the LS method, a smooth functionϕ is used to represent a signed distance function whose magnitude equals the shortest distance from the interface. The function ϕ
(
x
→
,
t
), at a point with position vector x
→ and at a time instant t, assumes values as the following:
(16)
ϕ
(
x
→
,
t
)
=
{
>
0
in
the
liquid
region
=
0
at
the
interface
<
0
in
the
gas
region
.
Liquid regions are regions in which ϕ
(
x
→
,
t
)
>
0, while gas regions are regions in which ϕ
(
x
→
,
t
)
<
0. The interface is implicitly represented by the set of points in which ϕ
(
x
→
,
t
)
=
0. One of the advantages of the LS method is its simplicity, especially when computing the curvature κ of the interface. Typically the level set function ϕ
(
x
→
,
t
) is maintained as the signed distance to the interface; that is, ϕ
(
x
→
,
t
)
=
-
d in the gas and ϕ
(
x
→
,
t
)
=
+
d in the liquid, where ϕ
(
x
→
,
t
) is the shortest distance from the point x
→ to the interface at time t. Two-phase cells of the LS function were shown in Figure 5.Figure 5
Two-phase cells of the LS function.The LS function data corresponding to a delta-shaped region are shown in Figure6. All the LS values are located at the cell center and assigned as the shortest distance to the interface. The LS function is initialized as a distance function because of its important property; namely, |
∇
ϕ
|
=
1, which can be used to make a number of simplifications.Figure 6
Level set function values corresponding to a delta-shaped region over a square grid.After initialization, the LS function is moved with the flow field according to the following advection equation:(17)
∂
ϕ
∂
t
+
(
V
→
·
∇
)
ϕ
=
0
.
Since the LS function is smooth and continuous, the discretization of (17) is much more straightforward and some simple advection schemes can be used. However, in order to reduce numerical errors, the level set function must be reinitialized, which can be achieved by obtaining a steady-state solution of the following equation [29]:
(18)
∂
ϕ
∂
t
+
(
w
→
·
∇
)
ϕ
=
ϕ
0
ϕ
0
2
+
h
2
,
where ϕ
0 is the LS function of the previous time step, t the artificial time, h the grid spacing, and w
→ the propagating velocity normal to the interface with unity magnitude, given by
(19)
w
→
=
ϕ
0
ϕ
0
2
+
h
2
(
∇
ϕ
|
∇
ϕ
|
)
.
After the reinitialization process, the level set function will return to a distance function. In order to guarantee mass conservation, the LS functions must be redistanced by calculating the distance from the cell center to the reconstructed interface before being used.


#### 2.3.3. Interface Reconstruction

There are two purposes of the interface reconstruction: one is to calculate the VOF fluxes across each computational cell with an interface and the other is to redistance the LS function for achieving mass conservation [29]. The interface within each cell is approximated by a straight line segment the orientation of which is given by the normal vector. The properly oriented interface is then located in the cell such that the area (volume) is determined from the VOF function.In the CLSVOF method, the interface normal vector and the curvature can be calculated using LS function in all two-phase cells given by(20)
n
→
=
∇
ϕ
|
∇
ϕ
|
=
∇
ϕ
,
κ
=
∇
·
∇
φ
|
∇
φ
|
.
This is different from the usual discontinuous VOF functions. The orientation angle of the interface is then defined as
(21)
α
=
ta
n
-
1
(
n
y
n
x
)
(
0
<
α
≤
2
π
)
,
where α is the angle that the outward pointing unit surface normal makes with the positive x-axis. There are 16 possible cases for the interface shape in the piecewise linear interface construction algorithm. For n
x
>
0 and n
y
>
0, the multitude of possible interface configurations is reduced and there exist 4 cases to be considered, as shown in Figure 7. The line segment is moved along the normal direction to fit the shadow area (volume) with the VOF value in the cell.Figure 7
Four configurations for the interface reconstruction in computational cell.The dashed area (volume) can be calculated by the followingn-sided area (volume) formula:
(22)
A
x
y
=
1
2
∑
i
=
1
n
(
x
i
y
i
+
1
-
x
i
+
1
y
i
)
,
for the two-dimensional case, and
(23)
V
r
z
=
1
6
∑
i
=
1
n
(
r
i
+
r
i
+
1
)
(
r
i
z
i
+
1
-
r
i
+
1
z
i
)
,
for the axisymmetric case. Once the calculated area (volume) matches the VOF value at the cell, the coordinates of endpoints of the line segment are determined and the reconstruction of the interface is completed. Then, the fluxes for the VOF advection can be evaluated based on the reconstructed interface [29]. Details of this procedure can be found in [31, 32].


#### 2.3.4. Reinitialization of the Level Set Function

At each time step after finding the updated LS functionϕ
n
+
1 and the VOF function, F
n
+
1, the LS function must be reinitialized to the exact signed normal distance from the reconstructed interface by coupling the LS function to the volume fraction in order to achieve mass conservation. The reinitialization of the LS function includes initial determination of the sign of the LS function and the subsequent calculation of the shortest distance from the cell centers to the reconstructed interface through a geometric process.For the two-dimensional case, the sign of the LS function,S
ϕ, is given by [33]
(24)
S
ϕ
=
sign
⁡
(
F
-
0.5
)
,
where sign denotes a function that returns the sign of the numeric argument. F is the VOF function for the two-dimensional case.Next, the magnitude of the LS function is determined to find the closest point on an interfacial cell to the neighboring cell centers. Generally, there are two cases for all the interfacial cells. One isF
=
0 or 1 for single-phase cells, and the other is 0
<
F
<
1 for interfacial cells. As shown in Figure 8(a), if the cell (
i
′
,
j
′
) is filled with liquid, the shortest distances are calculated simply by connecting the centers of the neighboring cells to the corners or face centroids of cell (
i
′
,
j
′
). In Figure 8(b), the shortest point on the shadowed area to point A is its projection point onto the line segment within cell (
i
′
,
j
′
) rather than the top right corner at all. For a more general case, as shown in Figure 8(b), for points A and B, the nearest distance is from the cell center to the projection point, like point A in Figure 8(b); for points C and D, the nearest point is the endpoint of the segment, and for the other points, the closest points on cell (
i
′
,
j
′
) are either corners or face centroids [29]. The details of reinitialization of the LS function followed the algorithm presented by [29, 33].Schematic for reinitialization of the LS function.
(a)
(b)


## 2.1. Governing Equations

Generally, for two incompressible and immiscible fluids separated by a moving surface, the Navier-Stokes equations were formulated to describe the motion of both fluids. The whole fluid fills a domain, which may be decomposed into any number of subdomains filled with the individual phases. For the interface between domains, some quantities, such as velocities, are required to be continuous, whereas others, such as pressure, are required to have specific jumps [28]. The discontinuities in physical properties such as density, viscosity, and surface tension at the interface can be treated in different ways; accordingly, the governing equations are reformulated into two useful types: whole-domain formulations and jump-condition formulations [28, 29].


### 2.1.1. Navier-Stokes Equation: Whole-Domain Formulation

The whole fluid domain can be divided into several subdomains which are occupied by individual phases. At the interface, the discontinuous physical properties, such as density, viscosity, and surface tension, can be smoothed over a transition region of finite thickness. Thus, the whole flow domain can be described by a single set of momentum and continuity equations within the one-field formulation approach, where different fluid properties are considered in each individual phase. Appropriate stress conditions at the interface between different phases are enforced implicitly. With the assumption that the fluid properties are constant in both phases, the mass and momentum conservation equations for the incompressible Newtonian fluids for the liquid and air phases can be written as(1)
∇
·
V
⃑
=
0
,
(2)
ρ
(
∂
V
⃑
∂
t
+
∇
·
(
V
⃑
V
⃑
)
)
=
-
∇
P
+
2
∇
·
μ
s
+
ρ
g
⃑
+
F
b
⃑
,
where V
→ is the velocity, ρ the solution density, p the pressure, g
→ the acceleration of gravity, F
b
→ the body force, μ the dynamic viscosity, and s the strain rate tensor. There are two physical properties of the bulk phases that may have discontinuities across phase boundaries. The strain rate tensor s can be given by
(3)
s
=
1
2
[
(
∇
V
⃑
)
+
(
∇
V
⃑
)
T
]
.The effective densityρ and viscosity μ at each grid point can be expressed as
(4)
ρ
=
ρ
g
(
1
-
H
)
+
ρ
l
H
,
μ
=
μ
g
(
1
-
H
)
+
μ
l
H
,
where subscripts g and l represent gas and liquid, respectively; H is the Heaviside function to prevent numerical instability arising from steep density gradients [29, 30], defined as
(5)
H
=
{
1
ϕ
>
ε
,
0
ϕ
≤
ε
,
1
2
+
ϕ
2
ε
+
1
2
π
[
sin
(
π
ϕ
ε
)
]
|
φ
|
≤
ε
,
where ϕ is the level set function and ε is the interface numerical thickness, which we have taken in our simulations as 1.5
Δ
x (Δ
x refers to the size of a mesh cell). By using the smoothed Heaviside function, one effectively assigns the interface a fixed finite thickness of a small parameter of the order ε, over which the phase properties are interpolated [26].Equation (2) can be discretized as [29]
(6)
ρ
n
V
⃑
n
+
1
-
V
⃑
n
δ
t
=
-
∇
·
(
V
⃑
V
⃑
)
n
ρ
n
-
∇
P
n
+
1
+
∇
·
(
2
μ
s
)
n
+
ρ
n
g
⃑
n
+
F
b
n
⃑
,
where the superscripts n and n
+
1 denote the value of the variable at consecutive time steps and δ is a delta function concentrated at the interface. In the above equation, the pressure P is an implicit term. Other physical quantities, such as gravity, surface tension, and viscosity, are well approximated with t
n values.


### 2.1.2. Navier-Stokes Equation: Jump-Condition Formulation

The discontinuous quantities, including density, viscosity, and surface tension, at the domain interface are treated in a continuum, where another useful formulation of the momentum-balance equation, jump-condition formulation, can be used.The jump conditions at the interface include the pressure boundary condition written in tensor form as [29]
(7)
p
v
-
p
+
σ
κ
=
-
2
μ
n
k
∂
u
k
∂
n
for the normal direction, where σ is the surface tension coefficient and κ the curvature of the interface. And the velocity boundary condition is
(8)
μ
(
t
i
∂
u
i
∂
n
+
n
k
∂
u
k
∂
s
)
=
∂
σ
∂
s
for the tangential direction, where ∂
/
∂
s and ∂
/
∂
n are the surface and normal derivatives, respectively. For a constant surface tension coefficient, the term on the right-hand side of (8) will be zero.In the fluid domain, the Navier-Stokes equation can be written as(9)
∂
V
→
∂
t
ρ
+
∇
·
(
V
→
V
→
)
ρ
=
-
∇
p
+
2
∇
·
μ
s
+
ρ
g
→
and the continuity equation remains the same as (1). Equation (9) is then discretized as
(10)
V
→
n
+
1
-
V
→
n
δ
t
ρ
n
=
-
∇
·
(
V
→
V
→
)
n
ρ
n
-
∇
p
n
+
1
+
∇
·
(
2
μ
s
)
n
+
ρ
n
g
→
n
.
When proper jump conditions were applied at the interface, (13) could be solved by the same two-step projection method as discussed in [29].


## 2.1.1. Navier-Stokes Equation: Whole-Domain Formulation

The whole fluid domain can be divided into several subdomains which are occupied by individual phases. At the interface, the discontinuous physical properties, such as density, viscosity, and surface tension, can be smoothed over a transition region of finite thickness. Thus, the whole flow domain can be described by a single set of momentum and continuity equations within the one-field formulation approach, where different fluid properties are considered in each individual phase. Appropriate stress conditions at the interface between different phases are enforced implicitly. With the assumption that the fluid properties are constant in both phases, the mass and momentum conservation equations for the incompressible Newtonian fluids for the liquid and air phases can be written as(1)
∇
·
V
⃑
=
0
,
(2)
ρ
(
∂
V
⃑
∂
t
+
∇
·
(
V
⃑
V
⃑
)
)
=
-
∇
P
+
2
∇
·
μ
s
+
ρ
g
⃑
+
F
b
⃑
,
where V
→ is the velocity, ρ the solution density, p the pressure, g
→ the acceleration of gravity, F
b
→ the body force, μ the dynamic viscosity, and s the strain rate tensor. There are two physical properties of the bulk phases that may have discontinuities across phase boundaries. The strain rate tensor s can be given by
(3)
s
=
1
2
[
(
∇
V
⃑
)
+
(
∇
V
⃑
)
T
]
.The effective densityρ and viscosity μ at each grid point can be expressed as
(4)
ρ
=
ρ
g
(
1
-
H
)
+
ρ
l
H
,
μ
=
μ
g
(
1
-
H
)
+
μ
l
H
,
where subscripts g and l represent gas and liquid, respectively; H is the Heaviside function to prevent numerical instability arising from steep density gradients [29, 30], defined as
(5)
H
=
{
1
ϕ
>
ε
,
0
ϕ
≤
ε
,
1
2
+
ϕ
2
ε
+
1
2
π
[
sin
(
π
ϕ
ε
)
]
|
φ
|
≤
ε
,
where ϕ is the level set function and ε is the interface numerical thickness, which we have taken in our simulations as 1.5
Δ
x (Δ
x refers to the size of a mesh cell). By using the smoothed Heaviside function, one effectively assigns the interface a fixed finite thickness of a small parameter of the order ε, over which the phase properties are interpolated [26].Equation (2) can be discretized as [29]
(6)
ρ
n
V
⃑
n
+
1
-
V
⃑
n
δ
t
=
-
∇
·
(
V
⃑
V
⃑
)
n
ρ
n
-
∇
P
n
+
1
+
∇
·
(
2
μ
s
)
n
+
ρ
n
g
⃑
n
+
F
b
n
⃑
,
where the superscripts n and n
+
1 denote the value of the variable at consecutive time steps and δ is a delta function concentrated at the interface. In the above equation, the pressure P is an implicit term. Other physical quantities, such as gravity, surface tension, and viscosity, are well approximated with t
n values.


## 2.1.2. Navier-Stokes Equation: Jump-Condition Formulation

The discontinuous quantities, including density, viscosity, and surface tension, at the domain interface are treated in a continuum, where another useful formulation of the momentum-balance equation, jump-condition formulation, can be used.The jump conditions at the interface include the pressure boundary condition written in tensor form as [29]
(7)
p
v
-
p
+
σ
κ
=
-
2
μ
n
k
∂
u
k
∂
n
for the normal direction, where σ is the surface tension coefficient and κ the curvature of the interface. And the velocity boundary condition is
(8)
μ
(
t
i
∂
u
i
∂
n
+
n
k
∂
u
k
∂
s
)
=
∂
σ
∂
s
for the tangential direction, where ∂
/
∂
s and ∂
/
∂
n are the surface and normal derivatives, respectively. For a constant surface tension coefficient, the term on the right-hand side of (8) will be zero.In the fluid domain, the Navier-Stokes equation can be written as(9)
∂
V
→
∂
t
ρ
+
∇
·
(
V
→
V
→
)
ρ
=
-
∇
p
+
2
∇
·
μ
s
+
ρ
g
→
and the continuity equation remains the same as (1). Equation (9) is then discretized as
(10)
V
→
n
+
1
-
V
→
n
δ
t
ρ
n
=
-
∇
·
(
V
→
V
→
)
n
ρ
n
-
∇
p
n
+
1
+
∇
·
(
2
μ
s
)
n
+
ρ
n
g
→
n
.
When proper jump conditions were applied at the interface, (13) could be solved by the same two-step projection method as discussed in [29].


## 2.2. Interface Tracking

The problem in the work is the formation of crater-like solution blowup with large flow distortions and topological changes, for which the Eulerian-based methods are better suited [29]. Two Eulerian-based methods, the volume-of-fluid (VOF) method and the level set (LS) method, have been widely used in this field. The two methods both use phase functions to track the interface implicitly: volume fraction for the VOF method and distance function for the LS method. One of the advantages of the two methods is that they can easily deal with flow problems with large topological changes and interface deformations such as liquid ligament breakup and bubble merging and bursting, which is particularly well suited to this study. However, each method has its own weaknesses. For example, the characteristics of level set determine its flexibility to describe changes in topology and boundary capture, but in the process of solving a serious loss of quality, resulting in decreased accuracy, while VOF lacks accuracy on the normal and curvature calculations. In order to eliminate such weaknesses of the two methods, the coupled level set and volume-of-fluid (CLSVOF) method was employed in this study.


## 2.3. CLSVOF Method

The function of the CLSVOF method is to combine the advantages of the LS method with the VOF method. A flow chart for the CLSVOF algorithm is shown in Figure2. Generally, the interface is reconstructed from the VOF function and the interface normal vector computed from the LS function. On the basis of the above reconstructed interface, the level set functions are redistanced to achieve mass conservation. By combining the advantages of VOF and LS, the CLSVOF method is capable of computing the normal and curvature more accurately while maintaining mass conservation. The coupling of the LS and VOF methods occurs at the interface reconstruction and the redistancing of the level set function.Figure 2
A flow-chart of CLSVOF method.


### 2.3.1. Volume-of-Fluid Method

In the volume-of-fluid method, the interface is tracked by the VOF function, which is defined as the liquid volume fraction in a cell with its value between zero and one in a surface cell and at zero and one in air and liquid, respectively. One has(11)
F
(
x
⃑
,
t
)
=
{
1
in the fluid
0
<
F
<
1
at
the interface
0
in the void.
An example for the VOF functions representing a delta-shaped region is shown in Figure 3.Figure 3
A sample VOF data on the mesh representing a triangular interface.In Figure3, the number in each cell represents the volume fraction occupied by the liquid. The void fraction F is introduced as the volumetric fraction of the liquid inside a control volume (cell), with the void fraction taking the values 0 for pure gas cell, 1 for pure liquid cell, and between 0 and 1 for a two-phase cell. The VOF functions can be written as
(12)
∂
F
∂
t
+
∇
·
(
V
→
F
)
=
0
.
Equation (12) is rewritten in the conservative form:
(13)
∂
F
∂
t
+
∇
·
(
V
→
F
)
=
F
(
∇
·
V
→
)
.
Equation (13) is discretized temporally and decomposed into two fractional steps [29]:
(14)
F
~
-
F
n
δ
t
+
∂
∂
x
(
u
F
n
)
=
F
~
∂
u
∂
x
,
F
n
+
1
-
F
~
δ
t
+
∂
∂
y
(
v
F
~
)
=
F
n
+
1
∂
v
∂
y
,
where F
~ is the intermediate VOF function. On the staggered grid, the VOF function, F, is located at the cell center and velocities, u and v, are stored at the cell edges, as shown in Figure 4. In Figure 4, the x- and y-velocities are located at the vertical and horizontal cell faces, respectively, and the pressure, the VOF function, and the level set function are stored at the cell centers.Figure 4
Diagram of the discrete variables,u, p, F, and ϕ, in relation to the computational cell.Discretizing the above equations spatially and integrating over a computational cell (i
,
j) yield [29]
(15)
F
~
i
,
j
=
F
i
,
j
n
δ
x
i
δ
y
j
-
δ
t
δ
y
j
(
flu
x
i
+
(
1
/
2
)
,
j
-
flu
x
i
-
(
1
/
2
)
,
j
)
δ
x
i
δ
y
j
-
δ
t
δ
y
j
(
u
i
+
(
1
/
2
)
,
j
-
u
i
-
(
1
/
2
)
,
j
)
,
F
i
,
j
n
+
1
=
F
~
i
,
j
δ
x
i
δ
y
j
-
δ
t
δ
x
i
(
flu
x
i
,
j
+
(
1
/
2
)
-
flu
x
i
,
j
-
(
1
/
2
)
)
δ
x
i
δ
y
j
-
δ
t
δ
x
i
(
v
i
,
j
+
(
1
/
2
)
-
v
i
,
j
-
(
1
/
2
)
)
,
where flu
x
i
±
(
1
/
2
)
,
j
=
(
u
F
n
)
i
±
(
1
/
2
)
,
j and flu
x
i
,
j
±
(
1
/
2
)
=
(
v
F
~
)
i
,
j
±
(
1
/
2
). They denote VOF fluxes across the edges of the computational cell.


### 2.3.2. Level Set Method

In the LS method, a smooth functionϕ is used to represent a signed distance function whose magnitude equals the shortest distance from the interface. The function ϕ
(
x
→
,
t
), at a point with position vector x
→ and at a time instant t, assumes values as the following:
(16)
ϕ
(
x
→
,
t
)
=
{
>
0
in
the
liquid
region
=
0
at
the
interface
<
0
in
the
gas
region
.
Liquid regions are regions in which ϕ
(
x
→
,
t
)
>
0, while gas regions are regions in which ϕ
(
x
→
,
t
)
<
0. The interface is implicitly represented by the set of points in which ϕ
(
x
→
,
t
)
=
0. One of the advantages of the LS method is its simplicity, especially when computing the curvature κ of the interface. Typically the level set function ϕ
(
x
→
,
t
) is maintained as the signed distance to the interface; that is, ϕ
(
x
→
,
t
)
=
-
d in the gas and ϕ
(
x
→
,
t
)
=
+
d in the liquid, where ϕ
(
x
→
,
t
) is the shortest distance from the point x
→ to the interface at time t. Two-phase cells of the LS function were shown in Figure 5.Figure 5
Two-phase cells of the LS function.The LS function data corresponding to a delta-shaped region are shown in Figure6. All the LS values are located at the cell center and assigned as the shortest distance to the interface. The LS function is initialized as a distance function because of its important property; namely, |
∇
ϕ
|
=
1, which can be used to make a number of simplifications.Figure 6
Level set function values corresponding to a delta-shaped region over a square grid.After initialization, the LS function is moved with the flow field according to the following advection equation:(17)
∂
ϕ
∂
t
+
(
V
→
·
∇
)
ϕ
=
0
.
Since the LS function is smooth and continuous, the discretization of (17) is much more straightforward and some simple advection schemes can be used. However, in order to reduce numerical errors, the level set function must be reinitialized, which can be achieved by obtaining a steady-state solution of the following equation [29]:
(18)
∂
ϕ
∂
t
+
(
w
→
·
∇
)
ϕ
=
ϕ
0
ϕ
0
2
+
h
2
,
where ϕ
0 is the LS function of the previous time step, t the artificial time, h the grid spacing, and w
→ the propagating velocity normal to the interface with unity magnitude, given by
(19)
w
→
=
ϕ
0
ϕ
0
2
+
h
2
(
∇
ϕ
|
∇
ϕ
|
)
.
After the reinitialization process, the level set function will return to a distance function. In order to guarantee mass conservation, the LS functions must be redistanced by calculating the distance from the cell center to the reconstructed interface before being used.


### 2.3.3. Interface Reconstruction

There are two purposes of the interface reconstruction: one is to calculate the VOF fluxes across each computational cell with an interface and the other is to redistance the LS function for achieving mass conservation [29]. The interface within each cell is approximated by a straight line segment the orientation of which is given by the normal vector. The properly oriented interface is then located in the cell such that the area (volume) is determined from the VOF function.In the CLSVOF method, the interface normal vector and the curvature can be calculated using LS function in all two-phase cells given by(20)
n
→
=
∇
ϕ
|
∇
ϕ
|
=
∇
ϕ
,
κ
=
∇
·
∇
φ
|
∇
φ
|
.
This is different from the usual discontinuous VOF functions. The orientation angle of the interface is then defined as
(21)
α
=
ta
n
-
1
(
n
y
n
x
)
(
0
<
α
≤
2
π
)
,
where α is the angle that the outward pointing unit surface normal makes with the positive x-axis. There are 16 possible cases for the interface shape in the piecewise linear interface construction algorithm. For n
x
>
0 and n
y
>
0, the multitude of possible interface configurations is reduced and there exist 4 cases to be considered, as shown in Figure 7. The line segment is moved along the normal direction to fit the shadow area (volume) with the VOF value in the cell.Figure 7
Four configurations for the interface reconstruction in computational cell.The dashed area (volume) can be calculated by the followingn-sided area (volume) formula:
(22)
A
x
y
=
1
2
∑
i
=
1
n
(
x
i
y
i
+
1
-
x
i
+
1
y
i
)
,
for the two-dimensional case, and
(23)
V
r
z
=
1
6
∑
i
=
1
n
(
r
i
+
r
i
+
1
)
(
r
i
z
i
+
1
-
r
i
+
1
z
i
)
,
for the axisymmetric case. Once the calculated area (volume) matches the VOF value at the cell, the coordinates of endpoints of the line segment are determined and the reconstruction of the interface is completed. Then, the fluxes for the VOF advection can be evaluated based on the reconstructed interface [29]. Details of this procedure can be found in [31, 32].


### 2.3.4. Reinitialization of the Level Set Function

At each time step after finding the updated LS functionϕ
n
+
1 and the VOF function, F
n
+
1, the LS function must be reinitialized to the exact signed normal distance from the reconstructed interface by coupling the LS function to the volume fraction in order to achieve mass conservation. The reinitialization of the LS function includes initial determination of the sign of the LS function and the subsequent calculation of the shortest distance from the cell centers to the reconstructed interface through a geometric process.For the two-dimensional case, the sign of the LS function,S
ϕ, is given by [33]
(24)
S
ϕ
=
sign
⁡
(
F
-
0.5
)
,
where sign denotes a function that returns the sign of the numeric argument. F is the VOF function for the two-dimensional case.Next, the magnitude of the LS function is determined to find the closest point on an interfacial cell to the neighboring cell centers. Generally, there are two cases for all the interfacial cells. One isF
=
0 or 1 for single-phase cells, and the other is 0
<
F
<
1 for interfacial cells. As shown in Figure 8(a), if the cell (
i
′
,
j
′
) is filled with liquid, the shortest distances are calculated simply by connecting the centers of the neighboring cells to the corners or face centroids of cell (
i
′
,
j
′
). In Figure 8(b), the shortest point on the shadowed area to point A is its projection point onto the line segment within cell (
i
′
,
j
′
) rather than the top right corner at all. For a more general case, as shown in Figure 8(b), for points A and B, the nearest distance is from the cell center to the projection point, like point A in Figure 8(b); for points C and D, the nearest point is the endpoint of the segment, and for the other points, the closest points on cell (
i
′
,
j
′
) are either corners or face centroids [29]. The details of reinitialization of the LS function followed the algorithm presented by [29, 33].Schematic for reinitialization of the LS function.
(a)
(b)


## 2.3.1. Volume-of-Fluid Method

In the volume-of-fluid method, the interface is tracked by the VOF function, which is defined as the liquid volume fraction in a cell with its value between zero and one in a surface cell and at zero and one in air and liquid, respectively. One has(11)
F
(
x
⃑
,
t
)
=
{
1
in the fluid
0
<
F
<
1
at
the interface
0
in the void.
An example for the VOF functions representing a delta-shaped region is shown in Figure 3.Figure 3
A sample VOF data on the mesh representing a triangular interface.In Figure3, the number in each cell represents the volume fraction occupied by the liquid. The void fraction F is introduced as the volumetric fraction of the liquid inside a control volume (cell), with the void fraction taking the values 0 for pure gas cell, 1 for pure liquid cell, and between 0 and 1 for a two-phase cell. The VOF functions can be written as
(12)
∂
F
∂
t
+
∇
·
(
V
→
F
)
=
0
.
Equation (12) is rewritten in the conservative form:
(13)
∂
F
∂
t
+
∇
·
(
V
→
F
)
=
F
(
∇
·
V
→
)
.
Equation (13) is discretized temporally and decomposed into two fractional steps [29]:
(14)
F
~
-
F
n
δ
t
+
∂
∂
x
(
u
F
n
)
=
F
~
∂
u
∂
x
,
F
n
+
1
-
F
~
δ
t
+
∂
∂
y
(
v
F
~
)
=
F
n
+
1
∂
v
∂
y
,
where F
~ is the intermediate VOF function. On the staggered grid, the VOF function, F, is located at the cell center and velocities, u and v, are stored at the cell edges, as shown in Figure 4. In Figure 4, the x- and y-velocities are located at the vertical and horizontal cell faces, respectively, and the pressure, the VOF function, and the level set function are stored at the cell centers.Figure 4
Diagram of the discrete variables,u, p, F, and ϕ, in relation to the computational cell.Discretizing the above equations spatially and integrating over a computational cell (i
,
j) yield [29]
(15)
F
~
i
,
j
=
F
i
,
j
n
δ
x
i
δ
y
j
-
δ
t
δ
y
j
(
flu
x
i
+
(
1
/
2
)
,
j
-
flu
x
i
-
(
1
/
2
)
,
j
)
δ
x
i
δ
y
j
-
δ
t
δ
y
j
(
u
i
+
(
1
/
2
)
,
j
-
u
i
-
(
1
/
2
)
,
j
)
,
F
i
,
j
n
+
1
=
F
~
i
,
j
δ
x
i
δ
y
j
-
δ
t
δ
x
i
(
flu
x
i
,
j
+
(
1
/
2
)
-
flu
x
i
,
j
-
(
1
/
2
)
)
δ
x
i
δ
y
j
-
δ
t
δ
x
i
(
v
i
,
j
+
(
1
/
2
)
-
v
i
,
j
-
(
1
/
2
)
)
,
where flu
x
i
±
(
1
/
2
)
,
j
=
(
u
F
n
)
i
±
(
1
/
2
)
,
j and flu
x
i
,
j
±
(
1
/
2
)
=
(
v
F
~
)
i
,
j
±
(
1
/
2
). They denote VOF fluxes across the edges of the computational cell.


## 2.3.2. Level Set Method

In the LS method, a smooth functionϕ is used to represent a signed distance function whose magnitude equals the shortest distance from the interface. The function ϕ
(
x
→
,
t
), at a point with position vector x
→ and at a time instant t, assumes values as the following:
(16)
ϕ
(
x
→
,
t
)
=
{
>
0
in
the
liquid
region
=
0
at
the
interface
<
0
in
the
gas
region
.
Liquid regions are regions in which ϕ
(
x
→
,
t
)
>
0, while gas regions are regions in which ϕ
(
x
→
,
t
)
<
0. The interface is implicitly represented by the set of points in which ϕ
(
x
→
,
t
)
=
0. One of the advantages of the LS method is its simplicity, especially when computing the curvature κ of the interface. Typically the level set function ϕ
(
x
→
,
t
) is maintained as the signed distance to the interface; that is, ϕ
(
x
→
,
t
)
=
-
d in the gas and ϕ
(
x
→
,
t
)
=
+
d in the liquid, where ϕ
(
x
→
,
t
) is the shortest distance from the point x
→ to the interface at time t. Two-phase cells of the LS function were shown in Figure 5.Figure 5
Two-phase cells of the LS function.The LS function data corresponding to a delta-shaped region are shown in Figure6. All the LS values are located at the cell center and assigned as the shortest distance to the interface. The LS function is initialized as a distance function because of its important property; namely, |
∇
ϕ
|
=
1, which can be used to make a number of simplifications.Figure 6
Level set function values corresponding to a delta-shaped region over a square grid.After initialization, the LS function is moved with the flow field according to the following advection equation:(17)
∂
ϕ
∂
t
+
(
V
→
·
∇
)
ϕ
=
0
.
Since the LS function is smooth and continuous, the discretization of (17) is much more straightforward and some simple advection schemes can be used. However, in order to reduce numerical errors, the level set function must be reinitialized, which can be achieved by obtaining a steady-state solution of the following equation [29]:
(18)
∂
ϕ
∂
t
+
(
w
→
·
∇
)
ϕ
=
ϕ
0
ϕ
0
2
+
h
2
,
where ϕ
0 is the LS function of the previous time step, t the artificial time, h the grid spacing, and w
→ the propagating velocity normal to the interface with unity magnitude, given by
(19)
w
→
=
ϕ
0
ϕ
0
2
+
h
2
(
∇
ϕ
|
∇
ϕ
|
)
.
After the reinitialization process, the level set function will return to a distance function. In order to guarantee mass conservation, the LS functions must be redistanced by calculating the distance from the cell center to the reconstructed interface before being used.


## 2.3.3. Interface Reconstruction

There are two purposes of the interface reconstruction: one is to calculate the VOF fluxes across each computational cell with an interface and the other is to redistance the LS function for achieving mass conservation [29]. The interface within each cell is approximated by a straight line segment the orientation of which is given by the normal vector. The properly oriented interface is then located in the cell such that the area (volume) is determined from the VOF function.In the CLSVOF method, the interface normal vector and the curvature can be calculated using LS function in all two-phase cells given by(20)
n
→
=
∇
ϕ
|
∇
ϕ
|
=
∇
ϕ
,
κ
=
∇
·
∇
φ
|
∇
φ
|
.
This is different from the usual discontinuous VOF functions. The orientation angle of the interface is then defined as
(21)
α
=
ta
n
-
1
(
n
y
n
x
)
(
0
<
α
≤
2
π
)
,
where α is the angle that the outward pointing unit surface normal makes with the positive x-axis. There are 16 possible cases for the interface shape in the piecewise linear interface construction algorithm. For n
x
>
0 and n
y
>
0, the multitude of possible interface configurations is reduced and there exist 4 cases to be considered, as shown in Figure 7. The line segment is moved along the normal direction to fit the shadow area (volume) with the VOF value in the cell.Figure 7
Four configurations for the interface reconstruction in computational cell.The dashed area (volume) can be calculated by the followingn-sided area (volume) formula:
(22)
A
x
y
=
1
2
∑
i
=
1
n
(
x
i
y
i
+
1
-
x
i
+
1
y
i
)
,
for the two-dimensional case, and
(23)
V
r
z
=
1
6
∑
i
=
1
n
(
r
i
+
r
i
+
1
)
(
r
i
z
i
+
1
-
r
i
+
1
z
i
)
,
for the axisymmetric case. Once the calculated area (volume) matches the VOF value at the cell, the coordinates of endpoints of the line segment are determined and the reconstruction of the interface is completed. Then, the fluxes for the VOF advection can be evaluated based on the reconstructed interface [29]. Details of this procedure can be found in [31, 32].


## 2.3.4. Reinitialization of the Level Set Function

At each time step after finding the updated LS functionϕ
n
+
1 and the VOF function, F
n
+
1, the LS function must be reinitialized to the exact signed normal distance from the reconstructed interface by coupling the LS function to the volume fraction in order to achieve mass conservation. The reinitialization of the LS function includes initial determination of the sign of the LS function and the subsequent calculation of the shortest distance from the cell centers to the reconstructed interface through a geometric process.For the two-dimensional case, the sign of the LS function,S
ϕ, is given by [33]
(24)
S
ϕ
=
sign
⁡
(
F
-
0.5
)
,
where sign denotes a function that returns the sign of the numeric argument. F is the VOF function for the two-dimensional case.Next, the magnitude of the LS function is determined to find the closest point on an interfacial cell to the neighboring cell centers. Generally, there are two cases for all the interfacial cells. One isF
=
0 or 1 for single-phase cells, and the other is 0
<
F
<
1 for interfacial cells. As shown in Figure 8(a), if the cell (
i
′
,
j
′
) is filled with liquid, the shortest distances are calculated simply by connecting the centers of the neighboring cells to the corners or face centroids of cell (
i
′
,
j
′
). In Figure 8(b), the shortest point on the shadowed area to point A is its projection point onto the line segment within cell (
i
′
,
j
′
) rather than the top right corner at all. For a more general case, as shown in Figure 8(b), for points A and B, the nearest distance is from the cell center to the projection point, like point A in Figure 8(b); for points C and D, the nearest point is the endpoint of the segment, and for the other points, the closest points on cell (
i
′
,
j
′
) are either corners or face centroids [29]. The details of reinitialization of the LS function followed the algorithm presented by [29, 33].Schematic for reinitialization of the LS function.
(a)
(b)


## 3. Computational Domain and Boundary Conditions

According to the characteristic of crater-like Taylor cone formation, a two-dimensional axisymmetric geometry model was established. Figure9 illustrates the schematic view of the flow domain used in the two-dimensional simulation. In this case, an 18 wt% polyvinyl alcohol (PVA, M
-
w
=
30000 g/mol)/distilled water solution was put into a custom-made quartz circular cylinder chamber with diameter D
=
40.0 mm and height H
=
35 mm. A gas tube with internal diameter d
=
2 mm was mounted in the center bottom of the chamber. The distance was 100 mm between the collecting electrode and the surface of solution. The simulations were performed using the following parameters: the gas is set as the standard values of the air, the mass density of PVA ρ
=
1100 Kg/m3, the viscosity of PVA solution μ
=
0.87 Pa·S, the surface tension σ
=
0.04 N/m, and the gas pressure P
=
4
–
50 KPa. In our experiments, the liquid surface is higher than the jet for 4 mm at the initial state. The pressure on the other free surface is set as 0 Pa. No-slip boundary conditions are used at all walls. The direction of gravity is along the vertical direction as shown in Figure 9.Figure 9
Illustration of the experimental zone.The pressure-implicit with splitting of operators (PISO) pressure-velocity coupling scheme was used to calculate the transient two-phase fluid flow, and the iterative time step is set to 10−6 S to ensure the accurate results of crater formation.


## 4. Results and Discussions

### 4.1. The Formation of Crater-Like Taylor Cone

When the inlet pressure is set as 4 kPa, the transient simulated results of 0–216 ms are shown in Figure10. At the gas pressure, the simulated shape of Taylor cone was like a bubble. With the increase of the gas pressure, the shape of Taylor cone changed. When the gas pressure was over 10 kPa, the frequency of bubble occurrence decreased. The transient shapes of Taylor cone were much more crater-like blowup. The above simulated results were consistent with the experimental results. In the simulated results, there are many tiny bubbles or drops which might be stretched into the charged jets in the experiments.Figure 10
The simulated results of crater-like Taylor cone formation at the gas pressure 4 kPa.


### 4.2. The Effect of Gas Pressure on the Formation Period of Crater-Like Taylor Cone

Both our experimental results and numerical simulation results showed that there existed a formation period of crater-like Taylor cone. The simulated results also showed that the crater-like Taylor cone formation period also changed with the increase of gas pressure. In this case, the Taylor cone, which underwent a life-cycle process of growing and bursting, is not steady. When the gas pressure was 4 kPa, there was a big bubble, as Taylor cone, on the solution surface and the life-cycle period of Taylor cone was about 216 ms. When the gas pressure was 10 kPa, the formation period dramatically reduced to about 68 ms. And the occurrence of bubble will be significantly reduced, as shown in Figure11.Figure 11
The simulated results of crater-like Taylor cone formation at the gas pressure 10 kPa.Along with the continued increase of gas pressure, the crater-like Taylor cone formation period decreased. Additionally, the occurrence of bubble also decreased with increasing gas pressure. For example, the simulated results of 25 kPa were shown in Figure12. The life-cycle period of crater-like Taylor cone formation was 52 ms and bubble did not appear during the period, which showed that the crater-like Taylor cone formed completely. When the air pressure reached 16, 35, and 50 kPa, the formation periods obtained from the simulated results were 56, 40, and 28 ms, respectively.Figure 12
The formation of crater while the inlet pressure is 25 kPa.The reason for the decrease of the formation period might be that there was more amount of gas passed through the solution at the same period of time that the gas pressure was increased. The experimental results agree well with the simulated results, as shown in Figure13, which indicates that the numerical model is validated to study the formation of crater-like Taylor cone.Figure 13
The relationship between gas pressure and the formation period.


### 4.3. The Effect of Gas Pressure on the Height of Crater-Like Taylor Cone

Besides the life-cycle period, the effect of gas pressure on the height of crater-like Taylor cone was also studied. Here the height of crater-like Taylor cone represents the maximum height from the solution surface to the tip of the blowup. The solution’s depth between the solution surface and the tip of gas tube was kept constant, and the applied voltage was 0 kV. The simulations were performed at different gas pressures. According to the simulated results, the height of crater-like Taylor cone increased with the increase of gas pressure, which agreed well with the experimental results. However, the relationship between the gas pressure and the height of Taylor cone did not remain constant along the increase of gas pressure. When the gas pressure was over 30 kPa, the growth rate of height has decreased or stopped, as shown in Figure14.Figure 14
The effect of gas pressure on the height of crater-like Taylor cone.


### 4.4. The Effect of Solution Depth on the Height of Crater-Like Taylor Cone

In the process of crater-like Taylor cone electrospinning, the solution depth between the solution surface and the tip of gas tube is another important factor that affects the height of Taylor cone. According to our previous experiments, the morphologies of nanofibers can also be influenced by the solution depth. In order to further verify the validation of the numerical method, the height of crater-like Taylor cone was calculated at different gas pressure. In our simulations, the inlet pressure was fixed at 25 kPa and the applied voltage 0 kV. The numerical results showed that the height of crater-like Taylor cone decreased with the increase of solution depth, which is in keeping with the experimental results, as shown in Figure15. The reason for this might be that the higher the solution depth, the greater the weight of the overlying fluid, which requires more gas pressure to form the crater-like Taylor cone.Figure 15
The effect of solution depth on the height of crater-like Taylor cone.


## 4.1. The Formation of Crater-Like Taylor Cone

When the inlet pressure is set as 4 kPa, the transient simulated results of 0–216 ms are shown in Figure10. At the gas pressure, the simulated shape of Taylor cone was like a bubble. With the increase of the gas pressure, the shape of Taylor cone changed. When the gas pressure was over 10 kPa, the frequency of bubble occurrence decreased. The transient shapes of Taylor cone were much more crater-like blowup. The above simulated results were consistent with the experimental results. In the simulated results, there are many tiny bubbles or drops which might be stretched into the charged jets in the experiments.Figure 10
The simulated results of crater-like Taylor cone formation at the gas pressure 4 kPa.


## 4.2. The Effect of Gas Pressure on the Formation Period of Crater-Like Taylor Cone

Both our experimental results and numerical simulation results showed that there existed a formation period of crater-like Taylor cone. The simulated results also showed that the crater-like Taylor cone formation period also changed with the increase of gas pressure. In this case, the Taylor cone, which underwent a life-cycle process of growing and bursting, is not steady. When the gas pressure was 4 kPa, there was a big bubble, as Taylor cone, on the solution surface and the life-cycle period of Taylor cone was about 216 ms. When the gas pressure was 10 kPa, the formation period dramatically reduced to about 68 ms. And the occurrence of bubble will be significantly reduced, as shown in Figure11.Figure 11
The simulated results of crater-like Taylor cone formation at the gas pressure 10 kPa.Along with the continued increase of gas pressure, the crater-like Taylor cone formation period decreased. Additionally, the occurrence of bubble also decreased with increasing gas pressure. For example, the simulated results of 25 kPa were shown in Figure12. The life-cycle period of crater-like Taylor cone formation was 52 ms and bubble did not appear during the period, which showed that the crater-like Taylor cone formed completely. When the air pressure reached 16, 35, and 50 kPa, the formation periods obtained from the simulated results were 56, 40, and 28 ms, respectively.Figure 12
The formation of crater while the inlet pressure is 25 kPa.The reason for the decrease of the formation period might be that there was more amount of gas passed through the solution at the same period of time that the gas pressure was increased. The experimental results agree well with the simulated results, as shown in Figure13, which indicates that the numerical model is validated to study the formation of crater-like Taylor cone.Figure 13
The relationship between gas pressure and the formation period.


## 4.3. The Effect of Gas Pressure on the Height of Crater-Like Taylor Cone

Besides the life-cycle period, the effect of gas pressure on the height of crater-like Taylor cone was also studied. Here the height of crater-like Taylor cone represents the maximum height from the solution surface to the tip of the blowup. The solution’s depth between the solution surface and the tip of gas tube was kept constant, and the applied voltage was 0 kV. The simulations were performed at different gas pressures. According to the simulated results, the height of crater-like Taylor cone increased with the increase of gas pressure, which agreed well with the experimental results. However, the relationship between the gas pressure and the height of Taylor cone did not remain constant along the increase of gas pressure. When the gas pressure was over 30 kPa, the growth rate of height has decreased or stopped, as shown in Figure14.Figure 14
The effect of gas pressure on the height of crater-like Taylor cone.


## 4.4. The Effect of Solution Depth on the Height of Crater-Like Taylor Cone

In the process of crater-like Taylor cone electrospinning, the solution depth between the solution surface and the tip of gas tube is another important factor that affects the height of Taylor cone. According to our previous experiments, the morphologies of nanofibers can also be influenced by the solution depth. In order to further verify the validation of the numerical method, the height of crater-like Taylor cone was calculated at different gas pressure. In our simulations, the inlet pressure was fixed at 25 kPa and the applied voltage 0 kV. The numerical results showed that the height of crater-like Taylor cone decreased with the increase of solution depth, which is in keeping with the experimental results, as shown in Figure15. The reason for this might be that the higher the solution depth, the greater the weight of the overlying fluid, which requires more gas pressure to form the crater-like Taylor cone.Figure 15
The effect of solution depth on the height of crater-like Taylor cone.


## 5. Conclusion

In this paper we suggested a numerical approach, CLSVOF method, to numerically simulate the formation of crater-like Taylor cone electrospinning process. Numerical simulation was performed for two-dimensional uncompressed flow in axisymmetic coordinates. The numerical results showed that the formation period of crater-like Taylor cone decreased with the increasing gas pressure. The height of crater-like Taylor cone increased with the increase of gas pressure, but the height decreased as the solution depth increased. The numerical results are consistent with the experimental results. The numerical method could be helpful to understand the mechanism of electrospinning process and improve the production rate of nanofibers.



---

*Source: 635609-2014-06-11.xml*